1971 年的一天，汉森先生说：“今天我要教你如何求一个数的平方根。 比以下方法更快的方法长除法我应该以我应该教的方式。 我就用一句话来教你。 而且你永远不会忘记。”然后他在黑板上写下
猜测
并强调它。 当笑声平息后，他解释了如何将猜测除以原始数字，平均，重复。 他可能提到牛顿-拉夫逊。
牛顿求方程根或零点的方法使用方程的导数来加速该过程。 10 的平方根看起来就像找到一个使 X^2 == 10 或 X^2 – 10 == 0 的 X 值。 X^2 – 10 的导数是 2X，查找 X 的连续值的循环核心是
X新 = 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 迭代；步骤出现在 while 循环中，其计算结果为下一个 X 值。 微分函数在简单多项式上效果很好，在本例中给出 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 的斜率接近于 0 时，下一个点需要更远才能改变 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() 函数可能应该被称为模量（），但与非复数版本匹配，因为有一些奇怪的-寻找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 设置为复数值，初始起点设置为实数轴之外的点。 完成后，日志显示
完成后，x ^ 2 - V 的根（对于 V = -100+0i）是 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.Three basins. Larger dots where the basins touch take longer to converge.Three basins. Larger dots where the basins touch take longer to converge.Three basins. Larger dots where the basins touch take longer to converge.三个盆。 盆地接触的点越大，需要更长的时间才能汇聚。
使用 JSL 复数实现将图像大小从 16K 点加倍到 65K 点大约需要 10 分钟。 我想制作 10,000 个 4Kx2K 点的位图，所以我需要另一种方法。
贝奥武甫星团
A贝奥武甫星团是一堆旧计算机，大多没有显示器或键盘，但通过 LAN 连接到控制计算机。 在我的例子中，控制计算机使用 JSL 告诉集群要做什么并收集结果。 基于 Linux 的集群正在运行 JMP 创建的 Python 程序。 运行 JMP 的 Windows 计算机使用 RunProgram 中的 PLINK 在每台计算机上打开一个或多个 SSH 会话并提供所需的程序和参数。 Python 程序运行，然后将结果存储在 NFS 磁盘上，然后退出。 当 JMP 发现 RunProgram 终止时，它会收集结果并启动下一个 RunProgram。 10,000 个 RunPrograms 后就有足够的图像来制作一个五分钟的视频。
Ethernet hub in center, below power-strip. Not enough mice and monitors and keyboards, SSH is best.Ethernet hub in center, below power-strip. Not enough mice and monitors and keyboards, SSH is best.Ethernet hub in center, below power-strip. Not enough mice and monitors and keyboards, SSH is best.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 列表包含正在运行 VirtualBox 的 Linux 计算机，而 VirtualBox 正在运行此 JSL 脚本（上面第 9 行），请小心使用 shutdown 命令。
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 已完成时，将调用处理函数。 Python 程序在其输出的末尾写入“sentinel”，此代码会对其进行检查，以粗略地检查一切是否运行顺利。 Bindata 是 Python 程序在 NFS 磁盘上的 farm 目录中生成的二进制文件。 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 );
)
);launch 函数启动另一个 RunProgram。 它使用 findmachine 等待可用机器，并向 Python 命令提供自定义 .py 文件和自定义参数。 下面进一步创建 .py 文件。 下面第 5 行的抛出是健全性检查； 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" || //
""
),参数关键字使其参数可用于 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 方程编译到代码中）将参数烘焙到代码中，则您只需要所有 RunProgram 的一份代码副本（如果它是）完全由命令行参数控制。 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");贝奥武夫密码到此结束。 图片生成器本来可以使用它，但没有。 我希望视频从一帧到下一帧使用相同的颜色。 颜色代表吸引力盆地，盆地中任何点被吸引到的零都应该具有相同的颜色。 但方程在每一帧上都会发生一点变化，因此零点会四处移动。 图片生成器尝试将上一张图片与下一张图片上的零进行匹配。 有时他们会跳，但这已经足够了。
制作图片的 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
grey * 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年。
