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年。
