Eines Tages im Jahr 1971 sagte Herr Hanson: „Heute werde ich Ihnen beibringen, wie man die Quadratwurzel einer Zahl berechnet.“ Ein schnellerer Weg als der Langteilung So wie ich unterrichten soll. Ich werde es dir mit einem Wort beibringen. Und du wirst es nie vergessen.“ Dann schreibt er an die Tafel
Erraten
und unterstreicht es. Als das Gelächter verstummte, erklärte er, wie man die Vermutung in die ursprüngliche Zahl, den Durchschnitt und die Wiederholung aufteilt. Er hat es wahrscheinlich erwähnt Newton-Raphson .
Newtons Methode zum Finden von Nullstellen oder Nullstellen einer Gleichung verwendet die Ableitung der Gleichung, um den Prozess zu beschleunigen. Die Quadratwurzel aus 10 sieht so aus, als würde man einen Wert von X finden, der X^2 == 10 oder X^2 – 10 == 0 ergibt. Die Ableitung von X^2 – 10 ist 2X und das Herzstück der Schleife zum Finden aufeinanderfolgender Werte von X ist
Xnew = X - ( (X^2 - 10) / (2X) )
was, wenn man damit spielt, genau das ist erraten Formel und konvergiert wie folgt für eine anfängliche Schätzung von 2:
2 – erste Vermutung
3.5 – wenden Sie x = (x + 10/x)/2 an
3.17857142857143
3.16231942215088
3.16227766044414
3.16227766016838
3.16227766016838 – sechs Iterationen zu einer stabilen Antwort!
Hier ist ein JSL-Programm zur Demonstration.
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 );
Die erste Zeile definiert einen Ausdruck f, der dabei hilft, die Quadratwurzel von V zu finden. Die zweite Zeile, Schritt, ist die Newton-Raphson-Iteration für einen Zyklus; Schritt erscheint in der while-Schleife, wo er zum nächsten X-Wert ausgewertet wird. Die Ableitungsfunktion eignet sich hervorragend für einfache Polynome und ergibt in diesem Beispiel 2X. NameExpr() ruft die Formel von f ab, ohne sie auszuwerten, und Expr() lässt NameExpr() ausführen (um die Formel zu erhalten), wenn EvalExpr() Expr() findet. Schritt ist auf x - f(x)/f'(x) eingestellt.
Innerhalb der while(1)-Schleife werden Schritte ausgeführt, gedruckt und gezählt. Der Test überprüft die Differenz zwischen dem aktuellen und dem nächsten X-Wert. Wenn der Wert bei der aktuellen Größe von
0 -10,25
1 -10.0030487804878
2 -10.0000004646115
3 -10
Fertig, eine Wurzel von x ^ 2 - V für V = 100 ist -10
Denken Sie daran, dass es sowohl eine positive als auch eine negative Quadratwurzel gibt! Startwerte bestimmen, zu welcher Wurzel die Konvergenz führt. Beginnend bei -8 stieg der Wert auf -10,25 und pendelte sich dann bei -10 ein.
(Warum wird durch die Ableitung f' dividiert? Wenn die Steigung der Linie f nahe bei Null liegt, muss der nächste Punkt intuitiv weiter entfernt sein, um das Y-Ergebnis von f zu ändern, und die Division durch die Steigung bewirkt dies.)
Die obige JSL verarbeitet unmögliche Anfragen nicht gut. Wenn die JSL aufgefordert wird, sqrt(-100) zu finden, durchläuft sie eine Schleife und sucht nach einer komplexen Zahl, für deren Suche sie nicht ausgestattet ist.
Komplexe Zahlen
Im JSL-Skriptindex führt die Suche nach „complex“ zu zwei verschiedenen unausgegorenen Beispielen. Eine basiert auf JSL-Klassen und eine auf Namespaces und Matrizen. Ich denke, die Matrixversion kommt einer robusteren Lösung näher und ist wahrscheinlich schneller. Die Klassenversion ist verlockend, aber ich hatte Probleme bei der Verwendung (Geschwindigkeit, Speicherprobleme). Da es in JMP keine Benutzerüberlastungen für Operatoren gibt, wird keines von beiden einfach zu verwenden sein.
Diese JSL ist der vorherigen JSL sehr ähnlich. Es gibt einen großen Codeblock für die komplexe Zahlenarithmetik, mehr oder weniger die gleiche Formel, und jetzt verarbeitet er auch 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() )
)
)
}
);Oben wurde eine pow()-Funktion für x^2 oder x^3 usw. hinzugefügt, wobei x eine komplexe Zahl ist, eine abs()-Funktion, die vielleicht Magnitude() heißen sollte, aber mit der nicht-komplexen Version als abs übereinstimmt, einige ungerade -suchende zusätzliche Parameter für add() und sub() sowie eine Funktion „convertExpression()“. Großes Lob an stackexchange/Michael Hardy für die leicht verständliche pow()-Funktion!
Die Funktion „convertExpression“ wandelt einen JSL-Ausdruck mit den Operatoren +, -, *, / und ^ in einen äquivalenten Ausdruck um, der die komplexen Namespace-Funktionen aufruft. Die JSL-Binäroperatoren gehen ziemlich aggressiv vor, wenn es darum geht, ((a+b)+(c+d)) in add(a,b,c,d) umzuwandeln, was keine Binärfunktion ist ... sie hat vier Argumente. Aus diesem Grund gibt es zusätzliche Optionsargumente für die Funktionen complex:add und :sub.
Der nächste Teil von JSL verschiebt die Schleife in eine Funktion und nimmt die zukünftige Entwicklung vorweg:
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 nimmt einen Startpunkt x und gibt eine Liste mit der iterierten Lösung x und der Iterationsanzahl n zurück.
Und der Code zum Aufrufen der Funktion … Formel und Ableitung sind für komplexe Zahlen gleich.
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 ) ) );Die dritte Zeile ruft die Funktion „convertExpression“ auf, um die Operatoren +, -, *, /, ^ von den reellen Zahlenversionen in die komplexen Zahlenversionen zu ändern. V wird auf einen komplexen Wert gesetzt und der anfängliche Startpunkt wird auf einen Punkt außerhalb der reellen Zahlenlinie gesetzt. Wenn Sie fertig sind, wird das Protokoll angezeigt
Fertig, eine Wurzel von x ^ 2 - V für V = -100+0i ist 0+10i
[0 0]
V=-100 + 0i ist die reelle Zahl -100 (ohne imaginäre Komponente). 0 + 10i ist die reine imaginäre Antwort ohne reelle Komponente; 10i * 10i ist 100*(i^2) oder -100, wie erwartet. Die komplexen Wurzeln kommen in konjugierten Paaren vor; Newtons Methode fand einen basierend auf dem Ausgangspunkt. Der andere liegt bei -10i.
Lassen Sie uns den Code ändern, um die Iterationsanzahl n für eine Reihe von Punkten in der komplexen Ebene zu verfolgen und diese grafisch darzustellen, vielleicht als Konturdiagramm, und kubisch statt quadratisch verwenden:
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>Alle Punkte in der komplexen Ebene, die zur gleichen Wurzel iterieren, bilden ein Becken. Oben kombiniert der Graph Builder die realen und imaginären Werte zu einem einzigen Wert, um eine Wurzel zu identifizieren.
Das Bild zeigt die drei Becken in Rot, Cyan und Dunkelblau. Alle roten Startpunkte konvergieren zum gleichen Wert. Bei größeren Punktgrößen dauert die Konvergenz länger. Die grauen Flecken konvergierten nach 200 Iterationen nicht:
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.Drei Becken. Bei größeren Punkten, an denen sich die Becken berühren, dauert es länger, bis sie zusammenlaufen.
Das Verdoppeln der Bildgröße von 16.000 auf 65.000 Punkte dauert mit der JSL-Implementierung komplexer Zahlen etwa 10 Minuten. Ich möchte 10.000 Bitmaps mit 4Kx2K-Punkten erstellen, daher brauche ich einen anderen Ansatz.
Beowulf-Cluster
A Beowulf-Cluster ist ein Haufen alter Computer, meist ohne Displays oder Tastaturen, aber über ein LAN mit einem steuernden Computer verbunden. In meinem Fall verwendet der steuernde Computer JSL, um dem Cluster mitzuteilen, was er tun soll, und um die Ergebnisse zu sammeln. Der Linux-basierte Cluster führt Python-Programme aus, die JMP erstellt. Der Windows-Computer, auf dem JMP ausgeführt wird, verwendet PLINK von RunProgram, um eine oder mehrere SSH-Sitzungen auf jedem Computer zu öffnen und die erforderlichen Programme und Parameter bereitzustellen. Das Python-Programm wird ausgeführt, speichert die Ergebnisse auf einer NFS-Festplatte und wird dann beendet. Wenn JMP feststellt, dass das RunProgram beendet wird, sammelt es die Ergebnisse und startet das nächste RunProgram. 10.000 RunPrograms später sind genug Bilder vorhanden, um ein fünfminütiges Video zu erstellen.
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.Ethernet-Hub in der Mitte, unterhalb der Steckdosenleiste. Nicht genügend Mäuse, Monitore und Tastaturen, SSH ist am besten.
Python verfügt über einen nativen Datentyp für komplexe Zahlen, sodass einige der oben genannten JSL wegfallen. Die JSL sieht jetzt so aus; Sie werden gleich etwas Python sehen, das als String eingebettet ist.
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;
Das obige Bit enthält die Gleichung und einige Videogrößenparameter. Das folgende Bit beschreibt das Netzwerk. Hostschlüssel, Server, CPU, Ausführung, Gigs sind Indizes in der Condat-Liste (Verbindungsdaten) mit Informationen über den Cluster.
// 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...
Das Problem mit dem Daemon besteht immer noch, aber ich weiß jetzt, dass er die Maschinen mit dem Suffix -2 (usw.) ankündigt, und ich kann es beheben, ohne einen Bildschirm anzuschließen.
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
};Die folgende JSL durchsucht die Condat-Liste und konvertiert Namen wie „dell9010“ in ihre gepunktete IP-Adresse (oder ein IPV6-Äquivalent). Dadurch werden spätere Suchfehler verhindert. Der Cluster lief etwa drei Tage lang und ich bin mir immer noch nicht sicher, warum der Desktop meiner Frau den DNS-Suchdienst bereitstellt … Bevor ich mich an GetAddrInfo erinnerte, habe ich den kommentierten Ping-Code verwendet, der etwas umständlicher ist. Mein Haus ist jetzt mit IPV6 viel zufriedener als noch vor ein paar Jahren. Manchmal gehört die IPV6-Adresse in eckige Klammern.
// 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." );Im nächsten Teil geht es um Anmeldeinformationen: Der PPK ist der private Schlüssel der JMP-Maschine. Der öffentliche Teil dieses Schlüssels befindet sich auf jedem Computer im Cluster. Plink benötigt dieses private Bit, um die SSH-Verbindung zum Remote-Computer zu signieren. Der if(0)-Code ist da, damit ich auf allen Maschinen experimentieren kann. Am Ende des Tages könnte ich ihnen allen den Befehl zum Herunterfahren senden. Seien Sie vorsichtig mit dem Befehl „shutdown“, wenn die Condat-Liste die Linux-Maschine enthält, auf der die VirtualBox ausgeführt wird, auf der dieses JSL-Skript ausgeführt wird (Zeile 9 oben).
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" );
);Im nächsten Teil geht es darum, viele Kopien von RunProgram gleichzeitig auszuführen. ImgNum ist der Schlüssel für jedes dieser assoziativen Arrays und die Nummer des Bildes, an dem das Python-Programm arbeiten soll. ImgNumToRP enthält ein RunProgram-Objekt, das ausgeführt wird. Dieses Objekt wird beim Beenden nach imgNumToRPdone verschoben. imgNumToResult enthält die Textausgabe des Python-Programms, imgNumToMachine die ID der Maschine und imgNumToLabel eine Bezeichnung für das Bild. (Die Beschriftungen werden später angebracht, mehr dazu weiter unten.)
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
Die Unstick-Funktion behandelt die RunPrograms, die die ReadFunction nicht ein letztes Mal aufrufen, nachdem RunProgram beendet wurde.
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>Die Funktion „findmachine“ unten sucht nach einer Maschine, auf der nicht zu viele SSH-Sitzungen gestartet sind. Die Funktion unstick() wird regelmäßig aufgerufen, um nach den 10 % der Fälle zu suchen, die hängen bleiben.
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
)
);Die Prozessfunktion wird aufgerufen, wenn JMP erkennt, dass ein RunProgram abgeschlossen wurde. Die Python-Programme schreiben „sentinel“ am Ende ihrer Ausgabe und dieser Code prüft dies als grobe Überprüfung, ob alles reibungslos läuft. Bindata ist die Binärdatei, die die Python-Programme im Farmverzeichnis auf der NFS-Festplatte erzeugen. Windows und JMP sehen die NFS-Festplatte im Pfad \\VBOXSVR, Linux und die Python-Programme sehen sie jedoch im Pfad //mnt. Dieser Code sollte die Bilder aus den Binärdaten erzeugen, dies wurde jedoch in einen Nachbearbeitungsschritt unterteilt. Jetzt wird nur noch die Anzahl der für die Formel entdeckten Nullen angezeigt. Das Binärformat wird im Python-Code erklärt; JMP liest es und lädt es mit Blobtomatrix, das hier zweimal verwendet wird. Der Relaunch am Ende funktioniert wahrscheinlich, wenn er unkommentiert und benötigt wird; Ich habe vergessen, es zu kommentieren und brauchte es nicht.
// 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 );
)
);Die Startfunktion startet ein weiteres RunProgram. Es verwendet findmachine, um auf eine verfügbare Maschine zu warten, und stellt dem Python-Befehl eine benutzerdefinierte .py-Datei und benutzerdefinierte Parameter zur Verfügung. Die .py-Datei wird weiter unten erstellt. Der Wurf in Zeile 5 unten ist eine Überprüfung der geistigen Gesundheit; Die JSL ist so kompliziert, dass ich das Gefühl hatte, ich müsste etwas zusätzlichen Code schreiben, um auf Logikfehler zu achten.
// 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" || //
""
),Das Parameter-Schlüsselwort stellt sein Argument der Lesefunktion zur Verfügung. Auf diese Weise können Sie mehrere RunPrograms parallel ausführen und den Lesefunktionen mitteilen, um welche es sich handelt.
Parameter( evallist({imagenumber,machine}) ),Die Lesefunktion erfasst die Ausgabe und versucht (in etwa 90 % der Fälle erfolgreich), zu erkennen, dass das RunProgram abgeschlossen ist. Da das RunProgram immer noch beendet wird und noch ausgeführt wird, wenn dieser Code ausgeführt wird, wird das Handle nach imgNumToRPDone verschoben und nicht einfach entfernt. Das ermöglicht einen sauberen Abschluss. Parm empfängt unten den Parameterwert oben. Unten finden Sie eine weitere Überprüfung der geistigen Gesundheit (Wurf).
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;
);Die Funktion makePythonCode generiert die .py-Datei mit der integrierten oneStep-Gleichung. Das Python-Programm enthält eine init()-Funktion, die die Befehlszeile analysiert, um Informationen zur Bildgröße und zum Namen der Ausgabedatei zu erhalten. Daten werden auf zwei Arten übergeben: über Befehlszeilenparameter und in Gleichungen eingebettet. Möglicherweise hätten alle Informationen integriert werden können, aber wenn es keinen Grund gibt (z. B. das Kompilieren der oneStep-Gleichung in den Code), Parameter in den Code zu integrieren, benötigen Sie nur eine Kopie des Codes für alle RunPrograms, wenn dies der Fall ist wurde vollständig durch Befehlszeilenparameter gesteuert. Der Python-Code ist so geschrieben, dass er einen Teil des Problems lösen kann, wenn mehrere Computer an einem einzelnen Bild arbeiten sollen. Insgesamt ist es schneller, ein vollständiges Bild auf einmal zu erstellen und Teile nicht später neu zu kombinieren.
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!Die nächste Zeile ist das eingebaute Bit … OneStep ist eine JSL-Zeichenfolge einer gültigen Python-Anweisung, die xNew aus dem xOld-Wert und der Newton-Gleichung festlegt, wie der ###-Kommentar zwei Zeilen oben.
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 / iymaxHier ist der Python-Code, der eine Binärdatei mit Ganzzahlen und Doppelzahlen schreibt. Die Doubles folgen den ganzen Zahlen, um es JMP zu erleichtern, alles mit blobToMatrix zu erfassen.
# 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()
]\";Das Python-Programm wird mit einer eindeutigen Nummer im NFS-Verzeichnis (wie Windows es kennt) gespeichert. Der zurückgegebene Dateiname wird nicht verwendet; Die Linux-Version des Dateinamens wird von der Startfunktion erstellt und über Plink übergeben. (Wir sind jetzt wieder in 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);
);
Die folgende Schleife generierte die 10.000 Frames.
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)Die Formel und die Konvertierung von OneStep in eine Python-Anweisung – das sieht kompliziert aus, weil es tatsächliche Werte in den Parametern F1...F6 und E1...E6 ersetzt. Am Ende verbirgt sich der Wechsel von JSLs ^ zu Pythons **, der im Text und nicht im Ausdruck erfolgt, da JSL keine Möglichkeit hat, ** als Ausdruck darzustellen. Die Tasten F1...E6 werden als Ausdrücke behandelt. Sin(Bildnummer/200) ist ungefähr 10 komplette Sinuswellen von -1 bis 1. Die Sinuswellen haben aufgrund der unterschiedlichen Teiler unterschiedliche Frequenzen. Die Faktoren bewegen sich sanft von -2 bis 2 und die Exponenten bewegen sich sanft um +/- 1 um ihr Ziel 1...6. Die Gleichung ist also von einem Bild zum nächsten ziemlich ähnlich, variiert aber etwas unvorhersehbar.
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 pythonDas assoziative Array imgNumToLabel kann später zum Beschriften von Bildern verwendet werden. Vereinfachen Sie xOld->X...
imgNumToLabel[imagenumber]=substitute(regex(oneStep,"xNew = xOld - \((.*?)\) / ","\1")," ** ","^"," * ","·","xOld","X");
write("\!n",imgNumToLabel[imagenumber]);
//"/mnt/nfsc/farm/binary" || char(imagenumber) || ".generated"Der Test einer vorhandenen Datei wurde während der Entwicklung verwendet, um die Arbeit zu beschleunigen, wenn das Programm neu gestartet werden musste. Möglicherweise funktioniert es nicht ganz richtig, es wurde beim letzten Lauf (drei Tage!) nicht verwendet.
if(fileexists("\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated"),
write("Exists ","\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated");
continue();
);
makePythonCode(imagenumber,oneStep);Das sieht nach fehl am Platz befindlichem, aufgegebenem Code aus.
n = 0; // eventually this becomes the offset // NoInterleave offset is always 0
sentinels = 0; // NoInterleave "sentinel" might be all it returns
//write("Launching... ");Hier beginnen die 10.000 RunPrograms. 50 oder 60 werden gleichzeitig mit den oben genannten Condat-Werten ausgeführt.
launch( imagenumber );
write("Launched. ")
);Die Zusammenfassung: Warten Sie, bis die letzten paar RunPrograms beendet sind. Wenn Sie den Stick entfernen, werden möglicherweise ein oder zwei gefunden.
while( N Items( imgNumToRP ) > 0,
Wait( .1 );
unstick();
);
write("\!ndone");Die Beschriftungen werden später für den Bildgenerator gespeichert. JSL/Python oben hat bereits Binärdaten gespeichert, die noch in das PNG-Format konvertiert werden müssen.
savetextfile("$desktop/imgNumToLabel.txt",char(imgNumToLabel));Der Rpdone enthält 10.000 beendete RunPrograms. Lassen Sie es gehen.
imgNumToRPdone = 0;// clean out // NoInterleave picture generation happens after all count and basin are created
Der Bildgenerator kann jetzt oder später ausgeführt werden ...
include("pythonFarmDisplay.jsl");Das ist das Ende des Beowulf-Codes. Der Bildgenerator hätte es gebrauchen können, tat es aber nicht. Ich möchte, dass das Video von Bild zu Bild die gleichen Farben verwendet. Die Farben stellen Anziehungsbecken dar. Der Nullpunkt, von dem jeder Punkt im Becken angezogen wird, sollte die gleiche Farbe erhalten. Aber die Gleichung ändert sich bei jedem Bild ein wenig, sodass sich die Nullen verschieben. Der Bildgenerator versucht, die Nullen des vorherigen Bildes mit dem nächsten Bild abzugleichen. Manchmal springen sie, aber es ist gut genug.
Es folgt die JSL zum Erstellen der Bilder. Erstellen Sie zunächst eine Farbtabelle für die Becken. Dies ist die hellste mögliche Version der Farbe; Sie werden später dunkler, wenn die Anzahl der Iterationen groß ist.
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;
);Die zuvor gespeicherten Bildbeschriftungen müssen neu geladen werden, die Bildabmessungen müssen übereinstimmen. Dieser Code sollte derselbe sein wie zuvor und in eine gemeinsame Parameterdatei extrahiert worden sein. Jetzt ist es zu spät!
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;
Die Farben für die Becken sollen sich auf leicht unterschiedliche Bilder/Gleichungen übertragen lassen. Hier leben die Nullen des vorherigen Bildes.
prevzeros = []; // what came before the first one
Die NFS-Festplatte war viel zu klein. Hier stehen 309 GB unkomprimierte Pixelinformationen zur Verfügung.
//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 );
Die Binärdateien werden in zwei JSL-Matrizen entpackt: eine Anzahl, wie viele Iterationen das Pixel benötigt, und eine kleine Ganzzahl, die das Becken des Pixels identifiziert. Danach gibt es noch ein paar Binärdaten (die Bilanzbytes), die die Beckennummer in eine komplexe Koordinate umwandeln; Die komplexe Koordinate ist das, was von einem Bild zum nächsten abgeglichen wird.
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
);
Die KDTable findet die nächsten Nachbarn; Die komplexe Ebene ist eine schöne Nachbarschaft.
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] );
);
);Das <
/*
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,Das neue Bild hat mehr Nullen als das alte Bild, lege sie in das letzte Becken.
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
Die Nullen sind oben sortiert; Reihe 1 erhält Farbe 1 usw.
prevzeros = zeros;
Die Bilder mit hoher Anzahl werden verworfen, um das Video flüssiger zu machen.
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;
Das GrayArray ist die skalierte Dunkelheit aus den Zählungen
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
Die graue * Redmap[basinarray] (usw.) ist eine skalierte RGB-Farbkomponente
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;
Die Korrekturen für Über-/Unterlauf der Farbwerte 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;
Die +-Markierungen werden an den Nullstellen hinzugefügt … grob, aber funktioniert. Zeichnen Sie nicht am Rand, die Matrixindizierung erstellt ganze Zeilen oder Spalten für den Index Null!
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];Die Beschriftung wird unten hinzugefügt. Erstellen Sie eine kleine Bitmap des Etiketts und kopieren Sie es dann in das größere Bild.
if(1,
//label = char(imagenumber)||" zeros="||char(N Rows( zeros ))||" "||label;
labelpic = textbox(label,setwrap(3000),<<SETFONTSIZE></SETFONTSIZE>
Das Bild ist hier entstanden. Entscheiden Sie dann, ob es im Video erhalten bleiben soll, indem Sie es benennen.
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>Das ist es! Hier ist das Video. Ich habe beschlossen, nur die wissenschaftliche Version zu erstellen, mit flackernden dunklen Rahmen dort, wo sie passieren, und einer geschäftigen Bildunterschrift darunter. Eine hübsche Version könnte sanfter und noch hypnotischer sein. Es handelt sich um ein 4K-Video, wenn Sie über die entsprechende Hardware und Verbindungsgeschwindigkeit verfügen, andernfalls skaliert YouTube es für Sie herunter.
Während ich dies poste, verbleiben noch zwei Stunden des Jahres 2020. Ich freue mich sehr auf 2021.
