1971 年的一天，汉森先生说：“今天我要教你如何求一个数的平方根。 比之前更快的方法长除法我应该教的方式。 我会用一个词来教你。 而且你永远不会忘记。”然后他在黑板上写
猜测
并强调它。 笑声平息后，他解释了将猜测分为原始数字、平均数、重复数。 他大概提到了牛顿-拉夫森.
牛顿求方程的根或零点的方法使用方程的导数来加速这个过程。 10 的平方根看起来像是找到一个 X 的值，它使 X^2 == 10 或 X^2 – 10 == 0。 X^2 – 10 的导数是 2X，寻找 X 的连续值的循环的核心是
Xnew = X - ( (X^2 - 10) / (2X) )
如果你玩它，那正是猜测公式，并像这样收敛，初始猜测为 2：
2 – 初步猜测
3.5 – 应用 x = (x + 10/x)/2
3.17857142857143
3.16231942215088
3.16227766044414
3.16227766016838
3.16227766016838 – 六次迭代得到一个稳定的答案！
这是一个 JSL 程序来演示。
f = Expr( x ^ 2 - V );
step = Eval Expr( x - Expr( Name Expr( f ) ) / Expr( Derivative( Name Expr( f ), x ) ) );
V = 1e2; // find the sqrt of V
x = -8; // initial guess.
xnext = .; // set up loop
n = 0; // counter
While( 1,
xnext = step;
If( Abs( x - xnext ) < Abs( 1e-14 * x ),
Break()
);
x = xnext;
Write( "\!n", n, " ", x );
n = n + 1;
);
Write( "\!ndone, a root of " || Char( Name Expr( f ) ) || " for V = " || Char( v ) || " is ", x );
第一行定义了一个表达式 f，它将帮助找到 V 的平方根。 第二行 step 是一个循环的 Newton-Raphson 迭代； step 出现在 while 循环中，它计算为下一个 X 值。 Derivative 函数在简单多项式上效果很好，在本例中为 2X。 NameExpr() 从 f 获取公式而不对其进行评估，当 EvalExpr() 找到 Expr() 时，Expr() 使 NameExpr() 运行（以获取公式）。 步长设置为 x - f(x)/f'(x)。
在 while(1) 循环中，执行步骤并打印和计数。 测试是检查当前和下一个 X 值的差异；如果在 X 的当前大小下该值非常小，它可能不会更接近并且可能在两个非常接近的数字之间振荡，因此它会跳出 while 循环。
0 -10.25
1 -10.0030487804878
2 -10.0000004646115
3 -10
完成后，对于 V = 100，x ^ 2 - V 的根是 -10
请记住，有正负平方根！ 起始值确定收敛导致哪个根。 从 -8 开始超过 -10.25，然后稳定到 -10。
（为什么它除以导数 f' ？ 直观地说，当直线 f 的斜率接近于零时，下一个点需要更远才能改变 f 的 Y 结果，除以斜率即可。）
上面的 JSL 不能很好地处理不可能的请求。 如果被要求查找 sqrt(-100)，则 JSL 循环，查找它无法查找的复数。
复数
在 JSL 脚本索引中，搜索 complex 会产生两个不同的半生不熟的示例。 一种基于 JSL 类，一种基于命名空间和矩阵。 我认为矩阵版本更接近于一个强大的解决方案，并且可能更快。 课程版本很诱人，但我在使用它时遇到了麻烦（速度、内存问题）。 由于 JMP 没有用于运算符的用户重载，因此两者都不会易于使用。
此 JSL 与之前的 JSL 非常相似。 有一大块代码来处理复数算术，相同的公式，或多或少，现在它也处理 sqrt(-100)。
ns = New Namespace("complex", // adapted from scripting index example
{
make = Function( {a, b}, a || b ),
add = Function( {x, y, t=[0 0], u=[0 0], v=[0 0], w=[0 0]}, // an example has 5 terms
if( Is Number( x ),x = x || 0);
if( Is Number( y ),y = y || 0);
x + y + t + u + v + w;
),
sub = Function( {x, y, t=[0 0], u=[0 0], v=[0 0], w=[0 0]}, // an example has 3 terms. this sub,add trick is good enough for a demo.
if( Is Number( x ),x = x || 0);
if( Is Number( y ),y = y || 0);
x - y - t - u - v - w;
),
mul = Function( {x, y},
if( Is Number( x ),x = x || 0);
if( Is Number( y ),y = y || 0);
local:z = J( 1, 2 );
local:z[1] = x[1] * y[1] - x[2] * y[2];
local:z[2] = x[1] * y[2] + x[2] * y[1];
local:z;
),
div = Function( {x, y},
if( Is Number( x ),x = x || 0);
if( Is Number( y ),y = y || 0);
local:z = J( 1, 2 );
local:d = (y[1] ^ 2 + y[2] ^ 2);
local:z[1] = (x[1] * y[1] + x[2] * y[2]) / local:d;
local:z[2] = (x[2] * y[1] - x[1] * y[2]) / local:d;
local:z;
),
// https://math.stackexchange.com/questions/1397437/formula-for-raising-a-complex-number-to-a-power
// If a,b are real, then a+ib=r(cosθ+i*sinθ) where r^2=a^2+b^2 and tanθ=b/a, and (a + i*b)^N = r^N(cos(Nθ) + i*sin(Nθ))
pow = function( {c, n},
if( Is Number( c ),c = c || 0);
local:a = c[1];
local:b = c[2];
local:theta = ATan( local:b, local:a );
local:r = Sqrt( local:a ^ 2 + local:b ^ 2 );
local:r ^ n * (Cos( n * local:theta ) || Sin( n * local:theta ));
),
abs = function( {x},
Sqrt( x[1] ^ 2 + x[2] ^ 2 )
),
toChar = Function( {x,w=20,d=14}, // rounding (to 14 places) is not a good idea if the numbers are small ( 1e-9 will lose lots of precision)
Char( round(x[1],d),w ) || "+" || Char( round(x[2],d),w ) || "i"
),
convertExpression = function( {e},
substitute( Name Expr( e ),
Expr( Add() ), Expr( complex:add() ),
Expr( Subtract() ), Expr( complex:sub() ),
Expr( Multiply() ), Expr( complex:mul() ),
Expr( Divide() ), Expr( complex:div() ),
Expr( Power() ), Expr( complex:pow() )
)
)
}
);上面为 x^2 或 x^3 等添加了一个 pow() 函数，其中 x 是一个复数，一个 abs() 函数，它可能应该被称为幅度（）但与非复数版本匹配为 abs，有些奇怪- 寻找 add() 和 sub() 的额外参数，以及一个 convertExpression() 函数。 向 stackexchange/Michael Hardy 大喊大叫，以获得易于掌握的 pow() 函数！
convertExpression 函数将带有 +、-、*、/ 和 ^ 运算符的 JSL 表达式处理为调用复杂命名空间函数的等效表达式。 JSL 二元运算符非常积极地将 ((a+b)+(c+d)) 转换为 add(a,b,c,d) 这不是二元函数......它有四个参数。 这就是 complex:add 和 :sub 函数有额外选项参数的原因。
JSL 的下一位将循环移动到一个函数中，预测未来的发展：
findzero = function( {x},
{xnext, n},
xnext = .; // set up loop
n = 0; // counter
While( n < 200,
xnext = step;
If( complex:Abs( complex:sub( x, xnext ) ) < 1e-14 * complex:abs( x ),
Break()
);
x = xnext;
n = n + 1;
);
If( n == 200, x = [. .] );
Eval List( {x, n} );
);findzero 采用起点 x 并返回一个包含迭代解 x 和迭代计数 n 的列表。
以及调用函数的代码……复数的公式和导数是相同的。
f = Expr( x ^ 2 - V );
step = Eval Expr( x - Expr( Name Expr( f ) ) / Expr( Derivative( Name Expr( f ), x ) ) );
step = complex:convertExpression( Name Expr( step ) );
V = complex:make( -100, 0 ); // find the sqrt of V==-100
x = complex:make( 0, 1 ); // initial guess. Does not like 0i for a start value.
{x, n} = findzero( x );
Write( "\!ndone, a root of " || Char( Name Expr( f ) ) || " for V = " || complex:toChar( v, 10, 5 ) || " is ", complex:toChar( x, 10, 5 ) );
Eval( complex:convertExpression( Name Expr( f ) ) );第 3 行调用 convertExpression 函数将 +、-、*、/、^ 运算符从实数版本更改为复数版本。 V 设置为复数值，初始起点设置为远离实数线的点。 完成后，日志显示
完成后，对于 V = -100+0i，x ^ 2 - V 的根是 0+10i
[0 0]
V=-100 + 0i 是实数 -100（没有虚部）。0 + 10i 是没有实部的纯虚数答案；正如预期的那样，10i * 10i 是 100*(i^2) 或 -100。 复根以共轭对形式出现；牛顿的方法是根据起点找到一个。 另一个是-10i。
让我们修改代码以跟踪复平面中一堆点的迭代计数 n 并将其绘制为等值线图，并使用三次而不是二次：
f = Expr( x ^ 3 - V );
step = Eval Expr( x - Expr( Name Expr( f ) ) / Expr( Derivative( Name Expr( f ), x ) ) );
step = complex:convertExpression( Name Expr( step ) );
V = complex:make( -2, 0 ); // find the sqrt of V==-2
dt = New Table( "complex",New Column( "xre" ),New Column( "yim" ),New Column( "iterations" ),New Column( "basinr" ),New Column( "basini" ));
stepsize=1/32;
for(xreal=-2,xreal<=2,xreal+=stepsize,
for(yimag=-2,yimag<=2,yimag+=stepsize,
x = complex:make( xreal, yimag ); // initial guess. Does not like 0i for a start value.
{x, n} = findzero( x );
dt<<ADDROWS></ADDROWS>复平面上所有迭代到同一个根的点形成一个盆地。 上面，图形生成器将实值和虚值组合成一个值来识别根。
图像显示了红色、青色和深蓝色的三个盆地。 所有红色的起点都收敛到相同的值。 较大的点尺寸需要更长的时间才能收敛。 200 次迭代后灰色点没有收敛：
Three basins. Larger dots where the basins touch take longer to converge.三个盆。 盆地接触的较大点需要更长的时间才能收敛。
使用 JSL 复数实现将图像大小从 16K 点加倍到 65K 点需要 10 分钟左右。 我想制作 10,000 个 4Kx2K 点的位图，所以我需要另一种方法。
贝奥武夫集群
一个贝奥武夫集群是一堆旧电脑，大部分没有显示器或键盘，但通过局域网连接到控制计算机。 在我的例子中，控制计算机使用 JSL 来告诉集群要做什么并收集结果。 基于 Linux 的集群正在运行 JMP 创建的 Python 程序。 运行 JMP 的 Windows 计算机使用 RunProgram 中的 PLINK 在每台计算机上打开一个或多个 SSH 会话并提供所需的程序和参数。 Python 程序运行，然后将结果存储在 NFS 磁盘上，然后退出。 当 JMP 看到 RunProgram 终止时，它会收集结果并启动下一个 RunProgram。 10,000 个运行程序之后，就有足够的图像制作一个五分钟的视频。
Ethernet hub in center, below power-strip. Not enough mice and monitors and keyboards, SSH is best.以太网集线器在中心，电源板下方。 没有足够的鼠标、显示器和键盘，SSH 是最好的。
Python 有一个原生的复数数据类型，所以上面的一些 JSL 消失了。 JSL 现在看起来像这样；稍后您会看到一些 Python 嵌入为字符串。
F = Expr( // f and e parameters are driven by the frame animation time
f6 * xOld ^ e6 + f5 * xOld ^ e5 + f4 * xOld ^ e4 + f3 * xOld ^ e3 + f2 * xOld ^ e2 + f1 * xOld ^ e1 - 1
); // <<<<<<<<<<< supply equation here <<<<<<<<<<<<<<
scale = 2 / 1;
xresolution = 1920 * scale;
yresolution = 1080 * scale;
// keep aspect ratio by chopping shorter dimension...
xcenter = -1;
ycenter = 0;
delta = 3;
ytop = ycenter - delta * yresolution / xresolution;
xleft = xcenter - delta;
xright = xcenter + delta;
ybottom = ycenter + delta * yresolution / xresolution;
上面的位包括公式和一些视频大小参数。 下面的位描述了网络。 Hostkey、server、cpus、running、gigs 是关于集群信息的 condat（连接数据）列表的索引。
// named indexes for condat
hostkey = 1; // from a failed connection, this makes the connection work
server = 2; // the machine names
cpus = 3; // from the lscpu command, below
running = 4; // dynamic: how many are running on this cpu, right now
gigs = 5; // from the free command, below
// the avahi-daemon may need to be stop/start/status sometimes,
// unfortunately that is hard to do except at server keyboard if you can't connect...
守护进程的问题仍然存在，但我现在知道它正在为带有 -2（等）后缀的机器做广告，我可以在不连接显示器的情况下修复它。
condat = {
{"4e:c3:8a:3b:e7:d4:bd:0d:85:62:36:40:17:15:b4:a3", "dell790a.local", 4, 0, 22},
{"00:c2:37:09:dc:d7:cf:fb:33:e3:14:e3:d6:fb:39:be", "dell790b.local", 4, 0, 14},
{"54:c8:1a:6a:50:e3:51:02:ed:1b:f8:fa:bd:f4:e4:17", "dell7010a.local", 8, 0, 30},
{"66:0d:ca:3c:1f:11:eb:45:d4:92:d0:0b:ae:01:b8:6b", "dell7010b.local", 8, 0, 30},
{"94:84:2f:a5:56:2d:db:16:b5:80:a9:0d:87:07:ee:53", "dell9010.local", 8, 0, 30},
{"b2:b7:b3:c3:83:63:11:e7:56:e7:c9:c2:e6:fc:f0:c9", "thinkb.local", 8, 0, 30},
{"2b:ca:fc:45:f0:64:cd:ac:70:ce:93:73:e8:84:2d:75", "thinkc.local", 8, 0, 30},
{"62:b7:67:23:8d:86:fe:80:1a:62:ec:d8:41:1b:b4:1b", "thinkd.local", 8, 0, 30},
{"8e:54:ee:b0:25:67:b8:5c:4d:ad:5e:e9:71:87:fd:71", "netflix.local", 2, 0, 30} // host for vm
};下面的 JSL 查看 condat 列表，并将 dell9010 之类的名称转换为它们的点 IP 地址（或 IPV6 等效地址）。它可以防止以后查找失败；集群运行了大约三天，我仍然不确定为什么我妻子的桌面提供 DNS 查找服务……在我记得 GetAddrInfo 之前，我使用了注释的 ping 代码，这有点笨拙。 我的房子现在使用 IPV6 比几年前更满意了。 有时 IPV6 地址属于方括号。
// on my local network, the lookups for xxx.local fail randomly,
// and more often when busy. pre-lookup the dotted ip now.
//
// better way? Get Addr Info( "thinkd.local" )[3][4]
//
for( itry = 1, itry < 10, itry += 1,
nprocs = 0;
MAXCPU = 0;
good = 0;
for( iCondat = 1, iCondat <= N Items( condat ), iCondat += 1,
nprocs += condat[iCondat][cpus]; // nprocs will divide up the work (stepsize)
condat[iCondat][running] = 0;
MAXCPU = Max( MAXCPU, condat[iCondat][cpus] ); // MAXCPU loops needed later
if( Starts With( condat[iCondat][server], "[" ),
good += 1, // else
// txt = Run Program(
// Executable( "PING.EXE"/*path probably not needed*/ ),
// Options( {"-n 1", condat[iCondat][server]} ),
// ReadFunction( "text" )
// );
// ip = Regex( txt, "\[\[([^]]+)\]]\", "\0" ); // typ: [2603:6081:1a02:1200:8973:e7be:605c:f753]
try(
ip = "[" || Get Addr Info( condat[iCondat][server] )[3][4] || "]";
if( !Is Missing( ip ),
Write( "\!n", condat[iCondat][server] || " --> " || ip );
condat[iCondat][server] = ip;
good += 1;
, //
Write( "\!n", condat[iCondat][server] || " lookup failed ", Left( txt, 70 ) )
); //
, // catch
show(exception_msg,iCondat,condat[iCondat][server],Get Addr Info( condat[iCondat][server] ));
);
);
);
if( good == N Items( condat ),
Break(),
Wait( .1 )
);
);
Write( "\!n",good," of ", N Items( condat )," names resolved." );下一点是关于凭证的：ppk 是 JMP 机器的专用密钥。 该密钥的公共部分位于集群中的每台机器上。 Plink 需要该私有位来签署与远程机器的 SSH 连接。 if(0) 代码在那里，所以我可以在所有机器上进行实验。 在一天结束时，我可能会向他们发送所有关闭命令。 如果 condat 列表包含正在运行该 JSL 脚本的 VirtualBox 的 Linux 机器（上面的第 9 行），请小心关闭命令。
remoteuser = "c@";
localkey = " -i c:\Users\v1\Desktop\keys\private.ppk ";
// test code examples...
if( 0, // sequential test using command
//command = " \!"python3.6 -m pip install numpy";
//command = " \!"python3.6 -m pip --version";
//command = " \!"python3.6 -m pip install -U pip";
//command = " \!"python3.6 --version";
command = " \!"hostname -A";
command = " \!"top -b -n 2|head -n 15";
// check if the common nfs mount is set up
//command = " \!"ls /mnt/nfsc"; // fstab: closet.local:/mnt/nfs /mnt/nfsc/ nfs defaults 0 0
// cycle the power; -r means restart. for this to work, must be sudo and not need password
//command = " \!"sudo shutdown -r now"; // visudo: c ALL=(ALL) NOPASSWD: /sbin/poweroff, /sbin/reboot, /sbin/shutdown
//command = " \!"sudo shutdown now"; // visudo: c ALL=(ALL) NOPASSWD: /sbin/poweroff, /sbin/reboot, /sbin/shutdown
// info about number of cores
//command = " \!"lscpu | egrep '^(Model name|Thread|CPU\(s\):)'";
// info about memory
//command = " \!"free -h";
for( iCondat = 1, iCondat <= N Items( condat ), iCondat += 1,
result = runprogram(
executable( "plink" ),
options(
localkey ||
/**/ " -hostkey " || condat[iCondat][hostkey] || /* this line hides the host key */
" " || remoteuser || condat[iCondat][server] || command
),
write function( Function( {this}, this << Write( "n\r\n" ) ) ),
readfunction( "blob" )
);
result = Regex( Blob To Char( result ), "(\r\n)+", "\r\n", GLOBALREPLACE );
Write( "\!n\!n", iCondat, " ", condat[iCondat][server], " ", result );
Wait( 0 );
);
Write( "\!n\!ntest done" );
);下一点是关于同时运行多个 RunProgram 副本。 ImgNum 是每个关联数组的键，是 Python 程序应该处理的图像的编号。 ImgNumToRP 拥有一个正在运行的 RunProgram 对象。 该对象在终止时移动到 imgNumToRPdone。 imgNumToResult 保存 Python 程序的文本输出，imgNumToMachine 保存机器的 id，imgNumToLabel 保存图片的标签。 （标签稍后应用，更多内容如下。）
imgNumToRP = [=> ]; // holder for active RunProgram objects
imgNumToRPdone = [=> ]; // holder for RunProgram objects that need a moment to terminate
imgNumToRPresult = [=> ""]; // result strings from RunProgram objects
imgNumToMachine = [=> 0]; // which machine is this running on?
imgNumToLabel = [=>""]; // f(x) labels for graphs
nextImageNumberToMakePicture = 0;
// also: keep track of zeros from previous image
unstick 函数处理在 RunProgram 终止后最后一次没有调用 ReadFunction 的 RunProgram。
unstick = function({},{imagenumber},
for( imagenumber = imgNumToRP << first, !Is Empty( imagenumber ), imagenumber = imgNumToRPnext,
imgNumToRPnext = imgNumToRP << next( imagenumber );
if( (imgNumToRP[imagenumber]) << canread == 0 & (imgNumToRP[imagenumber]) << canwrite == 0 & (imgNumToRP[imagenumber]) << isreadeof == 1,
imgNumToRPresult[imagenumber] = "\!n" || char(imagenumber) || ": " || imgNumToRPresult[imagenumber] || "\!n"; // NoInterleave save this with count and basin
imgNumToRPdone[imagenumber] = imgNumToRP[imagenumber];
imgNumToRP << removeitem( imagenumber );
//if(imgNumToMachine[imagenumber]!=machine,throw("b machine not match"||char(machine)||","||char(imgNumToMachine[imagenumber])));
condat[imgNumToMachine[imagenumber]][running] = condat[imgNumToMachine[imagenumber]][running] - 1;
imgNumToMachine[imagenumber] = 0;
try(
write("\!nunsticking ",imagenumber);
process( imagenumber ); // must be after...
,
show(exception_msg);
);
Break(); // restart with <<FIRST></FIRST>下面的 findmachine 函数查找没有启动太多 SSH 会话的机器。 unstick() 函数会定期调用以查找 10% 的卡住情况。
findMachine = function( {},{iCondat,n}, // find a machine that is not filled up
n=0;
while( 1,
for( iCondat = 1, iCondat <= N Items( condat ), iCondat += 1,
if( condat[iCondat][running] < condat[iCondat][cpus],
Return( iCondat )
)
);
if(mod(n+=1,60)==0,
write("wait for machine...");
// look for stuck
unstick();
);
Wait( 1 ); // a full second is fine; the processes run for 10s of seconds
)
);当 JMP 检测到 RunProgram 已完成时，将调用 process 函数。 Python 程序在其输出的末尾写上“哨兵”，这段代码检查它作为粗略检查一切是否顺利运行。 Bindata 是 Python 程序在 NFS 磁盘上的农场目录中生成的二进制文件。 Windows 和 JMP 在 \\VBOXSVR 路径中看到 NFS 磁盘，但 Linux 和 Python 程序在 //mnt 路径中看到它。 此代码将从二进制数据中生成图片，但它被分离到后处理步骤中。 现在它只报告为公式发现的零的数量。 二进制格式在python代码中有说明； JMP 读取它并使用 blobtomatrix 加载它，这里使用了两次。 如果未注释和需要，最后重新启动可能会起作用；我忘记取消注释它并且不需要它。
// process output from one thread in the farm
process = function( {imagenumber},
{bindata, start, stop, balanceBytes, matdata, countArray,basinArray,zeros},
//show("process",imagenumber);
if( Contains( imgNumToRPresult[imagenumber], "sentinel" ),
// process a file. maybe not this one if it is out of order... FIXME...nextImageNumberToMakePicture
bindata = Load Text File( "\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated", blob() );
matdata = blobtomatrix(
bindata,
"uint",
2,
"little",
2
);
countArray = J( yresolution, xresolution, 0 );
basinArray = J( yresolution, xresolution, 0 );
countArray[0] = matdata[1::nrows(countArray)*ncols(countArray), 1];
basinArray[0] = matdata[1::nrows(countArray)*ncols(countArray), 2];
balanceBytes = length(bindata) - 4*nrows(countArray)*ncols(countArray);
zeros = blobtomatrix(blobpeek(bindata,length(bindata)-balanceBytes),"float",8,"little",4);
write("\!nimagenumber=",imagenumber, if(!ismatrix(zeros)," zeros ***"," zeros="||char(nrows(zeros))));
deletefile("\\VBOXSVR\nfsc\farm\testpython"||char(imagenumber)||".py"); // source no longer needed
sentinels += 1;//
, // else
Write( imgNumToRPresult[imagenumber] );
Write( "\!nrelaunch disabled" || char(imagenumber) );
//launch( imagenumber );
)
);启动函数启动另一个 RunProgram。 它使用 findmachine 等待可用的机器，并为 Python 命令提供自定义 .py 文件和自定义参数。 .py 文件在下面进一步创建。 下面第 5 行的 throw 是一个健全性检查； JSL 足够复杂，我觉得我需要编写一些额外的代码来观察逻辑错误。
// launch a thread in the farm
launch = function( {imagenumber},{command,machine},
machine = findMachine();
//show("launch",imagenumber,machine);
if(imgNumToMachine[imagenumber]!=0,throw("image number already on machine"||char(imagenumber)||","||char(machine)));
imgNumToMachine[imagenumber] = machine;
condat[machine][running] = condat[machine][running] + 1;
//command = "\[python3 /mnt/nfsc/farm/testpython]\"||char(imagenumber)||"\[.py]\";
command = "\[hostname -A;python3 /mnt/nfsc/farm/testpython]\"||char(imagenumber)||"\[.py]\";
imgNumToRP[imagenumber] = runprogram(
executable( "plink" ),
options(
localkey || /* -v */ " -hostkey " || condat[machine][hostkey] || " " || remoteuser || condat[machine][server] || " " || command || //
" -o " || Char( 0 ) || // offset
" -s " || Char( 1 ) || // step
" -x " || Char( xresolution ) || //
" -y " || Char( yresolution ) || //
" -t " || Char( ytop ) || //
" -l " || Char( xleft ) || //
" -r " || Char( xright ) || //
" -b " || Char( ybottom ) || //
" -f " || "/mnt/nfsc/farm/binary" || char(imagenumber) || ".generated" || //
""
),parameter 关键字使其参数可用于 read 函数。 这就是您可以并行拥有多个 RunProgram 并让读取函数知道它们是哪一个的方法。
Parameter( evallist({imagenumber,machine}) ),read 函数捕获输出并尝试（大约 90% 的成功率）检测 RunProgram 是否已完成。 因为当此代码运行时 RunProgram 仍在终止并仍在运行，所以句柄被移动到 imgNumToRPDone 而不仅仅是被删除。 这使它可以干净地完成。 下面的 Parm 接收上面的参数值。 下面还有另一个健全性检查（抛出）。
readfunction(
function( {this, parm}, {txt, imagenumber,machine},
{imagenumber,machine} = parm ;
//show("readfunc",imagenumber,machine);
txt = "";
//write("readfun("||char(imagenumber)||")... ");
While( this << canread, txt = txt || (this << read) );
if( Length( txt ) > 0,
//show(txt);
imgNumToRPresult[imagenumber] = imgNumToRPresult[imagenumber] || txt
);
if( this << isreadeof & imgNumToRP << Contains( imagenumber ), // this *almost* always works, about 1:100 fail to get the readeof in the readrunction
imgNumToRPresult[imagenumber] = "\!n" || char(imagenumber) || ": " || imgNumToRPresult[imagenumber] || "\!n"; // NoInterleave save this with count and basin
process( imagenumber );
imgNumToRPdone[imagenumber] = imgNumToRP[imagenumber];
imgNumToRP << removeitem( imagenumber );
if(imgNumToMachine[imagenumber]!=machine,throw("a machine not match"||char(machine)||","||char(imgNumToMachine[imagenumber])));
condat[machine][running] = condat[machine][running] - 1;
imgNumToMachine[imagenumber] = 0;
);
)
)
);
(imgNumToRP[imagenumber]) << writeeof;
);makePythonCode 函数生成包含 oneStep 方程的 .py 文件。 python 程序包含一个 init() 函数，该函数解析命令行以获取图片大小和输出文件名信息。 数据以两种方式传递：命令行参数和公式中的烘焙。 可能所有信息都已被烘焙，但如果没有理由（如将 oneStep 方程编译到代码中）将参数烘焙到代码中，那么您只需要所有 RunPrograms 的代码副本，如果它完全由命令行参数控制。 Python 代码是这样编写的，因此如果我希望多台计算机处理一个图像，它可以解决一部分问题。 一次完成一个完整的图像而不是稍后重新组合部分会更快。
makePythonCode = function({imagenumber,oneStep},{filename,python},
python =
"\[
import cmath, math, sys, getopt
import numpy as np
def main():
offset, step, ixdim, iydim, top, bottom, left, right, filename = init() # NoInterleave offset will be 0 and step will be 1
iymax = iydim - 1
ixmax = ixdim - 1
yrange = top - bottom
xrange = right - left
result = np.zeros( (int(math.ceil((iydim*ixdim-offset)/step)),2),dtype=np.int16) # count, basin # NoInterleave simpler, no change required though
iresult = 0
ix = offset # NoInterleave like iy, now begins at zero
iy = 0
while ix > ixmax: # normalize before loop starts
ix = ix - ixdim
iy = iy + 1
y = bottom + yrange * iy / iymax # after iy is normalized, before loop
roundedComplexToBasinId = {}
basinIdToBasinSize = {} # NoInterleave this might be useful, but needs to be appended to file (now a constant size for main chunk)
nextBasinId = 1 # reserve basin 1 for "did not converge"
while iy <= iymax:
x = left + xrange * ix / ixmax
xOld = complex(x,y) #xOld - (xOld ** 5 - 1) / (5 * xOld ** 4)
count = 0
try: ## catch zero divide, rarely
for irun in range(1,1000): # NoInterleave prev stopped at 100
###xNew = xOld - (xOld ** 5 - 1) / (5 * xOld ** 4)
]\" || // maintain the indent and newlines!下一行是固定的…OneStep 是有效 Python 语句的 JSL 字符串，它根据 xOld 值和牛顿方程设置 xNew，就像上面两行的 ### 注释一样。
oneStep || "\[
count += 1
delta = abs(xNew - xOld)
dist = abs(xOld)
if delta < 1e-15 * dist:
break
xOld = xNew
xOld = complex( round(xOld.real,10), round(xOld.imag,10) )
if irun >= 999:
basinId = 1 # arbitrary, but available always. did not converge.
else:
try: ## catch new basin discovered, rarely
basinId = roundedComplexToBasinId[xOld]
except KeyError as e:
nextBasinId += 1
roundedComplexToBasinId[xOld] = nextBasinId
basinId = nextBasinId
basinIdToBasinSize[basinId] = 0
### print(offset, count, xOld, basinId)
except ZeroDivisionError as e:
basinId = -1
### print(offset, e)
try:
result[iresult,0] = count
except:
print("oops",ix,iy)
if basinId >= 0:
try:
basinIdToBasinSize[basinId] += 1
except KeyError as e:
basinIdToBasinSize[basinId] = 1 # happens with basinId==1
result[iresult,1] = basinId
iresult += 1
ix = ix + step # advance
### offset = offset + step # needed for return value
while ix > ixmax: # re-normalize
ix = ix - ixdim
iy = iy + 1
y = bottom + yrange * iy / iymax这是编写整数和双精度二进制文件的python代码。 双精度数跟随整数，以便 JMP 使用 blobToMatrix 轻松获取所有内容。
# output results
fileobj = open(filename, mode='wb') # NoInterleave maybe prepend a descriptor
result.tofile(fileobj) # NoInterleave append after this constant size block
# NoInterleave append roundedComplexToBasinId and basinIdToBasinSize rather than printing them, they should be same dimension
# 4 8-byte doubles: real,imag,basin,size
temp = np.zeros( (len(roundedComplexToBasinId),4),dtype=np.float64)
itemp = 0
for iz in roundedComplexToBasinId.keys():
temp[itemp,0] = iz.real
temp[itemp,1] = iz.imag
temp[itemp,2] = roundedComplexToBasinId[iz]
temp[itemp,3] = basinIdToBasinSize[roundedComplexToBasinId[iz]]
itemp += 1
temp.tofile(fileobj)
fileobj.close()
print("sentinel"); # so JSL can test if this ran to the end...
def init():
offset = 0
step = 1
xdim = 3
ydim = 2
top = -1.0
bottom = 1.0
left = -1.0
right = 1.0
filename = "temp.binary"
options, remainder = getopt.getopt(sys.argv[1:], 's:x:y:o:t:b:l:r:f:')
for opt, arg in options:
if opt == '-o':
offset = int(arg) # offset is the starting point
elif opt == '-s':
step = int(arg) # step to the next point
elif opt == '-x':
xdim = int(arg) # xdim field int size
elif opt == '-y':
ydim = int(arg) # ydim field int size
elif opt == '-t':
top = float(arg) # top of field user coords
elif opt == '-b':
bottom = float(arg) # bot field user coords
elif opt == '-l':
left = float(arg) # left field user coords
elif opt == '-r':
right = float(arg) # right field user coords
elif opt == '-f':
filename = arg # write binary file
print(" xdim=",xdim,
" ydim=",ydim,
" offset=",offset,
" step=",step,
" top=",top,
" bottom=",bottom,
" left=",left,
" right=",right)
return offset, step, xdim, ydim, top, bottom, left, right, filename
if __name__ == "__main__":
main()
]\";Python 程序以唯一的编号保存到 NFS 目录（Windows 知道它）。 不使用返回的文件名； Linux 版本的文件名由启动函数创建并通过 Plink 传递。 （我们现在回到 JSL...）
filename = Save Text File( "\\VBOXSVR\nfsc\farm\testpython"||char(imagenumber)||".py", python ); // NoInterleave - unique file name for each picture because the equation changes
return(filename);
);
下面的循环生成了 10,000 帧。
for( imagenumber = 0, imagenumber < 10000, imagenumber += 1, //<<<<< frame loop <<<<<<
write("\!nloop",imagenumber);
// for X^5 - 1 there are 5 roots; 1 real and 4 complex
// the python expression for x(n+1) = x(n) - F / F' is
// xNew = xOld - (xOld ** 5 - 1) / (5 * xOld ** 4)公式，并将 OneStep 转换为 Python 语句——这看起来很复杂，因为它将实际值替换为参数 F1...F6 和 E1...E6。 隐藏在最后的是从 JSL 的 ^ 到 Python 的 ** 的变化，这是在文本中完成的，而不是表达式，因为 JSL 无法将 ** 表示为表达式。 F1...E6 作为表达式处理。 Sin( imagenumber / 200 ) 大约是从 -1 到 1 的 10 个完整的正弦波。 由于除数不同，正弦波的频率也不同。 因子从 -2 到 2 平滑移动，指数围绕 1...6 目标平滑移动 +/- 1。 因此，从一张图片到另一张图片的方程式非常相似，但变化有点不可预测。
FF = Derivative( Name Expr( F ), xOld );
// onestep will be inserted into python below
oneStep = substitute( // "xNew = xOld - (xOld ** 5 - 1) / (5 * xOld ** 4)"
char( // "xNew = xOld - (xOld ^ 5 - 1) / (5 * xOld ^ 4)"
substitute( // char() adds the parens around F and FF
Expr( xNew = xOld - F / FF ), // build the statement with symbols, then
Expr( F ), Name Expr( F ), // replace F and FF symbols with the content
Expr( FF ), Name Expr( FF ), // of the function and derivative
Expr( F6 ), round(2 * Sin( imagenumber / (200+5) ) / 2,3),
Expr( F5 ), round(2 * Sin( imagenumber / (200+7) ) / 2,3),
Expr( F4 ), round(2 * Sin( imagenumber / (200+11) ) / 2,3),
Expr( F3 ), round(2 * Sin( imagenumber / (200+13) ) / 2,3),
Expr( F2 ), round(2 * Sin( imagenumber / (200+17) ) / 2,3),
Expr( F1 ), round(2 * Sin( imagenumber / (200+19) ) / 2,3),
Expr( E6 ), round(6.0 + Sin( imagenumber / (200+23) ) / 2,3),
Expr( E5 ), round(5.0 + Sin( imagenumber / (200+29) ) / 2,3),
Expr( E4 ), round(4.0 + Sin( imagenumber / (200+31) ) / 2,3),
Expr( E3 ), round(3.0 + Sin( imagenumber / (200+37) ) / 2,3),
Expr( E2 ), round(2.0 + Sin( imagenumber / (200+41) ) / 2,3),
Expr( E1 ), round(1.0 + Sin( imagenumber / (200+47) ) / 2,3)
) // sub F,FF
), // char
"^", "**"
); // substitute ^ in JSL to ** in pythonimgNumToLabel 关联数组稍后可用于标记图像。 简化 xOld->X...
imgNumToLabel[imagenumber]=substitute(regex(oneStep,"xNew = xOld - \((.*?)\) / ","\1")," ** ","^"," * ","·","xOld","X");
write("\!n",imgNumToLabel[imagenumber]);
//"/mnt/nfsc/farm/binary" || char(imagenumber) || ".generated"在开发过程中使用现有文件的测试来加快程序必须重新启动时的速度。 它可能不太正常，它没有用于最后的运行（三天！）
if(fileexists("\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated"),
write("Exists ","\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated");
continue();
);
makePythonCode(imagenumber,oneStep);这看起来像是不合时宜的废弃代码。
n = 0; // eventually this becomes the offset // NoInterleave offset is always 0
sentinels = 0; // NoInterleave "sentinel" might be all it returns
//write("Launching... ");这就是 10,000 个 RunPrograms 开始的地方。 50 或 60 将使用上面的 condat 值同时运行。
launch( imagenumber );
write("Launched. ")
);总结...等待最后几个 RunPrograms 结束，unstick 可能会找到一两个。
while( N Items( imgNumToRP ) > 0,
Wait( .1 );
unstick();
);
write("\!ndone");稍后会为图像生成器保存标签。 上面的 JSL/Python 已经保存了仍然需要转换为 .png 格式的二进制数据。
savetextfile("$desktop/imgNumToLabel.txt",char(imgNumToLabel));Rpdone 持有 10,000 个已终止的 RunProgram。 放手吧。
imgNumToRPdone = 0;// clean out // NoInterleave picture generation happens after all count and basin are created
图像生成器可以现在运行，也可以稍后运行...
include("pythonFarmDisplay.jsl");Beowulf 代码到此结束。 图片生成器可以使用它，但没有。 我希望视频从一帧到下一帧使用相同的颜色。 颜色代表吸引力盆地，盆地中任何一点被吸引到的零点都应该是相同的颜色。 但是方程在每一帧上都会发生一点变化，所以零点在移动。 图片生成器尝试将前一张图片上的零点与下一张图片匹配。 有时他们会跳，但这已经足够了。
制作图片的 JSL 如下。 首先为盆地制作一个颜色表。 这些是颜色的最轻版本；由于迭代次数很大，它们稍后会变暗。
colors=
hlscolor(0,L,1)|| // dark red
hlscolor(.05,L,0.8)|| // red
hlscolor(.15,L,0.9)|| // orn
hlscolor(.35,L,1.0)|| // yel
hlscolor(.50,L,1.0)|| // grn
hlscolor(.65,L,1.0)|| // cyn
hlscolor(.80,L,0.7)|| // blu
hlscolor(.95,L,1.0)|| // mag
hlscolor(1.0,L,0.0); // gry
redmap=j(ncols(colors),1,0);
greenmap=j(ncols(colors),1,0);
bluemap=j(ncols(colors),1,0);
for(i=1,i<=ncols(colors),i+=1,
{r,g,b}=colortorgb(colors[i]);
m = max(r,g,b);
redmap[i] = r/m;
greenmap[i] = g/m;
bluemap[i] = b/m;
);之前保存的图片标签需要重新加载，图片尺寸需要匹配。 这段代码应该和前面的一样，并且应该已经被提取到一个通用的参数文件中。 现在太晚了！
imgNumToLabel=Include( "$desktop/pythonFarmDisplayLabels.jsl" );
scale = 2 / 1;
xresolution = 1920 * scale;
yresolution = 1080 * scale;
// keep aspect ratio by chopping shorter dimension...
xcenter = -1;
ycenter = 0;
delta = 3;
ytop = ycenter - delta * yresolution / xresolution;
xleft = xcenter - delta;
xright = xcenter + delta;
ybottom = ycenter + delta * yresolution / xresolution;
盆地的颜色应该在略有不同的图片/方程式中进行追踪。 这是上一张图片的零点所在的位置。
prevzeros = []; // what came before the first one
NFS 磁盘太小了。 这里有 309GB 的未压缩像素信息。
//filepat = "\\VBOXSVR\nfsc\farm\binary^imagenumber^.generated";
// manually moved the files (running out of space)...
filepat = "f:\newtonmethod\generated\binary^imagenumber^.generated";
countarray = J( yresolution, xresolution, 0 );
basinarray = J( yresolution, xresolution, 0 );
二进制文件将被解压缩成两个 JSL 矩阵：一个像素需要多少次迭代的计数，以及一个标识像素盆地的小整数。 之后还有一点二进制数据（平衡字节）将盆地编号带到一个复杂的坐标；复杂的坐标是从一张图片到下一张图片的匹配。
For( imagenumber = 1, imagenumber < 10000, imagenumber += 1,
file = Eval Insert( filepat );
If( File Exists( file ),
bindata = Load Text File( file, blob() );
matdata = Blob To Matrix( bindata, "uint", 2, "little", 2 );
countArray[0] = matdata[1 :: yresolution * xresolution, 1];
basinArray[0] = matdata[1 :: yresolution * xresolution, 2];
balanceBytes = Length( bindata ) - 4 * yresolution * xresolution;
zeros = Blob To Matrix( Blob Peek( bindata, Length( bindata ) - balanceBytes ), "float", 8, "little", 4 );
//Show( zeros );
if(zeros[1, 3]==2,
zeros = [9e99 9e99 1 0] |/ zeros; // python code should have put this in
);
// re-assign the basin numbers to line up with prevzeros
If( N Rows( prevzeros ) == 0,
prevzeros = zeros
);
KDTable 找到最近的邻居；复平面是一个很好的邻域。
kdt = KDTable( prevzeros[0, 1 :: 2] ); // 1::2 means the real and imag cols
newbasin = (1 :: N Rows( zeros )) + 10000;
toomany = N Rows( zeros ) > N Rows( prevzeros );
For( i = 1, i <= N Rows( zeros ), i += 1,
If( zeros[i, 3] != i,
Throw( "basin index assumption wrong" )
); // zeros are sorted by basin number
{indx, dist} = kdt << knearestrows( 2, zeros[i, 1 :: 2] );
If( N Cols( dist ) > 1 & dist[1] < dist[2] / 10, // assume 10X closer is a good match
newbasin[i] = indx[1];
kdt << removerows( indx[1] );//
, /*else*/ !toomany, // if there are toomany new zeros, don't assign them this way
newbasin[i] = indx[1];
kdt << removerows( indx[1] );
);
);<
/*
zerobasin newbasin
1 4 change 1s to 4s in basinArray and zeros
2 3 2s to 3s
3 10001 too many case! no basin assigned. becomes 7, below
4 2
5 5
6 6
7 1
*/
unassigned = Loc( newbasin > 10000 );
If( N Rows( unassigned ) > 0,新图片比旧图片有更多的零，把它们放在最后一个盆地。
newbasin[unassigned] = newbasin[unassigned] - 10000 + N Rows( zeros )
);
basinArray = Shape( newbasin[basinArray], N Rows( basinArray ), N Cols( basinArray ) );
zeros[0, 3] = Shape( newbasin[zeros[0, 3]], N Rows( zeros ), 1 );
ri = Rank Index( zeros[0, 3] );
zeros[0, 0] = zeros[ri, 0]; // zeros are sorted by basin number
零在上面排序；第 1 行获得颜色 1，依此类推。
prevzeros = zeros;
丢弃计数高的图片以使视频更流畅。
countarray[loc(countarray>200)] = 200; // 150 is typical big count, 999 is not converged
rangecount = Range( countarray );
countarray=countarray-3; // 3 or 4 is typical minimum count, make it bright (0)
countarray[loc(countarray<0)]=0;
GrayArray 是从计数中缩放的黑暗
grayArray = (1 - (countArray /*- rangeCount[1]*/) / (200 /*rangeCount[2] - rangeCount[1]*/)) ^ 8;
rangebasin = range(basinarray);
//basinArray = basinArray / (7);
nc=ncols(colors);
basinarray[loc(basinarray>nc)]=ncols(colors); // if these exist they could be mapped to gray, perhaps
basinarray[loc(basinarray<1)]=ncols(colors); // ditto
gray * redmap[basinarray] (etc) 是缩放的 rgb 颜色分量
red = grayArray :*</img> shape(redmap[basinarray], N Rows( basinArray ), N Cols( basinArray ) );//(1 + Sin( 2 * Pi() * basinArray + (0 * 2 * Pi() / 3) )) / 2;
green = grayArray :*</img> shape(greenmap[basinarray], N Rows( basinArray ), N Cols( basinArray ) );//(1 + Sin( 2 * Pi() * basinArray + (1 * 2 * Pi() / 3) )) / 2;
blue = grayArray :*</img> shape(bluemap[basinarray], N Rows( basinArray ), N Cols( basinArray ) );//(1 + Sin( 2 * Pi() * basinArray + (2 * 2 * Pi() / 3) )) / 2;
上溢/下溢 0..1 颜色值的修正…
red[Loc( red < 0 )] = 0;
red[Loc( red > 1 )] = 1;
green[Loc( green < 0 )] = 0;
green[Loc( green > 1 )] = 1;
blue[Loc( blue < 0 )] = 0;
blue[Loc( blue > 1 )] = 1;
+ 标记添加在零位置……粗略，但有效。 不要在边缘绘制，矩阵索引对零索引执行整行或整列！
crosssize = floor(yresolution / 100);
For( iz = 1, iz <= N Rows( zeros ), iz += 1,
real = Round( xresolution * Interpolate( Num( zeros[iz, 1] ), xleft, 0, xright, 1 ) );
imag = Round( yresolution * Interpolate( Num( zeros[iz, 2] ), ytop, 0, ybottom, 1 ) );
if(crosssize+5 < imag < nrows(red)-crosssize-5 & crosssize+5 < real < ncols(red)-crosssize-5,
red[imag - 1, real - crosssize :: real + crosssize] = 0.5;
red[imag + 1, real - crosssize :: real + crosssize] = 0.5;
red[imag - 2, real - crosssize :: real + crosssize] = 0.5;
red[imag + 2, real - crosssize :: real + crosssize] = 0.5;
red[imag - crosssize :: imag + crosssize, real - 1] = 0.5;
red[imag - crosssize :: imag + crosssize, real + 1] = 0.5;
red[imag - crosssize :: imag + crosssize, real - 2] = 0.5;
red[imag - crosssize :: imag + crosssize, real + 2] = 0.5;
green[imag - 1, real - crosssize :: real + crosssize] = 0.5;
green[imag + 1, real - crosssize :: real + crosssize] = 0.5;
green[imag - 2, real - crosssize :: real + crosssize] = 0.5;
green[imag + 2, real - crosssize :: real + crosssize] = 0.5;
green[imag - crosssize :: imag + crosssize, real - 1] = 0.5;
green[imag - crosssize :: imag + crosssize, real + 1] = 0.5;
green[imag - crosssize :: imag + crosssize, real - 2] = 0.5;
green[imag - crosssize :: imag + crosssize, real + 2] = 0.5;
blue[imag - 1, real - crosssize :: real + crosssize] = 0.5;
blue[imag + 1, real - crosssize :: real + crosssize] = 0.5;
blue[imag - 2, real - crosssize :: real + crosssize] = 0.5;
blue[imag + 2, real - crosssize :: real + crosssize] = 0.5;
blue[imag - crosssize :: imag + crosssize, real - 1] = 0.5;
blue[imag - crosssize :: imag + crosssize, real + 1] = 0.5;
blue[imag - crosssize :: imag + crosssize, real - 2] = 0.5;
blue[imag - crosssize :: imag + crosssize, real + 2] = 0.5;
red[imag, real - crosssize :: real + crosssize] = 1;
red[imag - crosssize :: imag + crosssize, real] = 1;
green[imag, real - crosssize :: real + crosssize] = 1;
green[imag - crosssize :: imag + crosssize, real] = 1;
blue[imag, real - crosssize :: real + crosssize] = 1;
blue[imag - crosssize :: imag + crosssize, real] = 1;
);
);
label = imgNumToLabel[imagenumber];标签添加在底部。 制作标签的小位图，然后将其复制到更大的图像中。
if(1,
//label = char(imagenumber)||" zeros="||char(N Rows( zeros ))||" "||label;
labelpic = textbox(label,setwrap(3000),<<SETFONTSIZE></SETFONTSIZE>
图像是在这里制作的。 然后通过命名决定是否将其保留在视频中。
img = New Image( "rgb", {red, green, blue} );
if(6<=nrows(zeros)<=8 & 2<=rangebasin[1]<=4 & 6<=rangebasin[2]<=8 & 3<=rangecount[1]<=4 & rangecount[2]<200,
img << saveimage( "f:/newtonmethod/png3/keeper" || Char( 1e6 + imagenumber ) || ".png", "png" );
,
img << saveimage( "f:/newtonmethod/png3/zapped" || Char( 1e6 + imagenumber ) || ".png", "png" );
);
//img<<SCALE></SCALE>而已！ 这是视频。 我决定只制作科学版本，在它们发生的地方有闪烁的黑框，下面有一个繁忙的标题。 一个漂亮的版本可能更流畅，甚至更催眠。 这是一个 4K 视频，如果你有硬件和连接速度，否则 youtube 会为你缩小它。
2020 年还剩两个小时。 我真的很期待2021年。
