Un giorno del 1971, il signor Hanson disse: “Oggi ti insegnerò come ottenere la radice quadrata di un numero. Un modo più veloce rispetto al divisione lunga modo in cui dovrei insegnare. Ti insegnerò con una parola. E non dimenticherai mai. Poi scrive alla lavagna
Indovinare
e lo sottolinea. Quando le risate si placarono, spiegò come dividere l'ipotesi nel numero originale, media, ripetizione. Probabilmente ha accennato Newton-Raphson .
Il metodo di Newton per trovare radici o zeri di un'equazione utilizza la derivata dell'equazione per accelerare il processo. La radice quadrata di 10 sembra trovare un valore di X che rende X^2 == 10 o X^2 – 10 == 0. La derivata di X^2 – 10 è 2X, e il cuore del ciclo per trovare valori successivi di X è
Xnuovo = X - ( (X^2 - 10) / (2X) )
che, se ci giochi, è esattamente il Indovinare formula, e converge in questo modo, per un'ipotesi iniziale di 2:
2 – ipotesi iniziale
3.5 – applica x = (x + 10/x)/2
3.17857142857143
3.16231942215088
3.16227766044414
3.16227766016838
3.16227766016838 – sei iterazioni per una risposta stabile!
Ecco un programma JSL da dimostrare.
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 );
La prima riga definisce un'espressione, f, che aiuterà a trovare la radice quadrata di V. La seconda riga, step, è l'iterazione di Newton-Raphson per un ciclo; step appare nel ciclo while dove restituisce il valore X successivo. La funzione Derivata funziona alla grande su polinomi semplici, dando 2X in questo esempio. NameExpr() ottiene la formula da f senza valutarla e Expr() esegue NameExpr() (per ottenere la formula) quando EvalExpr() trova Expr(). Il passo è impostato su x - f(x)/f'(x).
All'interno del ciclo while(1), i passi vengono eseguiti, stampati e contati. Il test controlla la differenza dei valori X correnti e successivi; se il valore è molto piccolo alla dimensione attuale di X, probabilmente non si avvicinerà e potrebbe oscillare tra due numeri molto vicini, quindi esce dal ciclo while.
0 -10.25
1-10.0030487804878
2-10.0000004646115
3-10
fatto, una radice di x ^ 2 - V per V = 100 è -10
Ricorda che esiste sia una radice quadrata positiva che una negativa! I valori iniziali determinano a quale radice porta la convergenza. A partire da -8 è andato oltre a -10,25, poi si è stabilizzato a -10.
(Perché divide per la derivata, f' ? Intuitivamente, quando la pendenza della retta f è prossima allo zero, il punto successivo dovrà essere più lontano per cambiare il risultato Y di f, e dividendo per la pendenza lo fa.)
Il JSL sopra non gestisce bene le richieste impossibili. Se viene richiesto di trovare sqrt(-100), il JSL esegue un ciclo, cercando un numero complesso che non è in grado di trovare.
Numeri complessi
Nell'indice di scripting JSL, la ricerca di complex produce due diversi esempi semicotti. Uno è basato su classi JSL e uno basato su spazi dei nomi e matrici. Penso che la versione a matrice sia più vicina a una soluzione robusta e probabilmente più veloce. La versione di classe è allettante, ma ho avuto problemi a usarla (velocità, problemi di memoria). Poiché JMP non ha sovraccarichi utente per gli operatori, nessuno dei due sarà facile da usare.
Questo JSL è molto simile al precedente JSL. C'è un grosso pezzo di codice per gestire l'aritmetica dei numeri complessi, la stessa formula, più o meno, e ora gestisce anche 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() )
)
)
}
);Sopra è stata aggiunta una funzione pow() per x^2 o x^3 ecc. dove x è un numero complesso, una funzione abs() che forse dovrebbe essere chiamata grandezza() ma corrisponde alla versione non complessa come abs, un po' strano -guardando parametri extra per add() e sub() e una funzione convertExpression(). Un grande ringraziamento a Stackexchange/Michael Hardy per una funzione pow() facile da comprendere!
La funzione convertExpression massaggia un'espressione JSL con gli operatori +, -, *, / e ^ in un'espressione equivalente che chiama le funzioni complesse dello spazio dei nomi. Gli operatori binari JSL sono piuttosto aggressivi nel trasformare ((a+b)+(c+d)) in add(a,b,c,d) che non è una funzione binaria... ha quattro argomenti. Questo è il motivo per cui ci sono argomenti di opzione extra sulle funzioni complex:add e :sub.
Il prossimo bit di JSL sposta il ciclo in una funzione, anticipando lo sviluppo futuro:
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 prende un punto di partenza x e restituisce una lista con la soluzione iterata x e il conteggio delle iterazioni n.
E il codice per chiamare la funzione...la formula e la derivata sono le stesse per i numeri complessi.
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 ) ) );La terza riga chiama la funzione convertExpression per modificare gli operatori +, -, *, /, ^ dalle versioni dei numeri reali alle versioni dei numeri complessi. V è impostato su un valore complesso e il punto di partenza iniziale è impostato su un punto fuori dalla retta dei numeri reali. Al termine, viene visualizzato il registro
fatto, una radice di x ^ 2 - V per V = -100+0i è 0+10i
[0 0]
V=-100 + 0i è il numero reale -100 (senza componente immaginaria). 0 + 10i è la pura risposta immaginaria senza componente reale; 10i * 10i è 100*(i^2) o -100, come previsto. Le radici complesse si verificano in coppie coniugate; Il metodo di Newton ne trovò uno basato sul punto di partenza. L'altro è a -10i.
Modifichiamo il codice per tenere traccia del numero di iterazioni n per un gruppo di punti nel piano complesso e nel grafico che, forse come grafico di contorno, e utilizziamo cubico invece di quadratico:
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>Tutti i punti nel piano complesso che iterano alla stessa radice formano un bacino. Sopra, il generatore di grafici combina i valori reali e immaginari in un unico valore per identificare una radice.
L'immagine mostra i tre bacini in rosso, ciano e blu scuro. Tutti i punti di partenza rossi convergono allo stesso valore. Le dimensioni dei punti più grandi impiegano più tempo a convergere. I punti grigi non convergevano dopo 200 iterazioni:
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.Tre vasche. I punti più grandi dove i bacini si toccano impiegano più tempo a convergere.
Raddoppiare la dimensione dell'immagine da 16K a 65K punti richiede circa 10 minuti con l'implementazione del numero complesso JSL. Voglio creare 10.000 bitmap di punti 4Kx2K, quindi ho bisogno di un altro approccio.
Ammasso Beowulf
UN Ammasso Beowulf è un mucchio di vecchi computer, per lo più senza display o tastiere, ma collegati su una LAN a un computer di controllo. Nel mio caso, il computer di controllo utilizza JSL per dire al cluster cosa fare e per raccogliere i risultati. Il cluster basato su Linux esegue i programmi Python creati da JMP. Il computer Windows che esegue JMP utilizza PLINK da RunProgram per aprire una o più sessioni SSH su ciascun computer e fornire il programma e i parametri necessari. Il programma Python viene eseguito, quindi memorizza i risultati su un disco NFS, quindi esce. Quando JMP vede il RunProgram terminare, raccoglie i risultati e avvia il successivo RunProgram. 10.000 RunProgram dopo ci sono abbastanza immagini per realizzare un video di cinque minuti.
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.Hub Ethernet al centro, sotto la presa multipla. Non abbastanza mouse, monitor e tastiere, SSH è il migliore.
Python ha un tipo di dati numerico complesso nativo, quindi parte del JSL sopra scompare. Il JSL ora ha questo aspetto; vedrai un po 'di Python incorporato come una stringa tra un po'.
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;
Il bit sopra include l'equazione e alcuni parametri di dimensione del video. Il bit seguente descrive la rete. Hostkey, server, cpus, running, gigs sono indici nell'elenco condat (dati di connessione) di informazioni sul 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...
Il problema con il demone esiste ancora, ma ora so che pubblicizza le macchine con un suffisso -2 (ecc.) e posso risolverlo senza collegare un display.
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
};Il JSL di seguito esamina l'elenco dei condat e converte i nomi come dell9010 nel loro indirizzo IP puntato (o equivalente IPV6). Previene gli errori di ricerca successivi; il cluster ha funzionato per circa tre giorni e non sono ancora sicuro del motivo per cui il desktop di mia moglie fornisce il servizio di ricerca DNS... Prima di ricordare GetAddrInfo ho usato il codice ping commentato, che è un po' più goffo. La mia casa è molto più felice con IPV6 ora rispetto a qualche anno fa. A volte l'indirizzo IPV6 appartiene tra parentesi quadre.
// 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." );La parte successiva riguarda le credenziali: ppk è la chiave privata della macchina JMP. La parte pubblica di quella chiave si trova su ogni macchina del cluster. Plink ha bisogno di quel bit privato per firmare la connessione SSH alla macchina remota. Il codice if(0) è lì così posso sperimentare su tutte le macchine. Alla fine della giornata, potrei inviare a tutti loro il comando di spegnimento. Fai attenzione con il comando shutdown se l'elenco condat contiene la macchina Linux che esegue VirtualBox che esegue questo script JSL (riga 9 sopra).
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" );
);La parte successiva riguarda l'esecuzione simultanea di molte copie di RunProgram. ImgNum è la chiave per ciascuno di questi array associativi ed è il numero dell'immagine su cui dovrebbe lavorare il programma Python. ImgNumToRP contiene un oggetto RunProgram in esecuzione. L'oggetto viene spostato in imgNumToRPdone mentre sta terminando. imgNumToResult contiene l'output di testo del programma Python, imgNumToMachine l'id della macchina e imgNumToLabel un'etichetta per l'immagine. (Le etichette vengono applicate in seguito, più sotto.)
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
La funzione unstick gestisce i RunPrograms che non chiamano ReadFunction un'ultima volta dopo che RunProgram termina.
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>La funzione findmachine qui sotto cerca una macchina che non ha troppe sessioni SSH avviate. La funzione unstick() viene chiamata periodicamente per cercare il 10% dei casi che si bloccano.
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
)
);La funzione di processo viene richiamata quando JMP rileva che un RunProgram è terminato. I programmi Python scrivono "sentinel" alla fine del loro output e questo codice lo controlla come un controllo grezzo che tutto funzioni senza intoppi. Bindata è il file binario che i programmi Python producono nella directory farm sul disco NFS. Windows e JMP vedono il disco NFS nel percorso \\VBOXSVR ma Linux ei programmi Python lo vedono nel percorso //mnt. Questo codice avrebbe prodotto le immagini dai dati binari, ma è stato separato in una fase di post-elaborazione. Ora riporta solo il numero di zeri scoperti per la formula. Il formato binario è spiegato nel codice Python; JMP lo legge e lo carica con blobtomatrix, qui utilizzato due volte. Il rilancio alla fine probabilmente funziona se non commentato e necessario; Ho dimenticato di decommentarlo e non ne avevo bisogno.
// 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 );
)
);La funzione di avvio avvia un altro RunProgram. Utilizza findmachine per attendere una macchina disponibile e fornisce un file .py personalizzato e parametri personalizzati al comando Python. Il file .py viene creato più avanti. Il lancio sulla riga 5 di seguito è un controllo di sanità mentale; il JSL è abbastanza complicato, mi sentivo come se avessi bisogno di scrivere del codice extra per controllare gli errori logici.
// 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" || //
""
),La parola chiave parameter rende disponibile il suo argomento alla funzione read. In questo modo è possibile avere più RunProgram in parallelo e far sapere alle funzioni di lettura quale sono.
Parameter( evallist({imagenumber,machine}) ),La funzione read acquisisce l'output e tenta (circa il 90% delle volte con successo) di rilevare che RunProgram è terminato. Poiché RunProgram è ancora in fase di chiusura ed è ancora in esecuzione quando viene eseguito questo codice, l'handle viene spostato in imgNumToRPDone anziché essere semplicemente rimosso. Ciò gli consente di finire in modo pulito. Parm, sotto, riceve il valore del parametro, sopra. C'è un altro controllo di sanità mentale (lancio) di seguito.
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;
);La funzione makePythonCode genera il file .py con l'equazione oneStep incorporata. Il programma python include una funzione init() che analizza la riga di comando per ottenere informazioni sulla dimensione dell'immagine e sul nome del file di output. I dati vengono passati in due modi: parametri della riga di comando e cotti in equazioni. Forse tutte le informazioni avrebbero potuto essere inserite, ma se non c'è motivo (come compilare l'equazione oneStep nel codice) per inserire i parametri nel codice, avresti bisogno solo di una copia del codice per tutti i RunProgram se era completamente controllato da parms della riga di comando. Il codice Python è scritto in modo che possa fare una parte del problema se volessi che più computer lavorassero su una singola immagine. Nel complesso è più veloce creare un'immagine completa in una sola volta e non ricombinare le parti in un secondo momento.
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!La riga successiva è il bit cotto... OneStep è una stringa di caratteri JSL di un'istruzione Python valida che imposta xNew dal valore xOld e dall'equazione di Newton, come il commento ### due righe sopra.
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 / iymaxEcco il codice Python che scrive un file binario di numeri interi e doppi. I doppi seguono gli interi per facilitare a JMP l'acquisizione di tutto con blobToMatrix.
# output results
fileobj = open(filename, mode='wb') # NoInterleave maybe prepend a descriptor
result.tofile(fileobj) # NoInterleave append after this constant size block
# NoInterleave append roundedComplexToBasinId and basinIdToBasinSize rather than printing them, they should be same dimension
# 4 8-byte doubles: real,imag,basin,size
temp = np.zeros( (len(roundedComplexToBasinId),4),dtype=np.float64)
itemp = 0
for iz in roundedComplexToBasinId.keys():
temp[itemp,0] = iz.real
temp[itemp,1] = iz.imag
temp[itemp,2] = roundedComplexToBasinId[iz]
temp[itemp,3] = basinIdToBasinSize[roundedComplexToBasinId[iz]]
itemp += 1
temp.tofile(fileobj)
fileobj.close()
print("sentinel"); # so JSL can test if this ran to the end...
def init():
offset = 0
step = 1
xdim = 3
ydim = 2
top = -1.0
bottom = 1.0
left = -1.0
right = 1.0
filename = "temp.binary"
options, remainder = getopt.getopt(sys.argv[1:], 's:x:y:o:t:b:l:r:f:')
for opt, arg in options:
if opt == '-o':
offset = int(arg) # offset is the starting point
elif opt == '-s':
step = int(arg) # step to the next point
elif opt == '-x':
xdim = int(arg) # xdim field int size
elif opt == '-y':
ydim = int(arg) # ydim field int size
elif opt == '-t':
top = float(arg) # top of field user coords
elif opt == '-b':
bottom = float(arg) # bot field user coords
elif opt == '-l':
left = float(arg) # left field user coords
elif opt == '-r':
right = float(arg) # right field user coords
elif opt == '-f':
filename = arg # write binary file
print(" xdim=",xdim,
" ydim=",ydim,
" offset=",offset,
" step=",step,
" top=",top,
" bottom=",bottom,
" left=",left,
" right=",right)
return offset, step, xdim, ydim, top, bottom, left, right, filename
if __name__ == "__main__":
main()
]\";Il programma Python viene salvato nella directory NFS (come lo sa Windows) con un numero univoco. Il nome del file restituito non viene utilizzato; la versione Linux del nome file viene creata dalla funzione di avvio e passata tramite Plink. (Siamo tornati in JSL ora...)
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);
);
Il ciclo sottostante ha generato i 10.000 fotogrammi.
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)La formula e la conversione di OneStep in un'istruzione Python: sembra complicato perché sostituisce i valori effettivi nei parametri F1...F6 e E1...E6. Nascosto alla fine è il passaggio da ^ di JSL a ** di Python, che viene eseguito nel testo, non nell'espressione, perché JSL non ha modo di rappresentare ** come espressione. F1...E6 sono gestiti come espressioni. Sin( imagenumber / 200 ) è circa 10 onde sinusoidali complete da -1 a 1. Le onde sinusoidali hanno frequenze diverse a causa dei diversi divisori. I fattori si spostano agevolmente da -2 a 2 e gli esponenti si spostano agevolmente attorno al loro obiettivo 1...6 di +/- 1. Quindi l'equazione da un'immagine all'altra è piuttosto simile, ma varia in modo un po' imprevedibile.
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 pythonL'array associativo imgNumToLabel può essere utilizzato successivamente per etichettare le immagini. Semplifica xVecchio->X...
imgNumToLabel[imagenumber]=substitute(regex(oneStep,"xNew = xOld - \((.*?)\) / ","\1")," ** ","^"," * ","·","xOld","X");
write("\!n",imgNumToLabel[imagenumber]);
//"/mnt/nfsc/farm/binary" || char(imagenumber) || ".generated"Il test per un file esistente è stato utilizzato durante lo sviluppo per velocizzare le cose quando il programma doveva essere riavviato. Potrebbe non funzionare correttamente, non è stato utilizzato per la corsa finale (tre giorni!)
if(fileexists("\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated"),
write("Exists ","\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated");
continue();
);
makePythonCode(imagenumber,oneStep);Sembra un codice abbandonato fuori posto.
n = 0; // eventually this becomes the offset // NoInterleave offset is always 0
sentinels = 0; // NoInterleave "sentinel" might be all it returns
//write("Launching... ");È qui che iniziano i 10.000 RunPrograms. 50 o 60 funzioneranno contemporaneamente utilizzando i valori condat, molto sopra.
launch( imagenumber );
write("Launched. ")
);Il riepilogo ... aspetta che finiscano gli ultimi RunPrograms, unstick potrebbe trovarne uno o due.
while( N Items( imgNumToRP ) > 0,
Wait( .1 );
unstick();
);
write("\!ndone");Le etichette vengono salvate per il generatore di immagini, in seguito. JSL/Python sopra dati binari già salvati che devono ancora essere convertiti in formato .png.
savetextfile("$desktop/imgNumToLabel.txt",char(imgNumToLabel));Il Rpdone detiene 10.000 RunProgram terminati. Lasciarlo andare.
imgNumToRPdone = 0;// clean out // NoInterleave picture generation happens after all count and basin are created
Il generatore di immagini può essere eseguito ora o successivamente...
include("pythonFarmDisplay.jsl");Questa è la fine del codice di Beowulf. Il generatore di immagini avrebbe potuto usarlo, ma non l'ha fatto. Voglio che il video utilizzi gli stessi colori da un fotogramma all'altro. I colori rappresentano i bacini di attrazione, lo zero da cui è attratto qualsiasi punto del bacino dovrebbe avere lo stesso colore. Ma l'equazione cambia un po' a ogni fotogramma, quindi gli zeri si spostano. Il generatore di immagini tenta di far corrispondere gli zeri dell'immagine precedente con l'immagine successiva. A volte saltano, ma è abbastanza buono.
Segue la JSL per fare le foto. Per prima cosa crea una tavola dei colori per i bacini. Queste sono la versione più leggera possibile del colore; vengono oscurati in seguito dal numero di iterazioni elevato.
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;
);Le etichette immagine salvate in precedenza devono essere ricaricate, le dimensioni dell'immagine devono corrispondere. Questo codice dovrebbe essere lo stesso di prima e dovrebbe essere stato estratto in un file di parametri comune. Troppo tardi adesso!
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;
I colori per i bacini dovrebbero seguire immagini/equazioni leggermente diverse. Qui è dove vivono gli zeri dell'immagine precedente.
prevzeros = []; // what came before the first one
Il disco NFS era troppo piccolo. Ci sono 309 GB di informazioni sui pixel non compressi qui.
//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 );
I file binari verranno decompressi in due matrici JSL: un conteggio del numero di iterazioni necessarie per il pixel e un piccolo numero intero che identifica il bacino del pixel. Ci sono un po' più di dati binari dopo (i byte di bilanciamento) che portano il numero del bacino a una coordinata complessa; la coordinata complessa è ciò che viene abbinato da un'immagine all'altra.
For( imagenumber = 1, imagenumber < 10000, imagenumber += 1,
file = Eval Insert( filepat );
If( File Exists( file ),
bindata = Load Text File( file, blob() );
matdata = Blob To Matrix( bindata, "uint", 2, "little", 2 );
countArray[0] = matdata[1 :: yresolution * xresolution, 1];
basinArray[0] = matdata[1 :: yresolution * xresolution, 2];
balanceBytes = Length( bindata ) - 4 * yresolution * xresolution;
zeros = Blob To Matrix( Blob Peek( bindata, Length( bindata ) - balanceBytes ), "float", 8, "little", 4 );
//Show( zeros );
if(zeros[1, 3]==2,
zeros = [9e99 9e99 1 0] |/ zeros; // python code should have put this in
);
// re-assign the basin numbers to line up with prevzeros
If( N Rows( prevzeros ) == 0,
prevzeros = zeros
);
KDTable trova i vicini più prossimi; il piano complesso è un bel quartiere.
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] );
);
);Il <
/*
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,La nuova immagine ha più zeri della vecchia immagine, mettili nell'ultima bacinella.
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
Gli zeri sono ordinati sopra; la riga 1 ottiene il colore 1, ecc.
prevzeros = zeros;
Le immagini con conteggi elevati vengono scartate per rendere il video più fluido.
countarray[loc(countarray>200)] = 200; // 150 is typical big count, 999 is not converged
rangecount = Range( countarray );
countarray=countarray-3; // 3 or 4 is typical minimum count, make it bright (0)
countarray[loc(countarray<0)]=0;
GrayArray è l'oscurità scalata dai conteggi
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
Il gray * redmap[basinarray] (etc) è un componente di colore rgb in scala
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;
Le correzioni per l'overflow/underflow dei valori di colore 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;
I marcatori + vengono aggiunti nelle posizioni zero ... grezzo, ma funziona. Non disegnare sul bordo, l'indicizzazione della matrice esegue intere righe o colonne per l'indice zero!
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];L'etichetta viene aggiunta in basso. Crea una piccola bitmap dell'etichetta, quindi copiala nell'immagine più grande.
if(1,
//label = char(imagenumber)||" zeros="||char(N Rows( zeros ))||" "||label;
labelpic = textbox(label,setwrap(3000),<<SETFONTSIZE></SETFONTSIZE>
L'immagine è fatta qui. Quindi decidi se verrà mantenuto nel video nominandolo.
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>Questo è tutto! Ecco il filmato. Ho deciso di realizzare solo la versione scientifica, con fotogrammi scuri tremolanti dove si verificano e una didascalia fitta sotto. Una bella versione potrebbe essere più fluida e ancora più ipnotica. È un video 4K, se hai l'hardware e la velocità di connessione, altrimenti YouTube lo ridimensionerà per te.
Mancano due ore del 2020 mentre pubblico questo. Non vedo davvero l'ora che arrivi il 2021.
