345094_de-DE

liliconcom

Beowulf, Newton und Herr Hanson

Eines Tages im Jahr 1971 sagte Herr Hanson: „Heute werde ich Ihnen beibringen, wie man die Quadratwurzel einer Zahl zieht. Ein schnellerer Weg als die lange Teilung 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 den Schätzwert in die ursprüngliche Zahl, den Durchschnitt und die Wiederholung aufteilt. Hat er wohl erwähnt Newton-Raphson .

Newtons Methode zum Finden von Wurzeln 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

Xneu = 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 – anfängliche Vermutung

3.5 – Wende 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 hilft, die Quadratwurzel von V zu finden. Die zweite Zeile, Schritt, ist die Newton-Raphson-Iteration für einen Zyklus; step erscheint in der While-Schleife, wo er zum nächsten X-Wert ausgewertet wird. Die Ableitungsfunktion funktioniert hervorragend bei einfachen Polynomen und ergibt in diesem Beispiel 2X. Der NameExpr() erhält die Formel von f, ohne sie auszuwerten, und der Expr() lässt den NameExpr() laufen (um die Formel zu erhalten), wenn EvalExpr() den Expr() findet. Schritt ist auf x - f(x)/f'(x) gesetzt.

Innerhalb der While(1)-Schleife werden Schritte unternommen und gedruckt und gezählt. Der Test prüft die Differenz zwischen dem aktuellen und dem nächsten X-Wert; Wenn der Wert bei der aktuellen Größe von X sehr klein ist, kommt er wahrscheinlich nicht näher und kann zwischen zwei sehr nahen Zahlen oszillieren, sodass er aus der While-Schleife ausbricht.

0 -10.25

1 -10.0030487804878

2 -10.0000004646115

3-10

fertig, eine Wurzel aus 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 überschlug er sich auf -10,25 und pendelte sich dann auf -10 ein.

(Warum dividiert es durch die Ableitung f'? Wenn die Steigung der Linie f nahe Null ist, muss der nächste Punkt intuitiv weiter entfernt sein, um das Y-Ergebnis von f zu ändern, und die Division durch die Steigung tut dies.)

Die obige JSL verarbeitet unmögliche Anforderungen nicht gut. Wenn Sie aufgefordert werden, sqrt(-100) zu finden, macht die JSL eine Schleife und sucht nach einer komplexen Zahl, für deren Suche sie nicht ausgestattet ist.

 

Komplexe Zahlen

Im JSL-Scripting-Index 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 Matrix-Version ist näher an einer robusten Lösung und wahrscheinlich schneller. Die Klassenversion ist verlockend, aber ich hatte Probleme damit (Geschwindigkeit, Speicherprobleme). Da JMP keine Benutzerüberladungen für Operatoren hat, wird keine von beiden einfach zu verwenden sein.

Diese JSL ist der vorherigen JSL sehr ähnlich. Es gibt einen großen Codeblock, um die komplexe Zahlenarithmetik zu verarbeiten, dieselbe Formel, mehr oder weniger, und jetzt verarbeitet sie 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() und eine convertExpression()-Funktion. Großes Lob an StackExchange/Michael Hardy für eine 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 sind ziemlich aggressiv, wenn es darum geht, ((a+b)+(c+d)) in add(a,b,c,d) umzuwandeln, was keine binäre Funktion ist ... es hat vier Argumente. Aus diesem Grund gibt es zusätzliche Optionsargumente für die Funktionen complex:add und :sub.

Das nächste Bit von JSL verschiebt die Schleife in eine Funktion und antizipiert die zukünftige Entwicklung:

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 Anzahl der Iterationen n zurück.

Und der Code zum Aufrufen der Funktion … die Formel und die 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 Versionen mit reellen Zahlen in die Versionen mit komplexen Zahlen zu ändern. V wird auf einen komplexen Wert gesetzt und der anfängliche Startpunkt wird auf einen Punkt außerhalb der Geraden reeller Zahlen gesetzt. Wenn Sie fertig sind, wird das Protokoll angezeigt

fertig, eine Wurzel aus 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 treten in konjugierten Paaren auf; Newtons Methode fand einen basierend auf dem Startpunkt. Der andere ist bei -10i.

Lassen Sie uns den Code ändern, um die Iterationszahl n für eine Reihe von Punkten in der komplexen Ebene zu verfolgen und dies 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 zu derselben 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. Die größeren Punktgrößen brauchen länger zum Konvergieren. Die grauen Flecken konvergierten nach 200 Iterationen nicht:

Three basins. Larger dots where the basins touch take longer to converge.Drei Becken. Größere Punkte, an denen sich die Becken berühren, brauchen länger zum Konvergieren.

Das Verdoppeln der Bildgröße von 16.000 auf 65.000 Punkte dauert mit der JSL-Implementierung für komplexe Zahlen etwa 10 Minuten. Ich möchte 10.000 Bitmaps mit 4Kx2K-Punkten erstellen, also brauche ich einen anderen Ansatz.

 

Beowulf-Cluster

EIN Beowulf-Cluster ist ein Haufen alter Computer, meist ohne Bildschirme oder Tastaturen, aber über ein LAN mit einem steuernden Computer verbunden. In meinem Fall verwendet der steuernde Computer JSL, um dem Cluster mitzuteilen, was zu tun ist, 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 das erforderliche Programm und die erforderlichen Parameter bereitzustellen. Das Python-Programm wird ausgeführt, speichert die Ergebnisse auf einer NFS-Festplatte und wird dann beendet. Wenn JMP sieht, dass das RunProgram beendet wird, sammelt es die Ergebnisse und startet das nächste RunProgram. 10.000 RunPrograms später gibt es genug Bilder, um ein fünfminütiges Video zu machen.

Ethernet hub in center, below power-strip. Not enough mice and monitors and keyboards, SSH is best.Ethernet-Hub in der Mitte, unter der Steckdosenleiste. Nicht genug Mäuse und Monitore und Tastaturen, SSH ist am besten.

Python hat einen nativen Datentyp für komplexe Zahlen, sodass ein Teil der obigen JSL wegfällt. 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 Parameter für die Videogröße. Das folgende Bit beschreibt das Netzwerk. Hostkey, Server, CPUs, Running, 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 einem Suffix -2 (usw.) ankündigt, und ich kann es beheben, ohne ein Display 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 wandelt die Namen wie dell9010 in ihre gepunktete IP-Adresse (oder das IPV6-Äquivalent) um. Es verhindert spätere Suchfehler; der Cluster lief etwa drei Tage und ich bin mir immer noch nicht sicher, warum der Desktop meiner Frau den DNS-Lookup-Service 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. Irgendwann 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 Abschnitt geht es um Berechtigungsnachweise: ppk ist der private Schlüssel des JMP-Rechners. Der öffentliche Teil dieses Schlüssels befindet sich auf jeder Maschine im Cluster. Plink benötigt dieses private Bit, um die SSH-Verbindung zum Remote-Computer zu signieren. Der if(0)-Code ist vorhanden, damit ich auf allen Maschinen experimentieren kann. Am Ende des Tages könnte ich ihnen allen den Shutdown-Befehl schicken. Seien Sie vorsichtig mit dem Befehl shutdown, wenn die Condat-Liste den Linux-Computer enthält, auf dem 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 ist 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 Etiketten werden später angebracht, mehr 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 verarbeitet 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 folgende findmachine-Funktion 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 beendet wurde. Die Python-Programme schreiben „sentinel“ an das Ende ihrer Ausgabe und dieser Code überprüft es als groben Check, ob alles reibungslos läuft. Bindata ist die Binärdatei, die die Python-Programme im Farmverzeichnis auf der NFS-Festplatte erstellen. Windows und JMP sehen die NFS-Festplatte im \\VBOXSVR-Pfad, aber Linux und die Python-Programme sehen sie im //mnt-Pfad. Dieser Code würde die Bilder aus den Binärdaten erzeugen, aber das wurde in einen Nachbearbeitungsschritt getrennt. Jetzt meldet es nur die Anzahl der für die Formel entdeckten Nullen. Das Binärformat wird im Python-Code erklärt; JMP liest es und lädt es mit Blobtomatrix, die hier zweimal verwendet wird. Der Relaunch am Ende funktioniert wahrscheinlich, wenn er unkommentiert und benötigt wird; Ich habe vergessen, es zu entkommentieren 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 Launch-Funktion startet ein weiteres RunProgram. Es verwendet findmachine, um auf einen verfügbaren Computer zu warten, und stellt dem Python-Befehl eine benutzerdefinierte .py-Datei und benutzerdefinierte Parameter bereit. Die .py-Datei wird weiter unten erstellt. Der Wurf in Zeile 5 unten ist eine Plausibilitätsprüfung; Die JSL ist kompliziert genug, ich hatte das Gefühl, ich müsste 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 haben und den Lesefunktionen mitteilen, um welches es sich handelt.

  Parameter( evallist({imagenumber,machine}) ),

Die read-Funktion erfasst die Ausgabe und versucht (etwa 90 % der Zeit erfolgreich), zu erkennen, dass das RunProgram beendet ist. Da das RunProgram immer noch beendet wird und weiterhin ausgeführt wird, wenn dieser Code ausgeführt wird, wird das Handle nach imgNumToRPDone verschoben, anstatt nur entfernt zu werden. Dadurch lässt es sich sauber abschließen. Parm unten erhält den Parameterwert oben. Es gibt eine weitere Plausibilitätsprüfung (Wurf) unten.

  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 makePythonCode-Funktion generiert die .py-Datei mit der eingebrannten oneStep-Gleichung. Das Python-Programm enthält eine init()-Funktion, die die Befehlszeile analysiert, um Informationen zur Bildgröße und zum Ausgabedateinamen zu erhalten. Daten werden auf zwei Arten übergeben: Befehlszeilenparameter und in Gleichungen gebacken. Möglicherweise hätten alle Informationen eingebrannt werden können, aber wenn es keinen Grund gibt (wie das Kompilieren der OneStep-Gleichung in den Code), um Parameter in den Code einzubacken, benötigen Sie nur eine Kopie des Codes für alle RunPrograms, falls dies der Fall ist wurde vollständig von Befehlszeilenparametern gesteuert. Der Python-Code ist so geschrieben, dass er einen Teil des Problems lösen kann, wenn ich wollte, dass mehrere Computer an einem einzigen Bild arbeiten. Es ist insgesamt schneller, ein vollständiges Bild auf einmal zu erstellen und Teile später nicht 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 eingebackene Bit … OneStep ist eine JSL-Zeichenfolge einer gültigen Python-Anweisung, die xNew aus dem xOld-Wert und der Newton-Gleichung setzt, wie der ###-Kommentar zwei Zeilen weiter 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 / iymax

Hier ist der Python-Code, der eine Binärdatei mit Ganzzahlen und Doubles schreibt. Die Doubles folgen den Ganzzahlen, 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 im NFS-Verzeichnis (wie Windows es kennt) mit einer eindeutigen Nummer 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 das Konvertieren von OneStep in eine Python-Anweisung – das sieht kompliziert aus, weil es tatsächliche Werte in die Parameter F1...F6 und E1...E6 einsetzt. Am Ende ist der Wechsel von JSLs ^ zu Pythons ** verborgen, was im Text und nicht im Ausdruck erfolgt, da JSL keine Möglichkeit hat, ** als Ausdruck darzustellen. Die F1...E6 werden als Ausdrücke behandelt. Sin( imagenumber / 200 ) ist ungefähr 10 vollständige Sinuswelle von -1 bis 1. Die Sinuswellen haben aufgrund der unterschiedlichen Teiler unterschiedliche Frequenzen. Die Faktoren bewegen sich fließend von -2 bis 2 und die Exponenten bewegen sich fließend um +/- 1 um ihr Ziel 1...6. Die Gleichung von einem Bild zum nächsten ist also 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 python

Das assoziative Array imgNumToLabel kann später verwendet werden, um Bilder zu kennzeichnen. xAlt->X vereinfachen...

 imgNumToLabel[imagenumber]=substitute(regex(oneStep,"xNew = xOld - \((.*?)\) / ","\1")," ** ","^"," * ","·","xOld","X");
 write("\!n",imgNumToLabel[imagenumber]);
 
 //"/mnt/nfsc/farm/binary" || char(imagenumber) || ".generated"

Der Test für eine vorhandene Datei wurde während der Entwicklung verwendet, um die Dinge zu beschleunigen, wenn das Programm neu gestartet werden musste. Es funktioniert möglicherweise nicht ganz richtig, es wurde nicht für den letzten Lauf (drei Tage!)

 if(fileexists("\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated"),
  write("Exists ","\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated");
  continue();
 );
 
 makePythonCode(imagenumber,oneStep);

Dies sieht aus wie ein fehl am Platz aufgegebener Code.

 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 laufen gleichzeitig mit den condat-Werten ganz oben.

 launch( imagenumber );
 write("Launched. ")
);

Die Zusammenfassung ... warten Sie, bis die letzten paar RunPrograms beendet sind, lösen Sie möglicherweise ein oder zwei.

while( N Items( imgNumToRP ) > 0,
 Wait( .1 );
 unstick(); 
);

write("\!ndone");

Die Labels 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));

Rpdone hä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 verwenden können, tat es aber nicht. Ich möchte, dass das Video von einem Frame zum nächsten die gleichen Farben verwendet. Die Farben stellen Anziehungsbecken dar, die Null, von der jeder Punkt im Becken angezogen wird, sollte die gleiche Farbe erhalten. Aber die Gleichung ändert sich bei jedem Frame ein wenig, also bewegen sich die Nullen. Der Bildgenerator versucht, die Nullen des vorherigen Bildes mit dem nächsten Bild abzugleichen. Manchmal springen sie, aber es ist gut genug.

Die JSL, um die Bilder zu machen, folgt. Erstellen Sie zunächst eine Farbtabelle für die Becken. Dies ist die hellstmögliche Version der Farbe; sie werden später abgedunkelt, da 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 Bildetiketten müssen neu geladen werden, die Bildabmessungen müssen übereinstimmen. Dieser Code sollte derselbe wie zuvor sein und sollte 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 über leicht unterschiedliche Bilder/Gleichungen laufen. Hier leben die Nullen des vorherigen Bildes.

prevzeros = []; // what came before the first one

Die NFS-Festplatte war viel zu klein. Hier gibt es 309 GB unkomprimierte Pixelinformationen.

//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 Zählung, wie viele Iterationen das Pixel benötigt, und eine kleine Ganzzahl, die das Becken des Pixels identifiziert. Danach gibt es etwas mehr Binärdaten (die Ausgleichsbytes), die die Beckennummer in eine komplexe Koordinate bringen; Die komplexe Koordinate ist das, was von einem Bild zum nächsten passt.

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 feine 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] );
   );
  );

Der <

/*
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, legen Sie 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 bekommt 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 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 Fixups für Über-/Unterlauf der 0..1 Farbwerte…

  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 macht 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];

Das Etikett wird unten hinzugefügt. Erstellen Sie ein kleines 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 wird hier gemacht. Entscheiden Sie dann, ob es im Video 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 mich entschieden, nur die wissenschaftliche Version zu machen, mit flackernden dunklen Rahmen, wo sie passieren, und einer geschäftigen Bildunterschrift darunter. Eine hübsche Version könnte glatter und noch hypnotischer sein. Es ist ein 4K-Video, wenn Sie über die Hardware und Verbindungsgeschwindigkeit verfügen, andernfalls wird YouTube es für Sie herunterskalieren.

Es sind noch zwei Stunden von 2020 übrig, während ich dies poste. Ich freue mich sehr auf 2021.

Newtons Methode zum Auffinden von Nullstellen, komplexen Zahlen, Anziehungsbecken, Video