Un giorno, nel 1971, il signor Hanson disse: “Oggi ti insegnerò come ottenere la radice quadrata di un numero. Un modo più veloce del lunga divisione modo in cui dovrei insegnare. Ti insegnerò con una parola. E non dimenticherai mai". Poi scrive alla lavagna
Indovina
e lo sottolinea. Quando la risata si è spenta, ha spiegato di dividere l'ipotesi nel numero originale, media, ripetizione. Probabilmente ha menzionato Newton-Raphson .
Il metodo di Newton per trovare le radici o gli 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.162277766016838 – 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 Derivative funziona alla grande su polinomi semplici, dando 2X in questo esempio. NameExpr() ottiene la formula da f senza valutarla e Expr() esegue l'esecuzione di 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 passaggi vengono eseguiti, stampati e contati. Il test sta verificando la differenza tra il valore X attuale e quello successivo; se il valore è molto piccolo alla dimensione attuale di X, probabilmente non si avvicinerà e potrebbe oscillare tra due numeri molto vicini, quindi interrompe il 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. Partendo da -8 si è superato a -10,25, poi si è stabilizzato a -10.
(Perché si divide per la derivata, f' ? Intuitivamente, quando la pendenza della linea f è vicina a zero, il punto successivo dovrà essere più lontano per modificare il risultato Y di f, e dividendo per la pendenza lo fa.)
Il JSL sopra non gestisce bene le richieste impossibili. Se viene chiesto di trovare sqrt(-100), il JSL esegue il loop, 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 semilavorati. Uno è basato su classi JSL e uno basato su namespace 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 JSL precedente. 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 magnitudine() ma corrisponde alla versione non complessa come abs, qualche dispari -guardando parametri extra per add() e sub() e una funzione convertExpression(). Grande ringraziamento a Stackexchange/Michael Hardy per una funzione pow() facile da afferrare!
La funzione convertExpression esegue il massaggio di 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. Ecco perché ci sono argomenti di opzione aggiuntivi 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 iniziale x e restituisce un elenco 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 iniziale 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 risposta immaginaria pura 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 ha trovato uno basato sul punto di partenza. L'altro è a -10i.
Modifichiamo il codice per tenere traccia del conteggio dell'iterazione n per un gruppo di punti nel piano complesso e grafico che, forse come diagramma di contorno, e utilizziamo cubico anziché 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 del piano complesso che ripetono la 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 per convergere. Le macchie grigie non convergevano dopo 200 iterazioni:
Three basins. Larger dots where the basins touch take longer to converge.Tre bacini. I punti più grandi in cui i bacini si toccano impiegano più tempo a convergere.
Il raddoppio della 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 ei parametri necessari. Il programma Python viene eseguito, quindi archivia i risultati su un disco NFS, quindi esce. Quando JMP vede il RunProgram terminato, raccoglie i risultati e avvia il RunProgram successivo. Dopo 10.000 RunPrograms 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.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 delle dimensioni 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 sta pubblicizzando 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 di condat e converte i nomi come dell9010 nel loro indirizzo IP puntato (o equivalente IPV6). Previene errori di ricerca in seguito; 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." );Il prossimo passo riguarda le credenziali: ppk è la chiave pivate della macchina JMP. La parte pubblica di quella chiave si trova su ogni macchina nel cluster. Plink ha bisogno di quel bit privato per firmare la connessione SSH alla macchina remota. Il codice if(0) è lì, quindi posso sperimentare su tutte le macchine. Alla fine della giornata, potrei inviare loro tutti il comando di spegnimento. Fai attenzione con il comando di spegnimento se l'elenco di 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" );
);Il prossimo passo 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 il programma Python dovrebbe lavorare. ImgNumToRP contiene un oggetto RunProgram in esecuzione. Quell'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 più avanti, 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 il termine di RunProgram.
unstick = function({},{imagenumber},
for( imagenumber = imgNumToRP << first, !Is Empty( imagenumber ), imagenumber = imgNumToRPnext,
imgNumToRPnext = imgNumToRP << next( imagenumber );
if( (imgNumToRP[imagenumber]) << canread == 0 & (imgNumToRP[imagenumber]) << canwrite == 0 & (imgNumToRP[imagenumber]) << isreadeof == 1,
imgNumToRPresult[imagenumber] = "\!n" || char(imagenumber) || ": " || imgNumToRPresult[imagenumber] || "\!n"; // NoInterleave save this with count and basin
imgNumToRPdone[imagenumber] = imgNumToRP[imagenumber];
imgNumToRP << removeitem( imagenumber );
//if(imgNumToMachine[imagenumber]!=machine,throw("b machine not match"||char(machine)||","||char(imgNumToMachine[imagenumber])));
condat[imgNumToMachine[imagenumber]][running] = condat[imgNumToMachine[imagenumber]][running] - 1;
imgNumToMachine[imagenumber] = 0;
try(
write("\!nunsticking ",imagenumber);
process( imagenumber ); // must be after...
,
show(exception_msg);
);
Break(); // restart with <<FIRST></FIRST>La funzione findmachine di seguito cerca una macchina che non ha avviato troppe sessioni SSH. 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 "sentinella" alla fine del loro output e questo codice lo verifica come un controllo grezzo che tutto funzioni senza intoppi. Bindata è il file binario che i programmi Python producono nella directory della farm sul disco NFS. Windows e JMP vedono il disco NFS nel percorso \\VBOXSVR ma Linux e i 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, usato due volte qui. 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ù sotto. Il lancio sulla linea 5 in basso è un controllo di integrità; il JSL è abbastanza complicato, mi sentivo come se dovessi scrivere del codice extra per cercare 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 parametro rende disponibile il suo argomento alla funzione read. Questo è il modo in cui puoi 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 terminazione ed è ancora in esecuzione durante l'esecuzione di 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 integrità (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 integrata. Il programma python include una funzione init() che analizza la riga di comando per ottenere le dimensioni dell'immagine e le informazioni sul nome del file di output. I dati vengono passati in due modi: parametri della riga di comando e equazioni integrate. Probabilmente tutte le informazioni potrebbero essere state integrate, 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 RunPrograms se era completamente controllato dai parametri della riga di comando. Il codice Python è scritto in modo che possa fare una parte del problema se volevo che più computer lavorassero su una singola immagine. Nel complesso è più veloce eseguire un'immagine completa in una 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 rendere facile per JMP prendere 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 conosce 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 c'è il cambiamento da ^ di JSL a ** di Python, che viene fatto nel testo, non nell'espressione, perché JSL non ha modo di rappresentare ** come espressione. Gli F1...E6 vengono gestiti come espressioni. Sin( imagenumber / 200 ) è circa 10 onde sinusoidali complete da -1 a 1. Le onde sinusoidali sono frequenze diverse a causa dei diversi divisori. I fattori si spostano senza intoppi da -2 a 2 e gli esponenti si spostano senza intoppi attorno al loro target 1...6 di +/- 1. Quindi l'equazione da un'immagine all'altra è abbastanza 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 in seguito 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 operazioni quando è stato necessario riavviare il programma. 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 verranno eseguiti contemporaneamente utilizzando i valori condat, molto sopra.
launch( imagenumber );
write("Launched. ")
);La conclusione... attendi che gli ultimi RunPrograms finiscano, 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 un secondo momento. JSL/Python sopra i dati binari già salvati che devono ancora essere convertiti in formato .png.
savetextfile("$desktop/imgNumToLabel.txt",char(imgNumToLabel));Il Rpdone contiene 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 in un secondo momento...
include("pythonFarmDisplay.jsl");Questa è la fine del codice Beowulf. Il generatore di immagini avrebbe potuto usarlo, ma non l'ha fatto. Voglio che il video usi gli stessi colori da un fotogramma all'altro. I colori rappresentano bacini di attrazione, lo zero da cui è attratto qualsiasi punto del bacino dovrebbe assumere lo stesso colore. Ma l'equazione cambia, un po', su ogni fotogramma, quindi gli zeri si spostano. Il generatore di immagini cerca di far corrispondere gli zeri dell'immagine precedente con l'immagine successiva. A volte saltano, ma è abbastanza buono.
Segue il JSL per realizzare le immagini. Per prima cosa crea una tabella dei colori per i bacini. Queste sono la versione più chiara 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 dell'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;
Si suppone che i colori per i bacini seguano immagini/equazioni leggermente diverse. È qui che risiedono 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 di quante iterazioni sono necessarie al pixel e un piccolo numero intero che identifica il bacino del pixel. Ci sono un po' più di dati binari dopo quello (i byte di equilibrio) 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
);
Il KDTable trova i vicini più vicini; 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 foto ha più zeri della vecchia foto, mettili nell'ultimo bacino.
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;
Il 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 grigio * redmap[basinarray] (ecc) è 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 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 segnalini + vengono aggiunti nelle posizioni zero... grezzo, ma funziona. Non disegnare sul bordo, l'indicizzazione della matrice esegue intere righe o colonne per l'indice di 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 è 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 video. Ho deciso di realizzare solo la versione scientifica, con fotogrammi scuri tremolanti dove accadono e una didascalia impegnativa sotto. Una versione carina potrebbe essere più liscia 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 al 2020 mentre pubblico questo. Non vedo l'ora che arrivi il 2021.
