0 1

将Simulink模型移植到C#

最近工作中需要将一个比较复杂的Simulink模型用C#实现,着实花费了一番功夫,下面将移植的方法总结一下。

 模型的离散化

Simulink的模型可以是连续的,用复杂的ode Solver求解,Solver会自动调整每次运算的时间间隔,以达到最精确的结果。这样的模型很难移植,因此为了程序编写方便,先将Simulink的Solver改为Discrete(离散),并且设定为Fixed-Step,这样离散化了的模型就容易实现多了。

如果模型中包括一些Continuous(连续)的库模块的话,需要把它们都手工转换为响应的离散的模块,例如:



转换之后的模型的输出和原始的模型的输出可能会出现比较大的差别,这就需要调整离散模型的step size,以及各个模块的参数,直到连续模型和离散模型的输出一致。

我们可以分别将两种模型的输出保存到Workspace中的变量中,然后进行一致性比较。因为两种模型的时间间隔是不同的,因此不能直接进行数值上的比较,而需要将连续模型的输出按照离散模型的时间重新取样(简单的插值运算即可)。

在实际转换中我发现连续模型的Transfer Fcn(传递函数)模块比较难找到对应的离散模型,需要花费一番功夫。

模型的转化也可以试试Tool->Control Design->Model Discretizer,让Simulink自动进行离散转化,我的模型不知道为什么转化失败,所以只好进行手工转化了。

各个模块的输入输出

下面的工作就是将各个模块移植到C#中,首先通过Sinks中的To Workspace模块,将被移植模块的输入输出全部保存到一个变量中。如下图蓝色模块所示,将Velocity的输入输出都保存到simout变量中,并且SaveFormat设置为Array。


运行模拟之后,运行下面的命令将simout中的结果保存到test.csv文件中,为了和C#程序的输出进行比较,我将输出精度设为最大,保证没有任何进度的损失:

dlmwrite('test.csv', simout, 'delimiter', ',', 'precision', 20)

C#程序的编写

剩下的事情就是编写一个程序让其输入输出和test.csv中的输入输出一致。因为可能需要对程序进行多次调试,和进行多组输入输出的比对,所以最好用单元测试软件NUnit将这些测试工作自动化。关于这一节,在今后的文章中详述。

没有延时效果的模块可以用一个函数对其输入输出进行描述:

doube[] Do(double[] inputs);

也就是说,得到此模块的输入之后,将输入传递给此模块,就得到其输出了。

但是例如积分器和延时器之类输入输出之间有时间差的模块就比较麻烦了。这样的模块如果其输出没有连接到输入的话(形成回路),那么可以先计算其输入,然后计算其输出,所以也可以用上面的函数。如果输出连接到自己的输入话,那么要先计算输出,计算模块的输出时要使用模块的即时输入和模块内部的状态,然后通过输出再计算输入,并将输入的值传递给模块,更新模块的状态。例如要计算如下的模型:


这个模型中的积分器的输出又回到了自己的输入,因此形成了回路。这时需要先计算积分器的输出值p1,然后将此值加1得到p2: p2 = p1 + 1,然后将p2输入到积分器,积分器对p2进行累加,以得到积分器的下一个输出。所以这个模型的输出为:

0, 1, 3, 7, 15, 31, 63, 127, ...

关于这方面的内容,今后再写文章详述。


TiddlyWiki的抽认卡插件

关于抽认卡

一个抽认卡(或称闪卡(英文:Flash Card))是一小块纸片,在学校里用来作辅助教学(主要是英语国家)。抽认卡可以用来记录词汇、历史事件时间、公式等等。使用抽认卡的目的主要是帮助记忆,提高学生的学习动机、以及在课堂上为学生提供多种感官的刺激。你在每个卡片上写下一个问题(同时在背面记下答案),用它们来测试自己,并根据你的测试、学习结果把它们进行排序、分组。

这种策略使得学习更有选择性,也就是:在一组的卡片越难,你就越要花时间复习这一组。于是,你就能节省很多学习时间。

※ 摘自维基百科

抽认卡程序

http://memorizable.com/ 就是采用类似抽认卡的方法帮助用户记忆的网站。他提供了javascript的源程序。

最近发现TiddlyWiki版本更新为2.5,并且集成了jQuery库,这两个我最喜欢的javascript程序最终结合到了一起,真令人兴奋。于是我将其稍作修改,写成了一个TiddlyWiki的插件,并且专门做了一个TiddlySite来演示其效果: http://hyry.dip.jp/kt/index.py?site=memory,关于插件的源程序可以在 {CodeLib} 中找到。同时TiddlySite也更新为2.5,调试正常之后将发布最新的代码。

这个抽认卡程序效果很酷,但是写得比较乱,如果有时间能用jQuery重写一个就好了。

 


超声波探伤系统

前两个月一直在忙超声波探伤系统的开发。它是为某个圆柱型钢材加工工厂的质量检测部门开发的系统。软件部分需要控制各种机械部件、电机、传感器,以完成自动超声波探伤。整个系统的构成相当复杂,程序需要从近30多个传感器获取数据,控制6个电机的运动,并且实时地从两块高速数据采集卡PCI-9812获取超声波数据,判断钢材中是否有缺陷,并绘制缺陷示意图。

在研究室里写完整个程序的基本结构,花了近一个月的时间。然后就出差到工厂里面实地调试了两个星期。在实际调试的时候才发现有许多东西并不是按照想像的那样运行,例如传感器可能出现误报;电机的Encoder信号可能出现噪声,以至于无法正确判断电机的位置;激光测距仪可能因为材料表面的反射条件而无法正确工作等等。而最令人头疼的是,工厂环境下的电气噪声和接触媒质(水)的不稳定会影响超声波探头的信号,使得缺陷反射的超声波信号和噪声混杂在一起,无法进行正确判断。看着搞装置的同事们想尽各种办法解决这些问题,我觉得我的实践经验值也突飞猛涨。将来有时间要好好写写这个系统开发过程中的经验教训。

客户对缺陷示意图的要求只要是2D的平面图就可以了。不过由于钢材是圆柱型的,平面图很难直观地表现缺陷的位置。于是我利用空闲时间用Python写了一个3D表现缺陷的小程序。下面是这个程序的界面截图,上方显示的是缺陷的平面图。可以看出平面图中蓝色区域标示了许多缺陷(红色的点),它们的实际位置在钢材的圆心部分,实际上是钢材中心的孔。

这个程序采用wxPython做界面,pyOpengl绘制3D图形,用wxPython的Plot模块绘制AScope(超声波的反射信号),用户可以用鼠标任意旋转、移动、放大缺陷,选中某个缺陷之后,可以详细查看其AScope波形,应该能为缺陷的人工检定带来不少的便利。



NetLogo - 草、羊、狼的模型

NetLogo 是一个很奇妙的软件,它采用logo语言的语法,能够很方便地控制许许多多的海龟在一个世界中运动并交互,以此来观察个体和宏观的行为。用NetLogo自己的话来说,NetLogo是一个用来模拟自然现象和社会现象的可编程建模环境。它很适合为随时间而发展的复杂系统建立模型。它可以同时给成百上千的"行动者"(海龟)发布指令,分别对其进行操作。这样使我们能够很清楚地观察和研究微观级别的个体行为是如何产生宏观的行为模式。

NetLogo自己带了200多个示例,涵盖生物、物理、化学、社会、数学等多个领域。随便打开一个研究研究,就能不知不觉地度过几个小时。通过观察各种各样的参数对模型的行为影响,能使我们更深刻地理解它所要表达的内容。并且它用来建模的logo语言十分简单,通常一个模型的源程序都不超过100行。因此你可以很快地阅读完模型的源程序,理解其运作的原理并加以修改。

让我们打开一个叫做Wolf Sheep Predation 的示例模型看看吧。这是一个简单的模拟食物链的模型。 选择菜单File->Models Library,然后在对话框中找到Biology->Wolf Sheep Predation。

这个模型中有一种植物:草;两种动物:羊和狼。羊和狼都在草原上游荡。羊和狼每走一步都需要耗费一定的能量,而当其能量降低到0时,它就死亡了。羊每次行动之后,都试图吃掉它所在地方的草,以获取能量维持生命;同样,狼也会试图吃掉它所在地方的羊。羊和狼都按照一定的概率生产新的羊和狼。我们希望观察在这个模型中,随着时间的流逝羊和狼的数目会如何变化。

这个模型有两个版本,第一个是草资源无限丰富,羊不会被饿死;第二个版本中,羊吃掉附近的草之后,那块地方的草资源就为0了,必须等到一定时间才能再长出草来,因此羊会因为无草吃而饿死。

打开这个模型之后,我们会看到在Interface选项卡中有很多控件:按钮、开关、滑动条以及统计图。我们通过这些控件来和模型进行交互。

让我们来看看这些控件的具体含义。

  • Initial-number-sheep: 草原上的羊的初始数目
  • Initial-number-wolves: 草原上狼的初始数目
  • Sheep-gain-from-food: 羊吃草之后的能量增加值
  • Wolf-gain-from-food: 狼吃羊之后的能量增加值
  • Sheep-reproduce: 羊的出生概率
  • Wolf-reproduce: 狼的出生概率
  • Grass?: 是否包含草的生长计算,包括的话则是上面所说的第二种模
  • Grass-regrowth-time: 草重新生长的时间
  • Show-energy: 是否显示动物的能量值


现在按一下setup按钮,你会发现窗口中出现了绿草、肥羊和恶狼。羊和狼的数目是由Initial-number-sheep和Initial-number-wolves这两个滑动条所决定的,调整滑块,再按下setup按钮,就会按照新的设置初始化羊和狼的数目。



按下go按钮则开始运行模型。请注意go按钮的右下方有两个小箭头,这表示这个按钮是一个持续按钮,就是说按下go按钮之后,它将一直处在被按下的状态,同时羊和狼则会一直在草原上游荡并且捕食。再按一下go按钮则暂停运行。

现在让我们来看看羊有吃不完的草的情况。确定Grass?为OFF状态,依次按下setup和go按钮,并观察羊和狼的数量曲线的变化。结局是否有些出乎你的意料?我们让羊有无限的草吃,让狼有无限的羊吃,结果却令它们都灭绝了。通过观察数量曲线,你应该能很快发现奥妙所在吧。



让我们再看看草会被吃掉并且再生长的情况。确定Grass?为ON的状态,再次运行模型。这次你会发现草、羊和狼的数量相对稳定起来。这似乎告诉我们一个有趣的规律:越复杂的系统越稳定。



现在你可以着手调节各个参数,看看它们对生物们的数量的影响。并且完成下面几道思考题:

  • 当模型中只包含羊和狼的时候,能否找到一组参数配置,使得其数量保持平衡?
  • 设置 Grass?为ON,但是把狼的初始数目设置为0,这样就是一个羊和草的模型,这个模型稳定吗?
  • 当模型中包括草、羊和狼的时候,数量趋于稳定,但是我们可以观察到一定的变化周期,你能够通过调节参数扩大或者减小这种周期效应吗? 

Sudoku求解器

Sudoku中文叫做数独,是一种在日本很流行的数字谜题游戏。手工解决这种谜题需要用很多逻辑推算,我对这个一直不再行,只会几种最简单的逻辑,而且经常忘记和算错,所以一直想写一个求解程序来满足一下。

前两天用传播与分配法写了一个矩形谜题游戏的求解程序之后,感觉比较爽,所以今天就花了不少时间来做这个Sudoku solver 。这个程序同样用传播分配算法。传播算法就是用逻辑尽量减小每个方块的取值范围。在网上搜索到 Sudoku Solver by Logic这个网站,它试图尽量采用逻辑来解Sudoku,我从中学到了不少解Sudoku的高级逻辑,不过用程序来实现这些逻辑的效率反而不如让程序直接搜索。所以我写的这个求解器,只用最简单的逻辑进行传播计算,其余的就交给分配器搜索了。

由于我不太会用Javascript,所以这个程序是用Python写的,在服务器上运行,由于服务器是旧笔记本电脑(奔2),所以要等2,3秒才能求得答案。


用传播与分配算法解矩形谜题

谜题如下图所示,每个数字N表示一个矩形区域R, R所包含的小方格数正好为N。要求计算出所有数字对应的矩形区域,使得它们之间不互相覆盖,并且正好包含所有的方块。

        8
    
    8
      6
  
            
6
      6
     
    6
     9
   6
            
6
     8
     6
 
   6
     8
    
10
         6
  
   6
    6
  5
   6
     6
       
    6
      4
  

解决这种问题的最简单的办法就是回溯法,按照顺序把每个数字对应的矩形区域计算出来,如果当前计算出来的矩形区域和已经放置的矩形重叠的话,就尝试下一个候选矩形,直到全部区域被填满。题中一共有22个数字,每个数字对应的可能的矩形区域平均也有5,6种,这样下来,回溯法要花很长时间才能找到答案。

本文给出一种快速算法:传播与分配算法,此算法的基本知识,请参照文章:Constraint Programming in Oz

首先我们需要得到所有的数字所对应的所有可能矩形,我们称之为解空间,通过传播与分配算法,逐渐缩小这个解空间的取值范围,直到找到最终解。

解空间我们用英文status来表示,下面的程序计算出问题的初始解空间。

blocks列表是一个3个元素的的组元,(横坐标、纵坐标、此坐标上的数字)。

WIDTH = 12  
HEIGHT = 12
blocks = [
(0, 8, 10),(8, 4, 9),(7, 0, 8),(3, 1, 8),(5, 6, 8),
(7, 7, 8), (8, 9, 5),(11,9, 6),(4, 10,6),(3, 11,6),
(9, 11,4),(9,1, 6),(0, 3, 6),(6, 3, 6),(3, 4, 6),
(11,4, 6),(0, 6, 6),(10,6, 6),(2, 7, 6),(9, 8, 6),
(2, 9, 6),(6, 9, 6)
]

possibleBlock2函数计算出所有包含坐标为(x,y)的方块的,长宽为(w,h)的矩形区域,它返回的是有4个元素的组元,其意义分别为(矩形区域的左上点横坐标, 矩形区域的左上点纵坐标, 矩形区域的宽, 矩形区域的高)。

def possibleBlock2(x, y, w, h):
for i in range(-w+1, 1):
for j in range(-h+1, 1):
if x+i>=0 and y+j>=0 and x+i+w<=WIDTH and y+j+h<=HEIGHT:
yield (x + i, y + j, w, h)

possibleBlock函数计算出坐标为(x,y)的数字c所对应的所有可能的矩形区域。例如

>>> print list(possibleBlock(0, 8, 10))
[(0, 0, 1, 10), (0, 1, 1, 10), (0, 2, 1, 10), (0, 4, 2, 5),
(0, 5, 2, 5), (0, 6, 2, 5), (0, 7, 2, 5), (0, 7, 5, 2),
(0, 8, 5, 2), (0, 8, 10, 1)]
def possibleBlock(x, y, c):
for i in range(1, c+1):
if c%i == 0:
for p in possibleBlock2(x, y, i, c/i):
yield p

makeInitStatus函数通过收集所有的数字对应的所有可能的矩形区域,返回本问题的初始解空间status。status[i]为blocks[i]对应的所有可能的解。

def makeInitStatus(blocks):
status = []
for x, y, b in blocks:
status.append( list(possibleBlock(x,y,b)) )
return status

传播和分配的基本思想是通过传播操作,尽量缩小解空间的取值范围,在无法缩小的情况下,进行解空间的分配操作,也就是把解空间分解成两个子空间,然后分别对其进行传播和分配操作,直到找到解为止。

我们先来看看传播操作。对于本问题,传播操作就是剔除掉所有相互重叠的矩形。具体地说,就是当某个矩形R与某个数字N对应的所有的矩形都重叠的时候,从解空间中删除矩形R,下文称这种矩形为无效矩形。下面是这部分代码的说明。

rectangleOverlap函数判断矩形rect1和矩形rect2是否有重叠部分,有则返回True,无则返回False。

def rectangleOverlap(rect1, rect2):
x1, y1, w1, h1 = rect1
x2, y2, w2, h2 = rect2
if y1 + h1 - 1 < y2: return False
if y1 > y2 + h2 - 1: return False
if x1 + w1 - 1 < x2: return False
if x1 > x2 + w2 - 1: return False
return True

rectangleOverlapAll函数判断矩形rect和矩形列表rectlist中的所有矩形是否都重叠,都重叠时返回True。

def rectangleOverlapAll(rect, rectlist):
for rect1 in rectlist:
if not rectangleOverlap(rect, rect1):
return False
return True

removeOverlapRectangles函数从矩形列表rectlist1中删除掉与举行列表rectlist2中所有矩形重叠的矩形,如果删除了某个矩形的话,返回True,否则返回False。

def removeOverlapRectangles(rectlist1, rectlist2):
flag = False
for i, x in enumerate(rectlist1):
if rectangleOverlapAll(x, rectlist2):
rectlist1[i] = 0
flag = True
while 0 in rectlist1:
rectlist1.remove(0)
return flag

delImpossibleChoice函数完成对解空间status的传播操作,从解空间status中删除掉所有无效矩形,直到找不到更多的无效矩形。tocheck是一个待检查的集合,它的元素为一个2元组元(x, y),表示要删除status[x]中所有与status[y]重叠的矩形。当removeOverlapRectangles返回Ture时,表示status[x]发生了变化,因此需要重新检查(i, x),i为除了x之外的所有下标。当tocheck为空时,则表示无法继续进行传播操作了。此时的解空间的大小将明显的减少。例如本问题的初始解空间中共有个345矩形元素,经过一次传播操作之后剩下156个矩形元素。

def delImpossibleChoice(status):
tocheck = Set()
for i in range(len(status)):
for j in range(len(status)):
if i!=j: tocheck.add((i,j))
while len(tocheck) > 0:
x, y = tocheck.pop()
r = removeOverlapRectangles(status[x], status[y])
if r:
for i in range(len(status)):
if i!=x:
tocheck.add((i,x))

当传播操作无法继续时,进行分配操作,通过分配算法将当前的解空间分解为两个子解空间,然后再递归执行传播分配操作。下面是这部分程序的解释。

divideStatus函数对解空间status进行分配,得到两个子解空间status1和status2。分配的算法有很多种,这里采用的具体算法如下:
当某个数所对应的矩形列表长度为1时表示这个数所对应的矩形已经有解,当长度大于1时表示还多个候选解,找到候选解最少的下标index,然后把这组候选解[x0, x1, x2,...,xn]分为两部分:[x0], [x1, x2, ... , xn]。程序中的status1和status2是分解后的子解空间,显然status1和status2的集合就是status。

def divideStatus(status):
m = 100
index = -1
for i, x in enumerate(status):
if len(x) > 1 and m > len(x):
m = len(x)
index = i
status1 = deepcopy(status)
status2 = deepcopy(status)
status1[index] = [status[index][0]]
status2[index] = status[index][1:]
return status1, status2

isSolvedStatus函数判断status是否是最终解。当解空间status中所有的候选矩形数都为1时,则找到了最终解。

def isSolvedStatus(status):
for x in status:
if len(x)!= 1:
return False
return True

isImpossibleStatus函数判断status是否是无解解空间。 当解空间status中某组候选矩形数为0时,表示此空间无解。

def isImpossibleStatus(status):
for x in status:
if len(x) == 0:
return True
return False

solvePuzzle函数递归调用自己来实现传播分配操作,直到找到解。

def solvePuzzle(status):
delImpossibleChoice(status) #分配操作,尽可能地缩小解空间
if isSolvedStatus(status): #判断解空间是否为最终解
#printSolveStatus(status)
print status
return
if isImpossibleStatus(status): #判断解空间是否无解
return
#对解空间status进行分配操作,得到子解空间status1和status2
status1, status2 = divideStatus(status)
solvePuzzle(status1) #递归调用solvePuzzle搜索解空间status1中的解
solvePuzzle(status2) #递归调用solvePuzzle搜索解空间status2中的解

最后我们来运行一下程序:

>>> status = makeInitStatus(blocks)
>>> solvePuzzle(status)
[[(0, 7, 2, 5)], [(8, 2, 3, 3)], [(0, 0, 8, 1)], [(2, 1, 4, 2)],
[(5, 3, 1, 8)], [(6, 4, 2, 4)], [(8, 7, 1, 5)], [(11, 6, 1, 6)],
[(4, 5, 1, 6)], [(2, 11, 6, 1)], [(9, 10, 2, 2)], [(8, 0, 3, 2)],
[(0, 1, 2, 3)], [(6, 1, 2, 3)], [(2, 3, 3, 2)], [(11, 0, 1, 6)],
[(0, 4, 2, 3)], [(8, 5, 3, 2)], [(2, 5, 2, 3)], [(9, 7, 2,3)],
[(2, 8, 2, 3)], [(6, 8, 2, 3)]]

printSolveStatus函数将解打印成易读的形式:

def printSolveStatus(status):
chars = "abcdefghijklmnopqrstuvwxyz"
board = []
for i in range(HEIGHT):
board.append([" "]*WIDTH)
for cnt, rect in enumerate(status):
x, y, w, h = rect[0]
for i in range(x, x+w):
for j in range(y, y+h):
board[j][i] = chars[cnt]
for line in board:
print "".join(line)

最终的输出结果为:

cccccccclllp
mmddddnnlllp
mmddddnnbbbp
mmoooennbbbp
qqoooeffbbbp
qqssieffrrrp
qqssieffrrrh
aassieffgtth
aauuievvgtth
aauuievvgtth
aauuievvgkkh
aajjjjjjgkkh

在我的电脑上整个计算过程不过2,3秒。我用回溯法解了一下这个题,睡了一觉才等到答案,所以不知道到底花了多长时间。:-)

完整程序下载


rockbox语言文件转换器

rockbox的语言文件是二进制格式的,直接编辑其中的文字比较麻烦,所以我写了这个Python的小程序先把它转换为UTF-8的文本格式,用WindowXP的记事本编辑之后,再将其转换为二进制格式。运行这个程序需要安装Python,所以先到Python主页 下载一个Python的解释器,然后才能运行它。(点击此处 直接下载Python2.4.4 Windows版)

程序下载地址

解压上面的转换程序之后,启动Windows的命令行,CD到其所在目录,把需要编辑的chinese-simp.lng文件也复制到这个目录下。然后输入convert.py运行此程序。按照下面的例子输入:

D:\>convert.py
1: convert lng file to text file
2: convert text fle to lng file
3: update dictionary by text file
>>1
input lng file name>>chinese-simp.lng
output text file name>>chinese-simp.txt
use dictionary?(y/n)>>n  

这样就得到了一个chinese-simp.txt的文本文件,在此文件中每行是一条数据,其格式如下:

0|:-:|是|:-:|
1|:-:|否|:-:|

在每行的 |:-:|之后输入新的字符串即可,例如:

 1|:-:|否|:-:|不是

然后按照下面的例子把编辑后的chinese-simp.txt再转换为lng文件

D:\>convert.py
1: convert lng file to text file
2: convert text fle to lng file
3: update dictionary by text file
>>2
input text file name>>chinese-simp.txt
output lng file name>>chinese-simpnew.lng

这个chinese-simpnew.lng就是所要的结果了。


读写Rockbox的字体

前些天发现本博客提供的simsun字体中的标点符号上下颠倒了,在阅读文本的时候看上去很不舒服,所以一直想修正这个错误,今天终于找到一点空闲时间了。

由于bdf文件为文本文件,修改起来比较麻烦,所以我决定直接用程序读取rockbox的fnt字库格式。Rockbox网站上可以下载到完整的代码,因此分析fnt字库的格式也没有什么困难的。

有关fnt的读取都是在font.c中完成的,fnt格式在font.h中有介绍:

/*
* .fnt loadable font file format definition
*
* format len description
* ------------------------- ---- ------------------------------
* UCHAR version[4] 4 magic number and version bytes
* USHORT maxwidth 2 font max width in pixels
* USHORT height 2 font height in pixels
* USHORT ascent 2 font ascent (baseline) in pixels
* USHORT pad 2 unused, pad to 32-bit boundary
* ULONG firstchar 4 first character code in font
* ULONG defaultchar 4 default character code in font
* ULONG size 4 # characters in font
* ULONG nbits 4 # bytes imagebits data in file
* ULONG noffset 4 # longs offset data in file
* ULONG nwidth 4 # bytes width data in file
* MWIMAGEBITS bits nbits image bits variable data
* [MWIMAGEBITS padded to 16-bit boundary]
* USHORT offset noffset*2 offset variable data
* UCHAR width nwidth*1 width variable data
*/

根据这个格式,我写了一个小python的测试程序验证,其中的变量f是某个打开的字库文件:

version = f.read(4)
maxwidth, height, ascent, pad = unpack("HHHH",f.read(8))
firstchar, defaultchar, size, nbits, noffset, nwidth = unpack("LLLLLL", f.read(6*4))
if nbits >= 0xffdb:
nbits1 = (nbits/4+1)*4
else:
nbits1 = (nbits/2+1)*2

bits = f.read(nbits1)
offset = array("L", f.read(noffset*4))
width = array("B", f.read(nwidth))

其中bits储存的就是所有的点阵字形的信息,要找到某个文字的点阵字形的话,需要用到offset和width。fnt字库是按照文字的UCS-2编码的,所以首先要知道某个文字的UCS-2编码才能找到与其对应的字形信息。Rockbox中有完整的编码转换程序,可以参考一下。我用的办法是用notepad++ ,把文件编码改为UCS-2编码,然后输入我要查找的文字,然后改为16进制查看方式,就可以知道文字对应的UCS-2编码了。下面的Python程序找到文字对应的字形信息,参数char是一个整数,表示字符的UCS-2编码,例如逗号“,”的编码是0xff0c。

那么得到的字形信息是如何表示文字的二位点阵的呢?在Gray_draw.c中函数gray_mono_bitmap_part的前面有这么一段说明:

/* About Rockbox' internal monochrome bitmap format:
*
* A bitmap contains one bit for every pixel that defines if that pixel is
* foreground (1) or background (0). Bits within a byte are arranged
* vertically, LSB at top.
* The bytes are stored in row-major order, with byte 0 being top left,
* byte 1 2nd from left etc. The first row of bytes defines pixel rows
* 0..7, the second row defines pixel row 8..15 etc. */

根据这段说明我写了一个把rockbox 内部monochrome bitmap格式转换为二维点阵的程序。

01def readCharList(char):
02    bitmapOffset = offset[char]
03    srcBytes = width[char] * ((height + 7) / 8)
04    return bits[bitmapOffset:bitmapOffset+srcBytes]
05 
06def getBitmap(char):
07    charlist = readCharList(char)
08    w = width[char]
09    h = height
10    b = [[' ' for x in range(w)] for y in range(h)]
11    for i, c in enumerate(charlist):
12        for j in range(8):
13            v = (ord(c) >> j) & 0x0001
14            if v:
15                b[j+(i/w)*8][i%w] = '#'
16    return b
17 
18def printBitmap(bitmap):
19    for r in bitmap:
20        print "".join(r)
 

下面是一个显示的例子,printBitmap( getBitmap( 0x597d )),0x597d为“好”的UCS-2编码。
  #
# ######
# #
##### #
# # #
# # #
# #########
# # #
# # #
# #
# # #
# # #
# ##
有了上面这段程序,就很容易把标点符号纠正过来了。点击下载rockbox simsun 字库

分析-外星人计算Pi的程序

* 本文作于2002年,最早发表于网易论坛。现在用Google搜索一下,真还有不少转载的。

源程序

本文分析下面这个很流行的计算PI的小程序。下面这个程序初看起来似乎摸不到头脑,不过不用担心,当你读完本文的时候就能够基本读懂它了。

程序一:很牛的计算Pi的程序
int a=10000,b,c=2800,d,e,f[2801],g;
main() {
for(;b-c;)
f[b++]=a/5;
for(;d=0,g=c*2;c -=14,printf("%.4d",e+d/a),e=d%a)
for(b=c; d+=f[b]*a,f[b]=d%--g,d/=g--,--b; d*=b);
}

数学公式

数学家们研究了数不清的方法来计算PI,这个程序所用的公式如下:

          1          2          3                    k
pi = 2 + --- * (2 + --- * (2 + --- * (2 + ... (2 + ---- * (2 + ... ))...)))
3 5 7 2k+1

至于这个公式为什么能够计算出PI,已经超出了本文的能力范围。下面要做的事情就是要分析清楚程序是如何实现这个公式的。我们先来验证一下这个公式:

程序二:Pi公式验证程序
#include "stdio.h"
void main()
{
float pi=2;
int i;
for(i=100;i>=1;i--)
pi=pi*(float)i/(2*i+1)+2;
printf("%f\n",pi);
getchar();
}

上面这个程序的结果是3.141593。

程序展开

在正式分析程序之前,我们需要对程序一进行一下展开。我们可以看出程序一都是使用for循环来完成计算的,这样做虽然可以使得程序短小,但是却很难读懂。根据for循环的运行顺序,我们可以把它展开为如下while循环的程序:

程序三:for转换为while之后的程序
int a=10000,b,c=2800,d,e,f[2801],g;
main() {
int i;
for(i=0;i<c;i++)
f[i]=a/5;
while(c!=0)
{
d=0;
g=c*2;
b=c;
while(1)
{
d=d+f[b]*a;
g--;
f[b]=d%g;
d=d/g;
g--;
b--;
if(b==0) break;
d=d*b;
}
c=c-14;
printf("%.4d",e+d/a);
e=d%a;
}
}

注:for([1];[2];[3]) {[4];}的运行顺序是[1],[2],[4],[3]。如果有逗号操作符,例如:d=0,g=c*2,则先运行d=0,然后运行g=c*2,并且最终的结果是最后一个表达式的值,也就是这里的c*2。

下面我们就针对展开后的程序来分析。

程序分析

要想计算出无限精度的PI,我们需要上述的迭代公式运行无数次,并且其中每个分数也是完全精确的,这在计算机中自然是无法实现的。那么基本实现思想就是迭代足够多次,并且每个分数也足够精确,这样就能够计算出PI的前n位来。上面这个程序计算800位,迭代公式一共迭代2800次。

int a=10000,b,c=2800,d,e,f[2801],g;

这句话中的2800就是迭代次数。

由于float或者double的精度远远不够,因此程序中使用整数类型(实际是长整型),分段运算(每次计算4位)。我们可以看到输出语句 printf("%.4d",e+d/a); 其中%.4就是把计算出来的4位输出,我们看到c每次减少14( c=c-14;),而c的初始大小为2800,因此一共就分了200段运算,并且每次输出4位,所以一共输出了800位。由于使用整型数运算,因此有必要乘上一个系数,在这个程序中系数为1000,也就是说,公式如下:

               1          2          3                    k
1000*pi = 2k+ --- * (2k+ --- * (2k+ --- * (2k+ ... (2k+ ---- * (2k+ ... ))...)))
3 5 7 2k+1

这里的2k表示2000,也就是f[2801]数组初始化以后的数据,a=10000,a/5=2000,所以下面的程序把f中的每个元素都赋值为2000:

for(i=0;i<c;i++) f[i]=a/5; 

 

你可能会觉得奇怪,为什么这里要把一个常数储存到数组中去,请继续往下看。我们先来跟踪一下程序的运行:

while(c!=0)  假设这是第一次运行,c=2800,为迭代次数
{
d=0;
g=c*2; 这里的g是用来做k/(2k+1)中的分子
b=c; 这里的b是用来做k/(2k+1)中的分子
while(1)
{
d=d+f[b]*a; f中的所有的值都为2000,这里在计算时又把系数扩大了
a=10000倍。

这样做的目的稍候介绍,你可以看到输出的时候是d/a,所以这不影计算
g--;
f[b]=d%g; 先不管这一行
d=d/g; 第一次运行的g为2*2799+1,你可以看到g做了分母
g--;
b--;
if(b==0) break;
d=d*b; 这里的b为2799,可以看到d做了分子。
}
c=c-14;
printf("%.4d",e+d/a);
e=d%a;
}

只需要粗略的看看上面的程序,我们就大概知道它的确是使用的那个迭代公式来计算Pi的了,不过不知道到现在为止你是否明白了f数组的用处。如果没有明白,请继续阅读。

d=d/g,这一行的目的是除以2k+1,我们知道之所以程序无法精确计算的原因就是这个除法。即使用浮点数,答案也是不够精确的,因此直接用来计算800位的Pi是不可能的。那么不精确的成分在哪里?很明显:就是那个余数d%g。程序用f数组把这个误差储存起来,再下次计算的时候使用。现在你也应该知道为什么d=d+f[b]*a;中间需要乘上a了吧。

把分子扩大之后,才好把误差精确的算出来。d如果不乘10000这个系数,则其值为2000,那么运行d=d/g;则是2000/(2*2799+1),这种整数的除法答案为0,根本无法迭代下去了。

现在我们知道程序就是把余数储存起来,作为下次迭代的时候的参数,那么为什么这么做就可以使得下次迭代出来的结果为接下来的4位数呢?这实际上和我们在纸上作除法很类似:

        0142
/——------
7 / 1
10
7
---------------
30
28
---------------
20
14
---------------
60
.....

我们可以发现,在做除法的时候,我们通常把余数扩大之后再来计算,f中既然储存的是
余数,而f[b]*a;则正好把这个余数扩大了a倍,然后如此循环下去,可以计算到任意精
度。
这里要说明的是,事实上每次计算出来的d并不一定只有4位数,例如第一次计算的时候
,d的值为31415926,输出4位时候,把低四位的值储存在e中间,e=d%a,也就是5926。

最后,这个c=c-14不太好理解。事实上没有这条语句,程序计算出来的仍然正确。只是
因为如果迭代2800次,无论分数如何精确,最后Pi的精度只能够达到800。
你可以把程序改为如下形式尝试一下:

for(i=0;i<800;i++)
{
d=0;
g=c*2;
b=c;
while(1)
{
d=d+f[b]*a;
g--;
f[b]=d%g;
d=d/g;
g--;
b--;
if(b==0) break;
d=d*b;
}
// c=c-14; 不要这句话。
printf("%.4d",e+d/a);
e=d%a;
}

最后的答案仍然正确。

不过我们可以看到内循环的次数是c次,也就是说每次迭代计算c次。而每次计算后续位数的时候,迭代次数减少14,而不影响精度。为什么会这样,我没有研究。另外最后的e+d/a,和e=d/a的作用就由读者自己考虑吧。

 


Constraint Programming in Oz

本文通过几个简单的例子,介绍约束逻辑程序设计的基本知识,以及如何用Mozart-Oz来实现。

(本文中将把Constraint Programming翻译为约束程序设计,译者本人并不知道许多专有名词的译法,因此可能有谬误,请见谅)

介绍

1. 什么是约束程序设计?

我们还是从一个例子着手吧。
例1.
有如下的数学算式:SEND+MORE=MONEY,每个字母代表一位从0到9的数字,首位不能为0,每个字母所代表的数字不同,请问各个字母所代表的数字是多少?
这个问题有唯一解答:9567+1085=10652
这就是一个条件约束问题。解决此类问题最简单的方法可以描述为"产生并测试"。我们很容易看出这个问题
有8个字母,而每个字母又可以是从0到9,即有10种可选方案。所以总共的可选答案是8^10种。产生程序将产生这8^10种可能解答, 然后使用测试程序来测试所产生的解答的正确性。下面列出这个问题的约束条件:
1. 8个字母所代表的数字两两不同。
2. S M不能为零
3. S*1000+E*100+N*10+D + M*1000+O*100+R*10+E = M*10000 + O*1000+N*100+E*10+Y
只要产生的某种组合符合上述三个条件就可以说是找到答案了。
但是显然这种方法在多数情况下是行不通的。这是由于组合问题中通常都存在组合爆炸,例如上题中如果不是8个字母, 而是9个的话,那么可能的组合数目就是9^10,字母仅仅只多了一个,而要产生的排列却增加了3倍多。
本文将详细介绍如何编制解决此类问题的程序,不过我们不使用传统的程序设计语言,而是使用一种已经包含了完整解决方案的语言:OZ。

2.基本原理

首先来看看需要使用什么方法才能够减少测试的次数吧。我们从另外一个简单例子入手。典型的条件约束问题具有如下的形式:
它 包含n个变量:V1,......,Vn, 每个变量的取值范围分别是:D1,......Dn,其中D1,......Dn都是包含有限元素的集合,通常集合中的元素都可以使用整数来表示。目标是 找到一组满足约束条件C(V1,......,Vn)的解。下面是一个简单的例子:
例2.有15个变量V1,......,V15,他们的取值范围都是{1,......,15},约束条件为:V1
V14
很明显这个问题只有一个答案就是Vk=k。我们来看如何让电脑来解决这个问题吧。
首先使用通常的产生并测试的方法:
我们可以通过下面的方法来产生一系列的可能组合:
选择一个尚未赋值的变量,从该变量的值域中选择一个值赋给此变量,重复此步骤直到所有的变量都赋上了值。这个操作将生成一个树状结构:每个非叶节点都代表一个尚未完全赋值的状态,叶节点则代表了一种可能的组合。
满足约束条件的叶节点就是问题的解,不满足约束条件的叶节点就是失败点。在本文中我们将经常使用图形来表达这种树状结构:
其中蓝色圆形代表选择点,即非叶节点。红色方形代表失败点,绿色菱形代表解节点。为了方便起见,一棵所有叶节点都是失败点的子树使用红色三角形代表,而包括了解节点的子树则使用绿色三角形代表。

现在来看看前面的问题:一共有15个变量,每个变量的可能取值有15个,那么总共的可能组合就是15^15=437,893,890,380,859,275种, 就目前的计算机来说这是个天文数字了。

3. 交错产生和测试

前 面所介绍的"产生并测试"的方法,是在完全产生了一种组合之后,再来测试这个组合。我们可以把产生和测试交错进行,这样可以极大地提高效率。在每一次选择 之后都使用约束条件来判断这步选择是否合理,而不是完全产生了某种组合之后再来判断。例如在上面的例子中,当给V1和V2赋值之后,就可以立即使用V1
传播和分配(Propagate and Distribute)
到目前为止,我想你应该已经认识到使用简单的"产生并测试"的方法是行不通的,而"交错产生和测试"虽然能够提高效率,但是也有它的局限性:若使用这种方法来解例2,那么搜索树将包含917477个节点。如果变量数目增加到100又会如何呢?:)
那么使用什么方法呢?我们不希望测试过程只是被动的检测某种组合是否满足条件,而希望测试过程是积极的。这就是说通过测试,我们希望缩小变量的取值范围。在前面的例子中,V1和V2的取值范围都是1到15,但是由于V1
然而在通常情况下搜索是不可避免的,我们的目的是在每次搜索之前先通过约束条件尽量缩小变量的取值范围。这种方法在约束程序设计中被叫做"传播和分配(propagate and distribute)"。
传播通过简单的确定性的推导过程尽量减少可能的解的个数。当传播过程无法继续的时候,就需要进行分配。所谓分配就是进行非确定性的选择(下面将详细介绍,如果有何不明白请继续阅读),这样搜索树的尺寸能够尽量缩小,我们可以把传播的过程看作对搜索树的裁剪。

Mozart OZ

OZ 是一门程序设计语言,而Mozart则是它的实现系统之一。OZ是一门新兴的程序设计语言,因此它也包含了许多特别的技术。自然,约束程序设计是它的主要 功能之一,除此之外它有像Prolog一样的逻辑程序设计模式,有Lisp, SML式的函数式程序设计模式,另外还有并行计算,多线程,面向对象等模式。本文主要介绍它在约束程序设计方面的运用,至于其它方面将在必要的时候进行讲 解。

1. 安装

Mozart支持许多操作系统,你可以到www.mozart-oz.org下载相应的版本。除此之外你的电脑上必须安装Emacs(Xemacs),Mozart将使用这个编辑器作为开发界面。安装中如果出现什么问题,请参照mozart的帮助文档。
下面是运行界面,最上面的窗口就是我们写程序的地方。中间一部分(标有Oz compiler)用来显示编译信息,最下面一行是Emacs的提示栏,用来显示和输入一些编辑信息。

我 们先来看如何做HelloWorld程序,在最上面的窗口输入下面文字之后,按Ctrl+.然后按Ctrl+b就会把当前窗口中的程序送交编译器编译运 行,当然也可以选择OZ菜单中的feed buffer选项,另外还有Feed Region,Feed Line,Feed Paragraph等选项,Region表示当前的本选择的文字,Line表示光标所在的行,Paragraph表示一段程序(没有空行分隔)。我们看到 编译信息窗口出现了编译信息,accepted表示编译通过,可是运行的结果在哪里呢?运行的结果在一个尚未打开的buffer中,选择OZ-> Show/Hide->Emulator就可以看到结果了。

把上面的Show改为Browse,即{Browse 'Hello World!'}并feed buffer之后,将看到
窗口出现。此后我们都将使用OZ browser来显示结果。

2. 传播

OZ包含了整个约束程序设计所需要的功能,我们想来看看传播(propagate),即OZ是如何最大限度的缩小变量的取值范围的。
先输入如下程序(把上面的Oz Browser,先关闭,把程序buffer中的文字也请清空)

declare X Y Z
X::1#10
Y::1#10
Z::1#10
{Browse [X Y Z]}

然后用Feed Paragraph把上面程序段送交编译器。
declare 关键字用来声明变量,这里声明了三个变量标识符:X Y Z,OZ中所有的变量标识符都必须大写。X::1#10意思是说X的取值范围从1到10,包括1和10。::传播器是用来定义约束变量的取值范围的,而 1#10则是某个可能的取值空间。最后我们用Browse把三个变量的值显示出来,[X Y Z]把这三个变量做成一个列表,关于列表,如果你熟悉Prolog,Lisp的话应该不会陌生。出现的结果如下:

我们看到变量标识符后面就是此变量的取值范围。
此后我们每输入一句语句,就用Feed Line把它送交编译器,这样可以动态地看到程序运行的结果,Browse中的显示内容不再用图表示,我在前面加上>>前缀表示是Browser中的内容。

X<:Y
>> [X{1#9} Y{2#10} Z{1#10}]

输 入X<:Y之后我们惊奇的发现,Browser中的内容变化了,这也是Brwoser的强大功能之一,它并不只是简单的把结果显示出来,当变量发生 变化的时候,这些变化也会立刻在Browser中体现出来。X<:Y是一个一个传播器,当我们把这个约束条件送交编译器后,传播起作用了,同时缩小 了X Y的取值范围。

Y<:Z
>> [X{1#8} Y{2#9} Z{3#10}]

输入Y<:Z后我们发现X的取值范围也发生了变化,这也就是说X<:Y这个传播器也起了作用。

X+Y+Z=:6
>>[1 2 3]

我们发现这个约束条件加上之后,X Y Z的值就都确定下来了。

我们看到了传播器可以自己进行推导,尽量的减小变量的取值范围。不过这种推导也是非常复杂的,因此Mozart在这方面做了某些限制。例如Y=:2*X并不能够推导出Y必须是偶数,
除了上面的这些简单的传播器之外,Mozart还提供了许多很高级的约束传播器,等我们需要的时候再进行讲解。

3. 分配

当传播器无法在缩小变量的取值范围的时候,就需要进行分配了。所谓分配就是把变量组中的某个变量的取值范围分成两个部分,然后分别再对这两个部分进行传播和分配。
下面我们先来看一个完整的传播分配的例子:
我们要解决的问题如下:

X Y Z属于{0,1,...,7}
X+Y=3*Z
X-Y=Z

首先我们可以根据上面的约束条件写出如下的传播器:

X::1#7
Y::1#7
Z::1#7

上面三个传播器可以缩写为: [X Y Z] ::: 1#7,我们用列表把三个变量扩起来,然后用:::定义这个列表的取值范围,即列表中的每个元素的取值范围。剩下的两个等式的传播器如下:

X + Y =: 3*Z
X - Y =: Z

如果你把上述的传播器传送给编译器的话,

declare X Y Z
[X Y Z]:::1#7
X+Y=:3*Z
X-Y=:Z
{Browse [X Y Z]}

将得到如下结果:

>> [X{2#7} Y{1#6} Z{1#4}]

很显然,并没有找到答案。
因 此我们要使用分配器,把上述变量中的某个变量的取值范围分成两个部分,然后再分别对这两个部分进行传播。分配的方式有很多, 最简单的方式就是选择变量组中的第一个变量,假设这个变量的取值范围是a#b,那么分开的两组范围就是a和(a+1)#b。这种分配方式叫做naïve。 分配语句如下:

{FD.distribute naive [X Y Z]}

这句话表示我们要对变量组X Y Z的值域进行分配,分配方式是naïve。
我们先来看看实现程序:

declare
proc {Equations Sol}
X Y Z
in
Sol = solution(x:X y:Y z:Z)
% Propagate
[X Y Z] ::: 1#7
X + Y =: 3*Z
X - Y =: Z
% Distribute
{FD.distribute naive [X Y Z]}
end
{Explorer.all Equations}

这里我们声明了一个过程,过程的声明方式如下:

proc {Equations Sol}
end

Equations为过程名,Sol为过程的参数,这里这个参数是作为返回值用的。
X Y Z是这个过程的内部变量标识符,我们用一个结solution(x:X y:Y z:Z)把这三个变量包含起来,并赋予Sol,这样过程的返回值就是一组解了。
下面的分别定义传播器和分配器,%表示注释:

% Propagate
[X Y Z] ::: 1#7
X + Y =: 3*Z
X - Y =: Z
% Distribute
{FD.distribute naive [X Y Z]}

定义了上述函数之后,我们还必须把这个函数传递给专门的搜索引擎:

{Explorer.all Equations}
{Explorer.one Equations}

Explorer.all找Equations的所有解,而Explorer.one只搜索一个解。把上面的程序传送给编译器后,我们看到如下的结果:

我们看到Explorer显示的是一棵搜索树,按照上面树的节点的顺序一起双击,则节点的值在下面的Inspector

中显示出来。图中每个蓝色的圆形节点都表示一次分配的过程。我们看到节点1的值就是前面只进行传播之后的结果。由于得到这个结果之后无法再进行传播,因此分配器对变量组[X Y Z]中的第一个变量X进行分配,得到如下两组取值范围:

[X{2} Y{1#6} Z{1#4}]
[X{3#7} Y{1#6} Z{1#4}]

对[X{2} Y{1#6} Z{1#4}]进行传播最终得到一组解:2#solution(x:2 y:1 z:1),即图中的第二个节点,解都是由绿色的菱形表示。而节点3就是上面的第二组取值范围[X{3#7} Y{1#6} Z{1#4}]进行传播之后的结果[X{3#7} Y{1#5} Z{2#4}],我们看到Y和Z的值域发生了变化。因为传播器无法对[X{3#7} Y{1#5} Z{2#4}]这个取值范围进行传播以获得更小的取值范围,因此分配器再次进行分配:

[X{3} Y{1#5} Z{2#4}]
[X{4#7} Y{1#5} Z{2#4}]

传播器对第一组结果进行传播发现自相矛盾,因此我们得到一个红色的正方形节点,表示失败。对[X{4#7} Y{1#5} Z{2#4}]进行传播的结果为节点4。
如此下去当整个搜索树完成的时候,我们得到三个绿色的菱形节点,也就是三组解。
现 在我们知道了传播分配的基本原理,在OZ中提供各种各样的分配器。例如有些分配器是找到变量组中可能取值个数最少变量进行分配,这样有助于尽早确定这个变 量的值,而有的分配器是把变量的取值范围等分为两部分,这样做有助于传播器最大限度的对搜索树进行裁剪。不同的问题采用不同的分配器,得到结果最终是一样 的,不过效率就有很大的区别了。

传播与分配的原理

待续...


0 1