Graph Algorithm: Python Implementation

Because Numenta's algorithm require a graph algorithm(directed graph with weighted arcs), I find a lot of good links.

Dijkstra's algorithm for shortest paths
A comprehensive Python/Graph wiki
Simple directed graph essay(only DFS function)
Breath-first Search can be used as my Numenta Graph Algorithm implementation.

Read more ... ...

Python for Numerical Computing

Tentative NumPy Tutorial
NumPy for Matlab Users
Numpy Example List

Tuple T=(1,2,3,4), List L=[1,2,3,4]. List can be modified while Tuple cannot after it is created. So use Tuple as the key of dictionary.


Nothing else yet.

Read more ... ...

The Bigger Picture: A World Model

在这篇文章里，我们会站得远一点看，试着理解图片世界的一般特性，从而把简单的图片世界和复杂的真实世界联系起来。

图24:将图片”Cat”提交给HTM网络。工作在三个不同层片的不同节点的接受区域被表示出来。

考虑图.24的”Cat”。假设我们将这个图片提交给我们的网络来完成识别（译者注：指Numenta的HTM网络，HTM: Hierarchical Temporal Memory）。在这个图片里，我们可以在三个不同层次(level)的三个节点(node)上识别这个接受区域（receptive field）。顶层（top-level，也即第三层）的接受区域是整个一张图片。中间层(mid-level)的接受区域是一个中等大小的8x8像素块，底层（bottom-level）的接受区域则是比较小的4x4像素块。当这个图片被提交给HTM网络时，顶层节点会激活”Cat”状态，中间层节点会激活猫头（“Cat-head”）状态，底层节点则会激活直角（“Corner”）状态。简而言之，我们可以认为顶层，中间层和底层节点相信它们在这张图里分别看到了“Cat”，“Cat-head”和“Corner”。

图25：这幅图显示了视觉上的时间层次概念。箭头表示不同层次上的节点看到的不同状态的持续时间。与低层次相比，高层次的状态改变更慢一些。

现在让我们想象当这个“Cat”向上移动的时候会发生什么。顶层节点会依旧认为它看到的是“Cat”。在向上移动图片之后，底层节点会依次看到直角（“Corner”），转过90度的T字（rotated “T-junction”），垂直线（"vertical-line”）和空白（"blank”）。但是顶层节点在移动期间会始终相信自己看到的是“Cat”。中间层节点在此期间看到的事物也会发生一些变化。在一段时间内（这段时间比底层节点看到直角（“Corner”）的时间长，比顶层节点看到“Cat”的时间短），中间层节点会看到“Cat-head”。一段时间以后，也许是底层节点认为自己看到垂直线（"vertical-line”）的时候，中间层节点会认为自己看到了猫身（“cat-torso”）。

我们在这里观察到的现象是基于自然移动模式（natural movement pattern），也就是说，猫（Cat）不是随机地跳来跳去，那么较低层次上的节点切换状态的速度要比在较高低层次上的节点快一些。这个想法可以用图25表示。我们还知道与较低层次相比，较高层次包括更多的空间信息。这两者（译者注：即时间和空间）共同作用的结果是，与较低层次的节点相比，较高低层次的节点能够（通过较低层次的节点）表示更大的空间和时间跨度。在一个较高层次上表示的状态（state）被一个较大的空间跨度影响，与降低层次的状态相比将持续更长的时间。用图来表示空间组织和时间组织是相当困难的。我们尝试着画了图26。

图26：这幅图显示了空间和时间的抽象层次概念。越高的层次抽象越广的空间和越长的时间跨度。也就是说，高层次的状态变化比低层次慢一些。

我们可以上下倒转这个层次结构，把它想象成世界产生图像的过程（而不是把它想象成一个识别的层次结构）。我们的可视化的世界可以被想象成一个“生产模型”（“generative model”）采用一个层次结构来组成图像。这个模型工作在空间轴和时间轴上。当某个原因（cause）在较高层次上被触发，会在同时和在一些不同的地方，引起并展开许多的次级原因（sub-cause）。这些次级原因也会在更小的时间和空间跨度上，依次触发更细微的次级原因。所以当我们潜入到这个世界层次的下层的时候，会看到细微的空间和时间跨度。而在较低层次上产生的原因又反过来限制了在较高层次上变化的原因。在某一个顶层触发产生的原因的作用下，这种生成图像的机制会产生不同的图像。（译者注：意思是顶层原因会忽略掉一些下层原因的细节，所以会产生不同的图像。但是这些不同的图像归根结底都是由顶层原因所触发的。）

虽然我们用固定图像模式识别的例子引入了这个模型（译者注：指Numenta的HTM模型），这种多重跨度（multi-scale）组织结构在自然界是很普遍的。气候是一个比较好的例子。冬季就是一个高层次的状态影响整个国家。在冬季，纽约一般是大雪纷飞，同时在旧金山只会少几天阳光灿烂的日子。冬天，这个高层次的原因，在不同的地区产生了不同的次级状态。而且在冬季的一段日子里，下雪天和晴天此消彼长，相互交替出现。

让我们假设这个世界由一个时空交织的层次结构组成的，那么如果我们知道这个层次的结构，学习这样一个模型易如反掌。由于较高层次上的原因相比较较低层次上的原因持续较长的特性，这种层次结构使得我们可以利用时间来指导每一个节点来组成独特的群。因此，通过交替使用时间稳定的监督式学习（temporal stability based supervision）和一致性检测（coincidence detection），我们可以重建这个世界的模型。这就是HTM模型的一个基础思想。

译后记

这篇文章是”Numenta HTM learning Algorithms”, section 9 – A World Model的完整翻译。基本上概括了Hawkins在他“On Intelligence”一书里提到的人理解意义的关键过程。请注意，从脑科学的角度来讲，人的认知过程还是神秘的，所以Hawkins也承认他的书带有主观因素。然而，读一篇有启发性的文章总是令人愉悦的，我们不妨看看Hawkins和他的公司Numenta是如何将他的思想化为可以执行的算法。

正如这篇文章里所说的，理解意义的关键是“原因（cause）”和由原因引发的“次级原因（sub-cause）”是否是可以互相转化的。换句话说，HTM模型想要完成的是“一叶知秋”或者“窥一斑而知全豹”。当然“次级原因”的完整性也是不可或缺的要素。假如图片里只有一片秋叶，要问Numenta网络中分类器（Catalog）这片秋叶是从哪儿获得的，估计人的智能也无法完成这样的任务。

另外，时间也是一种理解和决定意义的关键。所有的有意义的事件在时间轴上连续变化，且绝大部分是可以预测的（量子层面上的变化可能需要排除，按照Einstein的话说，上帝也在扔色子决定量子的变化。但是我们仍能从概率的角度理解量子的变化。）。

又译后记

翻译文章是再创作。想当初，我在看傅雷译的“巴尔扎克全集”时，每每被傅氏的法式中文所震撼，并产生诸多感慨，希望有朝一日自己也能翻译些国外的文字，将国外的风貌和人情传递到国内来。然而自己终究在语言上少那么一根筋，之后去交大读了工程。

工作以后，有了闲暇时间，开始考虑自己下一步的方向。翻译这些技术性的小文章一方面是传递信息给国内的读者，另一方面是挖掘自己喜爱的题目（照我老妈的话说，不要忘了你的梦想）。现在我快找到了，是啊，有什么比研究人类自己更有趣的呢，你说呢:-)

Read more ... ...

Solving Every Sudoku Puzzles: Part 2

Contiune Solving Every Sudoku Puzzles: Part 2

数独(Sudoku)迷题的通用解法: Part 2

搜索

另一条路是通过搜索来得到答案：全面尝试所有的可能性知道我们恰巧得到一个可行的解答。这种方法的代码时候很少的几行，但是我们会面临另一个问题：有可能永远也算不完。让我们回到上面提到的hard puzzle，A2有4种可能（1679），A3有5种可能（12679）；那么现在一共有20种可能，如果我们连乘下去，对于这个迷题，我们会得到462838344192000000000000000000000000000 (或大约 4 ×10^38)种可能。你肯定吗？一定确定以及肯定，这里有61个待解开的方块，并且每一个这种方块都有4或5种可能。而且，事实上4^61<4x10^38<5^61。我们怎么来对付它呢？看来有两种选择。>首先，我们可以尝试蛮力法(brute force)。假设我们有一个非常智能的搜索算法可以用一条指令估计一个位置，而且我们有下一代的计算技术，就假设是一中1024核的10GHz处理器，我们买了一百万颗这样的处理器，当我们在购物的同时，还买了一台时间机器，帮助我们回到开天辟地的时候，让这个代码跑起来。就这道题目直到现在，我们大概可以计算过大概1%的可能值。

第二种选择是用某种方法使得每条机器指令处理估计一种以上的可能。这看起来不太可能，而幸运的是这种方法就是约束传递所作的。我们不用试过所有4 x 10^38种可能，因为我们只要试过一种，马上就可以消去很大一部分的可能值。例如，方块H7有两种可能，6和9。我们可以试试9，很快会发现有一个冲突，这意味着我们不仅消去了一种可能，而是4 x 10^38种可能的一半。

事实上，在解这个问题的时候我们只需要试25种可能值和9个方块（一共61个待解方块）；而约束传递解决了剩下来的问题。对于剩下来的95个hard puzzles，我们平均需要试64种可能值和不超过16个方块。

那什么是搜索算法呢？简单：先确定我们是不是已经得到结论了或者存在冲突，如果都不是，再选择一个待解方块并尝试它所有的可能值。一次一个，尝试给每一个方块赋所有的可能值，并且从已知的盘面开始搜索。换言之，我们搜索一个值d，使得我们可以成功的从将d赋值给方块中解出需要的解。这就是recursive search，我们将它称为depth-first搜索，因为我们（递归地）考虑values[s]=d是中所有，之后再考虑方块s的其它可能值。

为了防止bookkeeping nightmares，我们为每一次search递归调用复制一个新的拷贝。（译者注：bookkeeping nightmare不太理解，可能和实参形参有关）这样一来search tree的每一个branch都是独立的，不会相互干扰。（作者注：这就是为什么我选择将一个方块的可能值的集合设计成一个字符串’：我可以用简单有效的values.copy()来复制拷贝。而假如我用Python的set或list来实现这个可能值的集合，我不得不用copy.deepcopy(values)，而这个方法效率比较低。）另一种方法是保持values每一次改变的值，当我们碰到一次False的时候就将改变前的值恢复出来。这被称为backtracking search。这种方法在每一步都是单个改变（single change）的时候才有意义，而在我们的迷题的算法中有许多改变都来自于约束传递，这样的技巧会将我们引入复杂的泥潭。

所以，余下来的问题就是如何在搜索的每一步选择一个方块s来赋值。让我们回到上面的hard puzzle，假设我们选择B3，它有7种可能值（1256789），所以我们猜错的可能性是6/7。而，如果我们选择H7，它只有2种可能值（69），我们猜错的可能性是1/2。很明显，选择H7会更可能将我们引导到正确的答案，所以我们总是选择有最少可能值（或者其中一个）的待解方块s。

现在我们可以实现search函数了：

def search(values):
    "Using depth-first search and propagation, try all possible values."
    if values is False:
        return False ## Failed earlier
    if all(len(values[s]) == 1 for s in squares):
        return values ## Solved!
    ## Chose the unfilled square s with the fewest possibilities
    _,s = min((len(values[s]), s) for s in squares if len(values[s]) > 1)
    return some(search(assign(values.copy(), s, d))
                for d in values[s])

终于完成了！我们现在可以在理论上解决所有的数独迷题。在实际应用中，hard puzzles上的95个难题，我们的程序以每秒钟8个难题的速度解决；easy puzzle的容易题则是每秒钟30个。（假如我为了执行效率改写程序，这个程序可以快10倍，但是代码长度会增长到2到5倍。）然而，是否有可能存在这样一个迷题，我们的程序会用极长的时间去解它；我想这是不存在的。在零点几秒的时间里，这个程序可以解决全空数独迷题（81个未解方块），以及我在hardest sudoku上看到的五个迷题。特别的是，有一篇新文章描述了芬兰数学家Arto Inkala所说的“史上最难的数独迷题”；我的程序只用了0.013秒就解出来了（大约超过300次尝试）。

结论
你可以在这里看到Python源代码（100行），或者95个hard puzzle迷题的输出结果（1140行）。

译后记
读Peter Norvig的文章和代码，有如沐春风的感觉。另一篇”Spell Corrector”也值得一读，在Norvig的网站上，已有中文翻译。
Peter Norvig，现在是Google的技术主管。在UCB期间，合著影响广泛的“Modern AI”。

Read more ... ...

Solving Every Sudoku Puzzles: Part 1

数独(Sudoku)迷题的通用解法: Part 1

原文 Solving Every Sudoku Puzzle
Peter Norvig

在这篇文章里，我将提出一种数独(Sudoku)迷题的通用解法。它看上去相当简单（大约100行代码），并采用了两个思想：约束传递（constraint propagation）和搜索（search）。

迷题设定

首先，我们要确定问题符号。在查阅了一些数独(Sudoku)网站之后，我发现在符号上没有统一的规定，但大多数人喜欢行用A-I标记，列用1-9标记，将9个方块(squares)的集合（是这个方块所处的一行，一列，和3x3方块组成的区（box））称为一个unit，并将其他在同一个unit里面的方块称为这个方块的peers。所以，可以这样来命名方块：

 A1 A2 A3 A4 A5 A6 A7 A8 A9
 B1 B2 B3 B4 B5 B6 B7 B8 B9
 C1 C2 C3 C4 C5 C6 C7 C8 C9
---------+---------+---------
 D1 D2 D3 D4 D5 D6 D7 D8 D9
 E1 E2 E3 E4 E5 E6 E7 E8 E9
 F1 F2 F3 F4 F5 F6 F7 F8 F9
---------+---------+---------
 G1 G2 G3 G4 G5 G6 G7 G8 G9
 H1 H2 H3 H4 H5 H6 H7 H8 H9
 I1 I2 I3 I4 I5 I6 I7 I8 I9

（译者注：在这篇文章里，将数独里的最小单位称为方块square，3x3方块称为区box，9x9方块称为数独盘面Sudoku board,，游戏规则和中文命名看这里）。

我们可以在编程语言Python里实现这些：

def cross(A, B):
    return [a+b for a in A for b in B]

rows = 'ABCDEFGHI'
cols = '123456789'
digits   = '123456789'
squares  = cross(rows, cols)
unitlist = ([cross(rows, c) for c in cols] +
            [cross(r, cols) for r in rows] +
            [cross(rs, cs) for rs in ('ABC','DEF','GHI') for cs in ('123','456','789')])
units = dict((s, [u for u in unitlist if s in u])
             for s in squares)
peers = dict((s, set(s2 for u in units[s] for s2 in u if s2 != s))
             for s in squares)

（作者注：以上代码用到了generator expressions，需要Python2.4以上的版本）。
现在，A1的units和peers可以定义为：

>>> units['A1']
[['A1', 'B1', 'C1', 'D1', 'E1', 'F1', 'G1', 'H1', 'I1'],
 ['A1', 'A2', 'A3', 'A4', 'A5', 'A6', 'A7', 'A8', 'A9'],
 ['A1', 'A2', 'A3', 'B1', 'B2', 'B3', 'C1', 'C2', 'C3']]
>>> peers['A1']
set(['A2', 'A3', 'A4', 'A5', 'A6', 'A7', 'A8', 'A9',
     'B1', 'C1', 'D1', 'E1', 'F1', 'G1', 'H1', 'I1',
     'B2', 'B3', 'C2', 'C3'])

简而言之，方块A1有三个unit: 列1，行A，以及左上角区。方块A1有20个peers（行上8个，列上8个，区里且不包括行和列上4个）。所有其他的方块有同样数目的unit和peer。

接下来是设计数独盘面。像符号一样，在初始盘面上也没有统一的规定，但是一个较好的折中（common denominator）是一串字符，如1-9代表数字，0或句号或破折号代表空方块（但是在一些符号命名中，破折号被用于分隔区，而不是标记空方块）。空格以及其它的字符要被忽略。我们将这样的一个字符串称为grid。

我们想要向一个容易操作的数据结构读入这样一个字符串。有人可能认为一个9x9的矩阵是正确的数据结构。由于我们决定用像A1这样字符来标记方块，而不是用[0,0]，所以我们如果希望改用二维数组的话只能改符号标记了。并且，Python不直接支持二维数组（它支持数组的数组），但是它却支持像带有键A1那样的字典dictionaries（hash tables），所以我们将用dict来代表数独盘面。这个dict的键是方块的序号（如A1），而每一个键对应的值是这个方块的所有可能值。而如何来表示这样一个数字集合呢？我们可以用Python内置的set或list，但是我选择这些数字的字符串（我们将在下面看到为什么要这么做）。例如，表示将7填入方块A1，而A2为空，我们可以用一个dict[‘A1’ : 7, ‘A2’ : ‘123456789’, …]。

以下是将一个grid读入dict的程序:

def parse_grid(grid):
    "Given a string of 81 digits (or . or 0 or -), return a dict of {cell:values}"
    grid = [c for c in grid if c in '0.-123456789']
    values = dict((s, digits) for s in squares) ## To start, every square can be any digit
    for s,d in zip(squares, grid):
        if d in digits and not assign(values, s, d):
     return False
    return values

约束传递
函数parse_grid调用assign(values,s,d)来将数字d赋予方块s的值。所以，我们希望以values[s]=d来结束，但是我们还是想要用一些其他的方法来改变values。特别是我们希望在s的peers里消除所有可能的d。如果消除d导致其中的一个peer变成一种可能（译者注：即只有一个数字），我们称之为d2，那么我们希望从这个peer的peer消除所有可能的d2，（译者注：这里好像和代码有些出入？，代码里是消除d导致方块s的值变成一种可能），等等。以上被称为约束传递(constraint propagation)：将某个约束置于一个方块之上会引发蝴蝶效应，即将更多的约束递推到了其它方块上。

这里还有第二种约束传递。比如说我们刚刚将6从方块A1的所有可能值中消去。而假如我们看到包含A1的第一个unit里面所有的方块里，只有C1将6作为它的可能值。我们可以将6赋给C1。所以我们需要两个函数：assign(values,s,d) 和eliminate(values,s,d)，它们会互相调用来实现约束传递（我们将这种函数称为mutually recusive）。虽然不是很明显，但是我们有可能会犯一个错误：我们会试图将一个不符合数独规则的数赋给方块。在这种情况下，我们希望assign()和eliminate()返回False来指示错误。在通常情况下，每一个函数都会用传递来的约束稍稍改变一下数值，然后返回给另一个函数。以下是实现代码：

def assign(values, s, d):
    "Eliminate all the other values (except d) from values[s] and propagate."
    if all(eliminate(values, s, d2) for d2 in values[s] if d2 != d):
        return values
    else:
        return False

def eliminate(values, s, d):
    "Eliminate d from values[s]; propagate when values or places <= 2."
    if d not in values[s]:
        return values ## Already eliminated
    values[s] = values[s].replace(d,'')
    if len(values[s]) == 0:
        return False ## Contradiction: removed last value
    elif len(values[s]) == 1:
        ## If there is only one value (d2) left in square, remove it from peers
        d2, = values[s]
        if not all(eliminate(values, s2, d2) for s2 in peers[s]):
            return False
    ## Now check the places where d appears in the units of s
    for u in units[s]:
        dplaces = [s for s in u if d in values[s]]
        if len(dplaces) == 0:
            return False
        elif len(dplaces) == 1:
            # d can only be in one place in unit; assign it there
            if not assign(values, dplaces[0], d):
                return False
    return values

这里有一种有用的设计模式，好像从没有人提过。这个模式是：
如果你有两个mutually-recursive的函数分别影响一个对象的状态，请试着将所有的功能代码移到其中一个函数去。否则，你最后会发现有许多重复代码。

我是在许多年的Lisp编程后发现这个设计模式的，在Lisp里mutually-recursive函数很常见。看看我们怎么样在这个问题里应用这个模式：有人会认为assign()将包含赋值语句values[s]=d，而且会包含传递约束。你可以试着写这样一个函数。我想，你最后会发现你是在重复eliminate()里的代码。所以为了避免绕这么个弯子，我推论assign()函数做的就是消去方块s里除了d以外所有的数字，所以我将所有的功能代码写到了eliminate()里。

在我们探索更远之前，我们需要能够检验一下数独盘面的状态。以下就是printboard()的代码：

def printboard(values):
    "Used for debugging."
    width = 1+max(len(values[s]) for s in squares)
    line = '\n' + '+'.join(['-'*(width*3)]*3)
    for r in rows:
        print ''.join(values[r+c].center(width)+(c in '36' and '' or '')
                      for c in cols) + (r in 'CF' and line or '')
    print

现在我们可以开始解题了。我选了easy puzzles上的第一个问题，试了一下：

>>> grid = """
003020600
900305001
001806400
008102900
700000008
006708200
002609500
800203009
005010300"""

>>> printboard(parse_grid(grid))
4 8 3 9 2 1 6 5 7
9 6 7 3 4 5 8 2 1
2 5 1 8 7 6 4 9 3
------+------+------
5 4 8 1 3 2 9 7 6
7 2 9 5 6 4 1 3 8
1 3 6 7 9 8 2 4 5
------+------+------
3 7 2 6 8 9 5 1 4
8 1 4 2 5 3 7 6 9
6 9 5 4 1 7 3 8 2

在这个例子里，这个数独迷题完全被我们的约束传递解开了！只是通过赋予32个方块值，我们简单的约束传递规则就把剩下来的所有方块都填满了。但是，不是所有的题都是这么容易。接下来是hard puzzles里的第一个问题:

>>> grid = '4.....8.5.3..........7......2.....6.....8.4......1.......6.3.7.5..2.....1.4......'

>>> printboard(parse_grid(grid))
   4      1679   12679    139     2369    269      8      1239     5
 26789     3    1256789  14589   24569   245689  12679    1249   124679
  2689   15689   125689    7     234569  245689  12369   12349   123469
------------------------+------------------------+------------------------
  3789     2     15789    3459   34579    4579   13579     6     13789
  3679   15679   15679    359      8     25679     4     12359   12379
 36789     4     56789    359      1     25679   23579   23589   23789
------------------------+------------------------+------------------------
  289      89     289      6      459      3      1259     7     12489
   5      6789     3       2      479      1       69     489     4689
   1      6789     4      589     579     5789   23569   23589   23689

在这个例子里，我们离解出这个问题还差得很远。我们开始只有17个方块填了数字（这被认为是最少的可以达到唯一解的数目），在约束传递之后，只有3个方块被解了出来（虽然所有的方块都被消去了一些可能值）。

我们接下来要怎么做呢？我们可以尝试一些更加复杂的约束传递技巧，就像这里说的。比如说naked pairs技巧寻找在同一个unit里的两个方块，它们有两个相同的可能值。假设A1和A4都有2和6的可能值。我们可以推论2和6一定在A1和A4里（虽然我们不知道究竟哪个在哪个里面），而且我们可以将2和6从这个A行unit里面的其它方块里面消去。我们可以仅仅在在代码里加上几行，比如elif len(values[s]) == 2来实现这个功能。

类似的代码技巧是可行的，但是会导致代码量的膨胀（大概有二三十种技巧），而且我们永远不会知道我们是否能依靠这些技巧解出所有迷题。

Go to Part 2

Read more ... ...

Tips to make your blog beautiful

A website have intuitive HTML tutotial. This site provide HTML toolbar plugin for VIM.

Add <p> before every paragraph and </p> can be omitted.

Wood Xiao told me that use Shift-Enter instead of Enter when we need a return in blog.

I follow this post to represent expandable post summary.

I follow this postto do code prettify(google code prettify)

I will follow this post, or this Chinese post to insert LaTeX to my blog.



Python code

def cross(A, B):
    return [a+b for a in A for b in B]

Java code

class Voila {
public:
  // Voila
  static const string VOILA = "Voila";

  // will not interfere with embedded tags.
}

Read more ... ...