2009年1月20日星期二

Simulating N-Body Problem


理解多体问题的数值模拟

多体问题
物理上,多体问题(N-Body Problem)是指:3维空间中给n个质点在某一时刻的位置和速度,已知他们的质量,假定他们彼此之间只受万有引力作用,求他们后续任意时刻的位置和速度。数学上,多体问题是一个高阶的常微分方程初值问题,具体说来:
设位置函数向量x=(x1,...,xn),,质量mi
满足方程

其分量形式为

和初值条件

求x(t)
当n=2时(又称二体问题),有解析解,解曲线是平面上的二次曲线(《常微分方程教程》P132)。
n>2时目前没有解析解,虽然理论上得到了很多很漂亮的定性分析结果,但是实践中还要归诸数值模拟,也就是用差分代替微分。

多体问题的数值模拟
Part1:差分格式
参考http://en.wikipedia.org/wiki/Verlet_integration
数值模拟最朴素的想法是欧拉迭代,在每一个时刻,根据当前位置计算加速度,更新下一时刻的速度,根据当前速度更新下一时刻位置:(设时间步长为h)
a(t) = f(x(t))
v(t + h) = v(t) + a(t)h
x(t + h) = x(t) + v(t)h
这个方法误差很大,即使假设速度v没有产生误差,也很容易分析出(通过Taylor展开到二阶项),迭代一步产生的误差为O(h2)。

经典力学中的数值模拟常用的是Verlet方法,它利用两侧Taylor展开消去速度项,构造误差阶数更小的差分格式:
首先Taylor展开到4阶

注意到

并将-h带入上面的h,得

两式相加,消去1阶项,顺带把3阶项也消掉了(赚到了,呵呵)。得到:

定义差分格式为

显然一步误差O(h4)
然后来分析误差扩散。
(注:通常Big O Notation是不能做减法的,这里我们不加说明地假设它可以做减法,使用Lagrange余项的Taylor展开也许能解释清楚,我们有深究它。)
来归纳证明,k步误差有如下二次扩散关系

k=1、2时不难直接验证。假设直到k-1都成立,来看k时

其中假设f连续可微,x取值在一个紧集上,就有f是Lipschitz连续的,从而两个h2f项做差是高阶小,直接扔掉。从而上式成立。
因此,若模拟总时间T=kh固定,则总误差O(h2)


Part2:受力计算的分治算法
参考http://www.cs.berkeley.edu/~demmel/cs267/lecture26/lecture26.html
计算每个质点受的合力,按照定义,是其他质点对它的万有引力的矢量和。n个点,每个点计算(n-1)个分力,这样迭代一步的复杂度是O(n2),不能忍。
Barnes-Hut算法使用分治策略,将复杂度降低到O(nlogn),它的思想是,将较远的若干个点当做一个点,一起计算来自他们的分力。举例来说,假设x,y到点集A的质心的距离很远,A包含a个点,则按照定义计算x,y受A的力需要计算2a个分力,而若把A当做一个对象处理,预先计算它的质心和总质量,则完成同样的工作,仅需要a次计算。
算法用到了四叉树的结构(这是对于2维情形,3维使用八叉树)。
递归地将矩形分成四等分,落入第i象限的点属于当前节点的第i子节点,直到当前节点只有一个节点为止,依此从平面上的点集构造四叉树,如图——

通过适当选取划分轴,可以保证四叉树均衡,也称自适应四叉树,这一步复杂度O(nlogn)
得到四叉树后,后续遍历一次得到每个节点的质心和总质量。最后,计算每个质点的受力:
从根节点开始考虑每个子树,如果质点到该子树的根节点的质心足够远,则累加质点受到该子树(作为整体对应的质心和总质量)的万有引力,否则,质点受到该子树的合力等于该子树的每个子树的引力矢量和,而计算子树的引力又递归调用当前算法。由四叉树对应的几何空间划分,很容易看出,计算一个质点的受力,要遍历树上的节点个数不超过树的高度。若四叉树均衡,则复杂度O(nlogn),从而总复杂度也如此。

总结:
1、利用Taylor展开构造有效差分格式可以减少计算误差
2、利用分治思想改进算法可以降低计算复杂度

题外话:
计算数学看算法就看三件事:收敛性、复杂性、稳定性。这个例子中讨论了前两个,第三个稳定性主要是针对良好的问题的,也就是说,如果真实解对初值的依赖是连续的,则要求算法的数值解对初值的依赖也是连续的。但是多体问题是一个典型的混沌系统,问题本身是病态的,讨论它的算法稳定性也就意义不大了。
一直以来,我就特别看好物理的sense好的人,像上面这个分治算法的设计,物理背景就很强(berkeley的讲义上用仙女座星系做四叉树构造思想的导言),或者说,对问题域有透彻了解的人,往往能发明出新的、好的方法。

2009年1月13日星期二

Keypoints for Basic Functional Analysis

Keypoint for Basic Functional Analysis
基础泛函分析知识点整理
度量空间

聚点:x∈G',若x的任何去心邻域交G非空

闭集的等价定义:
1、余集是开集
2、闭包为自身
3、若其中一序列收敛,则极限点属于集合

闭包的等价定义
1、G的闭包=G∪G'
2、G的闭包是包含G的最小闭集

可分空间的子空间可分
对X'可分,则X可分
C[0,1]可分
l∞不可分

完备子集闭
完备空间的闭子空间完备
自反,则完备
离散度量空间总完备
lp、l∞、c0、c、C[a,b]完备

M是紧子集,若M的任何序列都有收敛到M中的子列
紧空间完备
紧子集是有界闭集
紧集中的闭集是紧的
点到紧集的距离总可达到

压缩映射:T:X->X成为压缩的,若存在0<a<1,d(Tx,Ty)<=ad(x,y)对于任意x,y∈X
不动点:x为T不动点,若Tx=x
不动点定理:X完备,存在m>0使得Tm压缩,则T有唯一不动点

连续映射:
开集的原像是开集、闭集的原像是闭集
把紧集映成紧集
与极限号可交换

线性空间

G是Hamel基,若G无关且spanG=X
有限维线性空间代数自反

实线性空间X,Z为子空间,p为X上次线性泛函,f属于Z*,f(x)<=p(x)对于任意x∈Z,则存在g属于X*使得g|Z=f,g(x)<=p(x)对于任意x∈X
线性空间X,Z为子空间,p为X上半范,|f(x)|<=p(x)对于任意x∈Z,则存在g属于X*使得g|Z=f,|g(x)|<=p(x)对于任意x∈X

赋范空间

度量具有平移不变性、比例缩放、球可线性运算

(en)为Schauder基,若对于任意x∈X,存在唯一λn使得x=∑λnen
X有Schauder基,则X可分
Schauder基无关,则dimX=∞

Reiz引理:X是赋范空间,Y是Z真子空间,Y是闭集。对于任意0<θ<1,存在z∈Z,||z||=1,使得ρ(z,Y)>θ

X是有限维赋范空间,(ei)是X的基,则存在a>0,使得a||x||1<=||x||对于任意x成立
有限维赋范空间是Banach空间
有限维赋范空间中有界闭集是紧集
赋范空间X,X有限维当且仅当单位闭球是紧的
赋范空间中,点到有限维子空间的距离总可达到

X严格凸,点到闭子空间的距离至多可达一处

X,Y赋范空间,T:X->Y线性,则三者等价:T连续,T在0处连续,T有界算子

Z为X子空间,f∈Z',则存在g∈X',g|Z=f,|g|=|f|
对于非零元x,存在f∈X',|f|=1,f(x)=|x|
|x|=max{f(x)|f∈X',|f|=1}
Y真闭,x=X-Y,则存在f∈X',|f|=1,f|Y=0,f(x)=ρ(x,Y)>0

T属于B(X,Y),定义共轭算子T×:Y'->X',T×f定义为fT,f∈Y'

自反,则完备
H自反,lp自反
l1、l∞、c0、C[a,b]不自反

Banach空间

赋范空间完备,当且仅当任何绝对收敛的级数都收敛
X,Y赋范空间,B(X,Y)完备当且仅当Y完备

X完备,Tn∈B(X,Y):
对于任意x∈X,Tnx有界,则Tn有界
Tn弱收敛到T,则T有界
T∈B(X,Y),Tn强收敛到T,当且仅当Tn有界且Tn在X的一个完全集上处处收敛到T

X,Y完备:
T∈B(X,Y),若T满,则T为开映射,若T双,则T逆连续
D(T)闭,T闭,则T∈B(X,Y)
X,Y赋范空间,T∈B(D(T),Y):
D(T)闭,则T闭
T闭,Y完备,则D(T)闭

Hilbert空间

范数由内积给出,当且仅当平行四边形等式成立|x+y|2+|x-y|2=2|x|2+2|y|2

内积空间中,点到完备凸集的距离唯一达到,特别到闭子空间对
内积空间中,点到子空间的距离可达到

任何集合的正交补都是闭子空间
X=MN,若任意x∈X,存在唯一m∈M,n∈N,使得x=m+n
X为Hilbert空间,M为闭子空间,则X=MM正交补
M完全当且仅当M正交补为{0}

Bessal不等式:(ei)为标准正交列,则∑|<x,ei>|2<=||x||2

完全标准正交集,可数则Schauder基,有限则Hamel基
Hilbert空间总有完全标准正交基
伴随算子:H1,H2为Hilbert空间,T∈B(H1,H2),T*∈B(H2,H1)使得<T*y,x>=<y,Tx>

点到集合的距离
有限维内积



2009年1月12日星期一

广播加密算法调研报告



广播加密算法调研报告

一、广播加密问题
1、问题的提法
一个系统有一个广播源和n个用户,订阅者S为用户集的一个子集,大小为k。广播源将消息m加密后广播给所有用户,使得S中的用户解密得到m,非S中的用户
无法解密得到m。其中如何加密解密,并使得密钥和密文长度尽量短,就是广播加密问题。相当于多把钥匙一把锁,锁可以制定,想让谁能开谁就能开,而且想让谁能开谁才能开。
注意到,不可以给每个订阅者单独加密,这样虽然可以把广播加密问题归结为已经解决的点对点加密问题,但是实际中不能接受:应为那样广播源加密的计算开销将随着k的增大而线性增大。这是广播加密问题和传统的加密问题的区别。
广播加密问题的形式化提法:
通过
Setup(用户总数 n)函数得到一个广播源的公钥PK,每个用户i的私钥di,i=1,..,n。广播源针对订阅者集合S加密消息m,即执行加密函数Encrypt(广播明文m, 订阅者集合S, 公钥 PK)得到密文Hdr,通过网络广播Hdr给所有用户。所有用户都有相同的解密函数Decrypt(密文Hdr,私钥d)用户i接收到Hdr后计算Decrypt(Hdr, di)Decrypt(Hdr, di)=m当且仅当i属于S。注意到,初始化公钥和私钥的算法与订阅者S无关,实际中订阅者会经常变化。我们希望开发出其中Setup、Encrypt、Decrypt的算法,满足一系列安全要求,并且PK、di关于n不太长,给定m后Hdr不太长。
关于广播正文m的一点注解:实际中广播内容是一段比较长的视频流stream,但是用此非对称加密算法加密stream可能无法满足实时性需要,因此,我们用计算复杂度比较低的某种对称加密算法加密stream,设对称密钥为K,然后用广播加密算法(非对称加密)加密此对称密钥K,然后广播非对称加密后的K和用K对称加密过的stream即可满足实时性的需要。因此我们下面讨论的广播正文m,实际就是这里的对称密钥K,因此,m只需要是某个随机的整数即可,与待广播的视频流stream无关。

2、广播加密问题的若干安全要求:
防破解:防止非S中的用户通过Hdr得到m。
盗版追踪:如果S中的某个i将di出卖给非S中的j,如何检测到并确定i的行为。
防恶意插播:如果另一个非法广播源在广播信道中广播内容,接收端如何识别,也就是如何验证广播源的身份。 

3、接收端无状态
[todo]

二、 广播加密算法
1、多项式插值
91Berkovits的《How to broadcast a secret1》提出这样一个基于拉格朗日多项式插值的广播加密算法,它利用了这样一个事实:一个n阶多项式能被平面上(横坐标互不相同的)n+1个点确定。
算法构造的思路如下(a special but typical case)
Setup(用户总数 n)
随机生成平面上
横坐标互不相同的n个点,分配(xi,yi)作为用户i的私钥di。不妨设所有xi都不为0。
不用公钥。
Encrypt(广播明文m, 订阅者集合S),S的大小为k
计算一个k阶多项式P,过下列k+1个点:(0,m), (xi,yi) ,i属于S
然后广播密文Hdr为一个长度为k的串/数组:(xi,yi),i不属于U,且(xi,yi)在P上,不妨设所有xi都不为0。
Decrypt(密文Hdr,私钥di)
如果用户i属于S,则i知道k+1个在P上的点(k个来自Hdr,还有一个是自己的私钥),从而确定P并能够求出P在0点的值,也就是密文m
否则i只知道k个在P上的点,自己的私钥(xi,yi)不在P上,无法求出P在0点的值。

实际加密的时候还要引入j个扰乱的点,可能是为了安全,为了计算方便,作者给出了矩阵描述的等价算法。


2、子集覆盖
3、配对
05年
Boneh等人的《Collusion Resistant Broadcast Encryption With Short Ciphertexts and Private Keys2》利用配对技术构造了这样一个广播加密算法,它利用了双线性DH假设。具体来说,给定误差限,给定p阶双线性群G和它上面的双线性映射e(相关定义和性质见附录),若h,g随机取值于G中 随机取值于中,则不存在这样一个(确定性的?且多项式时间的?)算法A使得
,其中
注意到,如果能在G中从g和计算出,则显然可计算出因为其中已知。在G中从g和计算出就是经典的离散对数计算,因此双线性DH假设比离散对数假设更强,如果离散对数多项式可解则双线性DH假设不再成立!(介于二者难度之间的还有一个DH假设,见附录)
算法构造的思路如下
(a special but typical case):
Setup(用户总数 n)
随机取值于中,令
公钥PK:
私钥di:
Encrypt(广播明文m, 订阅者集合S)
t随机取值于中,不妨设,广播密文为一个二元数组Hdr=(C0,C1),其中


Decrypt(密文Hdr,私钥di)


三、研究者
1、Naor
2、Boneh
四、应用
1、DRM
2、CA for IPTV
附、Mathematical Background


http://www.springerlink.com/index/XXBDHG4QH59RGCH4.pdf

http://crypto.stanford.edu/~dabo/papers/broadcast.pdf

2009年1月9日星期五

LLE Algorithm and Implement in Matlab

《数学建模案例分析》的大作业用LLE算法,但是原作者网站上提供的源代码有些问题,主要是因为不同版本的Matlab,内置函数eigs返回的特征向量的顺序不同:老版本对应的特征值是升序,而新版本的是降序。

在这个问题中,0总是特征值,对应的特征向量为(1,1,…,1),这不是我们要的(如果把它放进来,则用它给swiss_roll数据集降维总得到一条粗直线。)

如果你用的是Matlab6.5,则应该把line 65

Y = Y(:,2:d+1)'*sqrt(N);

改为

Y = Y(:,n-d:n-1)'*sqrt(N);

打了这个补丁之后,lle.m就可以正常工作了!

 

我认为LLE算法很漂亮,下面过一下它的设计与实现~原作者的paperlle做关键词即可在Google上搜到。

 

问题的提法:

ND维列向量X1,…,XN,希望通过映射得到Nd(<D)维列向量Y1,…,YN,要求保持邻域关系:原来离的近的点,映射过来也近。

LLE算法的思想:如果原像Xi能够表示成他邻域内点的线性组合,则像Yi也应该用相同的组合系数表示成对应像点的线性组合。即:局部线性嵌入。

算法骨架

1、  求近邻:计算每个Xi的邻居集合Si

直接用2范数下的K近邻,其中K作为算法的参数

2、  求权:改变组合系数W,极小化局部线性表出原像的方差

WN×N矩阵,Ei=Xi-{jSi}{WijXj}D维残差向量,极小化E(W)={i}{|Ei|^2},满足如下约束:归一化——W行和为1、局部性——对于任意iWij=0j不属于Si

3、  求像:改变Y1,…,YN,极小化用W重构Y1,…,YN的方差

Ei=Yi-{jSi}{WijYj}d维残差向量,极小化E(Y)={i}{|Ei|^2}

算法实现细节

1、  求近邻

XDN列的矩阵,X的第j列为输入向量Xj

向量XiXj的距离为sqrt(|Xi-Xj|^2),得到距离方阵distance,然后每行排序,取前k个,对应原下标就是近邻点的编号。在Matlab中实现如下

X2 = sum(X.^2,1);

distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;

[sorted,index] = sort(distance);

neighborhood = index(2:(1+K),:);

这里充分利用了矩阵运算技巧避免了编程中使用循环,不仅使代码紧凑,而且Matlab矩阵运算底层优化了性能。下面逐一解释——

a)       距离方阵

首先注意到这样一个事实|Xi-Xj|^2=|Xi|^2+|Xj|^2-2<Xi, Xj>,其中< , >为内积。

X2 = sum(X.^2,1)1N列的矩阵,第j个元素为Xj列的平方和,也就是|Xj|^2

repmat是平铺函数,repmat(X2,N,1)得到N×N方阵,每行都是X2,同理repmat(X2',1,N)每列都是X2''是转置的意思。X'*X得到N×N方阵,(i,j)位置元素为<Xi, Xj>。因此distance就是我们要的距离方阵的平方。

b)      省去开平凡

因为我们只需要K近邻,而不需要具体的距离值,而开平方是严格单调函数,保序,因此我们可以省去开平方,节省了计算量。

c)      反向索引

sort给矩阵每列排序,并返回置换前后下标的映射关系indexindex(i,j)distancej列第i小的元素的下标。我们要每列前K小的下标,也就是每个XiK近邻。取2:(1+K)是因为Xi到自己的距离总是0最小,排除掉。

       任意两个点至少做一次内积O(D),共O(N^2)个点对,故复杂度:O(DN^2)

2、  求权

首先,因为对于任意iWij=0j不属于Si,故W只需要存KN列即可。

其次,易见E(W)极小当且仅当每一求和项极小,因此我们依次计算W的每一列。固定列i,记x=Xiw=Wi列,ŋj=Xj,极小化|x-{j=1..K}{wjηj}|^2,满足归一化约束∑{ j=1..k }{wj}=1。用矩阵语言描述:记B=( ŋ1-x,…, ηk-x)D×K矩阵,G=B'BK×K方阵(讲义中称之为Gram方阵,半正定,在摄动意义下总可以假设它非奇异),e=(1,…,1)'K维单位列向量,则问题化为——

min |Bw|^2也就是min w'Gw(二次型)

s.t. e'w=1

用拉格朗日乘数法求此条件极值:做辅助函数F(w,λ)= w'Gw-λ(e'w -1)

对每个wj求偏导数令为0Gw=λe,反解出w=G^{-1}λe,代入到归一化约束得

λ=(e'G^{-1}e)^{-1},即最优解w=(e'G^{-1}e)^{-1} G^{-1}e

实际操作时,我们先解线性方程组Gw=e,然后再将解向量w归一化,易见得到的就是上述最优解。

Matlab中如下实现:

W = zeros(K,N);

for ii=1:N

   z = X(:,neighborhood(:,ii))-repmat(X(:,ii),1,K); % 计算B

   C = z'*z;                               %  计算G,复杂度O(DK^2)

   C = C + eye(K,K)*tol*trace(C);             % 必要时摄动

   W(:,ii) = C\ones(K,1);                  % 解方程Gw=e,高斯消去复杂度O(K^3)

   W(:,ii) = W(:,ii)/sum(W(:,ii));               % 解向量w归一化

end;

故总复杂度O(DNK^3)

3、  求像

将上一步得到的W视为N×N方阵,Yd×N矩阵,YjYj为降维映射下Xi的像。min {i}{|Yi-{j}{WijYj}|^2}

注意到Yi-{j}{WijYj}={j}{ Wij (Yi-Yj)},因此若Y为最优解,则所有Yi平移任一相同的向量也是最优解,为了定解,我们不妨假设∑{j}{Yj}=0。事实上,若∑{j}{Yj}=Z,则有∑{j}{Yj-Z/N}=0

此外,Y=0为平凡最优解,为了避免这种退化情形,我们不妨假设∑{j}{YjYj'}/N=I,即∑{j}{YajYbj}=Nδ(a,b),即Yd个行向量,都在半径为sqrt(N)N维单位球上,且两两正交。

验证:∑{i}{|Yi-{j}{WijYj}|^2}={i,j}{MijYi'Yj},其中M=(Mij)=(I-W)'(I-W)

={i}{Yi'Yi-2{j}{WijYi'Yj}+({j}{WijYj'})({k}{WikYk})}

={i}{Yi'Yi}

-2{i,j}{WijYi'Yj}

+{i,j,k}{WijWikYj'Yk}

={i,j}{δij Yi'Yi}

-{i,j}{WijYi'Yj}

-{i,j}{WjiYi'Yj}

+{j,k}{{i}{WijWik}Yj'Yk}

={i,j}{(δij-Wij-Wji+{k}{WkiWkj})Yi'Yj}

=

而在单位球上极小化二次型∑{i,j}{MijYi'Yj}当且仅当Y每行是M最小的d个特征值对应的特征向量,排除0对应的特征向量。(d=1的情形为Rayleigh-Ritz定理。)

这步复杂度取决于计算矩阵(部分)特征值/向量算法的复杂度。

2009年1月8日星期四

Playing with Nonlinear Dimensionality Reduction

这学期有一门课程《数学建模案例分析》,讲到了两种经典的非线性降维(Nonlinear Dimensionality Reduction)方法:lleisomap。问题是别人做过的,算法实现也基本都是现成的,我只是拿来"玩一玩"。实验有很多环节,最有趣的一个环节, 是给你698张人脸的图像(64×64灰度),通过isomap降维方法将每张脸当做一个点映到二维平面上,使得横坐标恰好反映人脸左右看的程度,纵坐标反映人脸上下看的程度。

如果你也对这个实验感兴趣,就往下读吧,很简单的~

实验环境Matlab6.5

 

实验步骤

步骤一:准备数据集和工具包

下载

人脸数据集:http://waldron.stanford.edu/~isomap/face_data.mat.Z

并解压缩

 

下载

isomap算法实现:http://waldron.stanford.edu/~isomap/code/

的所有代码

 

步骤二:

准备图片标记的人脸序号集:一共有698张人脸,都画在平面上太拥挤了,所以选了30个人脸(存入posesSelect.matks向量),选取的准则是:30个人脸的姿态尽量不同,也就是希望画在平面上尽量分散。事实上,face_data.mat数据集中,poses是一个2698列的矩阵,第j列就是第j张人脸的客观姿态。

绘制客观姿态分布图:

load face_data

load posesSelect

showFacesOnR2(images,poses,ks)

步骤三:降维

Isomap算法将4096维的人脸数据images降维到2维,并绘制在平面上

load face_data

load posesSelect

D=L2_distance(images,images,1);

options.dims = [2];

[Y, R, E] = IsomapII(D, 'k', 7, options);

showFacesOnR2(images,Y.coords{1},ks);

D是一个距离矩阵,ij列值表示人脸i和人脸j的距离,这里把一个人脸图像数据当做一个向量,使用2范数定义距离。

IsomapII是高性能算法,先把Dk=7近邻打成稀疏矩阵,然后用基于斐波那契堆的Dijkstra算法计算最短路,Dijkstra算法用C实现使用并且编译成了.dll文件为了提高效率。计算结果对我们有用的是Y.coords{1},它保存了降维后的结果,是2698列的矩阵。

观察计算结果发现,以中间那个正的人脸为中心,他左边的都在向左看,而且越是靠左的转动越明显。同理,他右面的都在向右看、上面的都在向下看、下面的在向上看。与客观姿态分布基本吻合。

 

实验细节:

showFacesOnR2.m

%把头像和姿态坐标画在平面上

function showFacesOnR2(images,poses,ks);

%normalize into 1:1

poses(1,:)=poses(1,:)/range(poses(1,:));

poses(2,:)=poses(2,:)/range(poses(2,:));

%draw all points

scatter(poses(1,:),poses(2,:),12,'o','filled');

xlabel('left-right pose');

ylabel('up-down pose');

hold on

%draw selected points

scatter(poses(1,ks),poses(2,ks),24,'ro');

hold on

%draw images on selected points

scale = 0.001;

x=zeros(64,64);

for p=1:size(ks,2)

    k=ks(p);

    for i=1:64

        x(:,i)=images((i-1)*64+1:i*64,k);

    end

    xc=poses(1,k);

    yc=poses(2,k);

    imshow(xc:scale:xc+64*scale,yc:-scale:yc-64*scale,x);

    hold on

end

return

 

 

FoxPro8: update with join

元旦那天,一个高中同学找来,要我改一个FoxPro8的程序,生成满足一定要求的抽样,听起来不复杂。我没碰过那么古老的数据库软件,问,FoxPro支持SQL么,他说支持一部分,我说那好吧我试试看。
最关键的一步是让FoxPro8的update语句支持跨表(只需要自然连接),这是SQL标准,但FoxPro8还不支持。Google了一通,找到了解决方案,虽然不是很理解,但是照猫画虎勉强拿来用了。


SET RELATION TO k INTO R IN S
这里k是表R和S的公共字段,都在k上建了索引


REPLACE myfield WITH myvalue IN S FOR FOUND("R")
这样做结果相当于SQL语句
UPDATE S, R
SET S.myfield = myvalue
WHERE S.k = R.k

非标准的SQL语法让我颇费周折……

一个历史太悠久的公司,不敢轻易升级业务系统的软件,导致越往后工作越没效率。
这也是为什么工业革命之后英国一度辉煌,但是后来还是被美国超过的缘故。
这也是Diane解释FreeWheel之所以能快速发展的原因之一:没有历史包袱,敢用最新最好的技术。