Un día en 1971, el Sr. Hanson dijo: “Hoy les voy a enseñar cómo obtener la raíz cuadrada de un número. Una forma más rápida que la división larga forma en que se supone que debo enseñar. Te voy a enseñar con una palabra. Y nunca lo olvidarás. Luego escribe en la pizarra.
Adivinar
y lo subraya. Cuando la risa se calmó, explicó acerca de dividir la conjetura en el número original, promedio, repetición. Probablemente mencionó Newton-Raphson .
El método de Newton para encontrar raíces o ceros de una ecuación usa la derivada de la ecuación para acelerar el proceso. La raíz cuadrada de 10 parece encontrar un valor de X que haga que X^2 == 10 o X^2 – 10 == 0. La derivada de X^2 – 10 es 2X, y el corazón del ciclo para encontrar valores sucesivos de X es
Xnuevo = X - ( (X^2 - 10) / (2X) )
que, si juegas con él, es exactamente el adivinar fórmula, y converge así, para una suposición inicial de 2:
2 – suposición inicial
3.5 – aplicar x = (x + 10/x)/2
3.17857142857143
3.16231942215088
3.16227766044414
3.16227766016838
3.16227766016838: ¡seis iteraciones para una respuesta estable!
Aquí hay un programa JSL para demostrar.
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 primera línea define una expresión, f, que ayudará a encontrar la raíz cuadrada de V. La segunda línea, paso, es la iteración de Newton-Raphson para un ciclo; step aparece en el ciclo while donde se evalúa al siguiente valor X. La función Derivada funciona muy bien en polinomios simples, dando 2X en este ejemplo. NameExpr() obtiene la fórmula de f sin evaluarla y Expr() hace que NameExpr() se ejecute (para obtener la fórmula) cuando EvalExpr() encuentra Expr(). El paso se establece en x - f(x)/f'(x).
Dentro del bucle while(1), se toman, imprimen y cuentan los pasos. La prueba está comprobando la diferencia de los valores X actuales y siguientes; si el valor es muy pequeño en el tamaño actual de X, probablemente no se acerque y podría oscilar entre dos números muy cercanos, por lo que sale del bucle while.
0 -10.25
1 -10.0030487804878
2 -10.0000004646115
3 -10
hecho, una raíz de x ^ 2 - V para V = 100 es -10
¡Recuerde que hay una raíz cuadrada tanto positiva como negativa! Los valores iniciales determinan a qué raíz conduce la convergencia. Comenzando en -8 sobrepasó a -10.25, luego se instaló en -10.
(¿Por qué se divide por la derivada, f'? Intuitivamente, cuando la pendiente de la línea f es cercana a cero, el siguiente punto deberá estar más lejos para cambiar el resultado Y de f, y dividir por la pendiente hace eso).
El JSL anterior no maneja bien las solicitudes imposibles. Si se le pide que encuentre sqrt (-100), el JSL realiza un bucle, buscando un número complejo que no está equipado para encontrar.
Números complejos
En el índice de secuencias de comandos de JSL, la búsqueda de complejos produce dos ejemplos diferentes a medias. Uno se basa en clases JSL y otro en espacios de nombres y matrices. Creo que la versión de matriz está más cerca de una solución robusta y probablemente más rápida. La versión de clase es tentadora, pero tuve problemas para usarla (velocidad, problemas de memoria). Dado que JMP no tiene sobrecargas de usuarios para los operadores, ninguno será fácil de usar.
Este JSL es muy similar al JSL anterior. Hay una gran cantidad de código para manejar la aritmética de números complejos, la misma fórmula, más o menos, y ahora también maneja 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() )
)
)
}
);Arriba se agregó una función pow() para x ^ 2 o x ^ 3, etc., donde x es un número complejo, una función abs () que tal vez debería llamarse magnitud () pero coincide con la versión no compleja como abs, algo extraño -buscando parámetros adicionales para agregar() y sub(), y una función convertExpression(). ¡Un gran saludo a stackexchange/Michael Hardy por una función pow() fácil de entender!
La función convertExpression transforma una expresión JSL con los operadores +, -, *, / y ^ en una expresión equivalente que llama a las funciones de espacio de nombres complejo. Los operadores binarios de JSL son bastante agresivos al convertir ((a+b)+(c+d)) en add(a,b,c,d) que no es una función binaria... tiene cuatro argumentos. Es por eso que hay argumentos de opción extra en las funciones complex:add y :sub.
El siguiente bit de JSL mueve el ciclo a una función, anticipando el desarrollo 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 toma un punto de partida x y devuelve una lista con la solución iterada x y el número de iteraciones n.
Y el código para llamar a la función... la fórmula y la derivada son las mismas para los números complejos.
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 tercera línea llama a la función convertExpression para cambiar los operadores +, -, *, /, ^ de las versiones de números reales a las versiones de números complejos. V se establece en un valor complejo y el punto de partida inicial se establece en un punto fuera de la recta numérica real. Cuando haya terminado, el registro muestra
hecho, una raíz de x ^ 2 - V para V = -100+0i es 0+10i
[0 0]
V=-100 + 0i es el número real -100 (sin componente imaginario). 0 + 10i es la respuesta imaginaria pura sin componente real; 10i * 10i es 100*(i^2), o -100, como se esperaba. Las raíces complejas ocurren en pares conjugados; El método de Newton encontró uno basado en el punto de partida. El otro está a -10i.
Modifiquemos el código para realizar un seguimiento del recuento de iteraciones n para un grupo de puntos en el plano complejo y grafique eso, tal vez como un gráfico de contorno, y use cúbico en lugar de cuadrático:
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>Todos los puntos en el plano complejo que iteran a la misma raíz forman una cuenca. Arriba, el generador de gráficos combina los valores reales e imaginarios en un solo valor para identificar una raíz.
La imagen muestra las tres cuencas en rojo, cian y azul oscuro. Todos los puntos de inicio rojos convergen al mismo valor. Los tamaños de punto más grandes tardan más en converger. Los puntos grises no convergieron después de 200 iteraciones:
Three basins. Larger dots where the basins touch take longer to converge.Tres cuencas. Los puntos más grandes donde se tocan las cuencas tardan más en converger.
Duplicar el tamaño de la imagen de 16 000 a 65 000 puntos lleva unos 10 minutos con la implementación de números complejos de JSL. Quiero hacer 10.000 mapas de bits de puntos de 4Kx2K, así que necesito otro enfoque.
Cúmulo de Beowulf
A Cúmulo de Beowulf es una pila de computadoras viejas, en su mayoría sin pantallas ni teclados, pero conectadas en una LAN a una computadora de control. En mi caso, la computadora de control usa JSL para decirle al clúster qué hacer y recopilar los resultados. El clúster basado en Linux ejecuta programas de Python que crea JMP. La computadora con Windows que ejecuta JMP usa PLINK de RunProgram para abrir una o más sesiones SSH en cada computadora y proporcionar el programa y los parámetros necesarios. El programa Python se ejecuta, luego almacena los resultados en un disco NFS y luego sale. Cuando JMP ve que RunProgram termina, recopila los resultados e inicia el siguiente RunProgram. 10,000 RunPrograms más tarde hay suficientes imágenes para hacer un video de cinco minutos.
Ethernet hub in center, below power-strip. Not enough mice and monitors and keyboards, SSH is best.Concentrador Ethernet en el centro, debajo de la regleta de alimentación. No hay suficientes ratones, monitores y teclados, SSH es lo mejor.
Python tiene un tipo de datos de número complejo nativo, por lo que parte del JSL anterior desaparece. El JSL ahora se ve así; verá algo de Python incrustado como una cadena en un momento.
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;
El bit anterior incluye la ecuación y algunos parámetros de tamaño de video. El bit a continuación describe la red. Hostkey, server, cpus, running, gigs son índices en la lista condat (datos de conexión) de información sobre el clúster.
// 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...
El problema con el daemon aún existe, pero ahora sé que anuncia las máquinas con un sufijo -2 (etc.) y puedo solucionarlo sin conectar una pantalla.
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
};El JSL a continuación revisa la lista de datos y convierte los nombres como dell9010 en su dirección IP punteada (o equivalente IPV6). Previene fallas de búsqueda posteriores; el clúster funcionó durante unos tres días y todavía no estoy seguro de por qué el escritorio de mi esposa proporciona el servicio de búsqueda de DNS... Antes de recordar GetAddrInfo, usé el código de ping comentado, que es un poco más torpe. Mi casa está mucho más feliz con IPV6 ahora que hace unos años. En algún momento, la dirección IPV6 pertenece entre corchetes.
// 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." );El siguiente bit es sobre las credenciales: el ppk es la clave privada de la máquina JMP. La parte pública de esa clave está en cada máquina del clúster. Plink necesita ese bit privado para firmar la conexión SSH a la máquina remota. El código if(0) está ahí para que pueda experimentar en todas las máquinas. Al final del día, podría enviarles a todos el comando de apagado. Tenga cuidado con el comando de apagado si la lista de condatos contiene la máquina Linux que ejecuta VirtualBox que ejecuta este script JSL (línea 9 anterior).
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" );
);El siguiente bit se trata de ejecutar muchas copias de RunProgram al mismo tiempo. ImgNum es la clave para cada una de estas matrices asociativas y es el número de la imagen en la que debería estar trabajando el programa de Python. ImgNumToRP contiene un objeto RunProgram que se está ejecutando. Ese objeto se mueve a imgNumToRPdone cuando finaliza. imgNumToResult contiene la salida de texto del programa Python, imgNumToMachine la identificación de la máquina e imgNumToLabel una etiqueta para la imagen. (Las etiquetas se aplican más adelante, más abajo).
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 función unstick maneja los RunPrograms que no llaman a ReadFunction por última vez después de que 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 siguiente función findmachine busca una máquina que no tenga demasiadas sesiones SSH iniciadas. La función unstick() se llama periódicamente para buscar el 10% de los casos que se atascan.
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 función de proceso se llama cuando JMP detecta que RunProgram ha finalizado. Los programas de Python escriben "centinela" al final de su salida y este código lo verifica como una verificación cruda de que todo funciona sin problemas. Bindata es el archivo binario que los programas Python producen en el directorio de la granja en el disco NFS. Windows y JMP ven el disco NFS en la ruta \\VBOXSVR pero Linux y los programas de Python lo ven en la ruta //mnt. Este código iba a producir las imágenes a partir de los datos binarios, pero eso se separó en un paso de procesamiento posterior. Ahora solo informa la cantidad de ceros descubiertos para la fórmula. El formato binario se explica en el código python; JMP lo lee y lo carga con blobtomatrix, usado dos veces aquí. El relanzamiento al final probablemente funcione si no se comenta y es necesario; Olvidé descomentarlo y no lo necesitaba.
// 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 función de inicio inicia otro RunProgram. Utiliza findmachine para esperar una máquina disponible y proporciona un archivo .py personalizado y parámetros personalizados al comando de Python. El archivo .py se crea más abajo. El lanzamiento en la línea 5 a continuación es un control de cordura; el JSL es lo suficientemente complicado, sentí que necesitaba escribir un código adicional para detectar errores lógicos.
// 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 palabra clave de parámetro hace que su argumento esté disponible para la función de lectura. Así es como puede tener varios RunPrograms en paralelo y hacer que las funciones de lectura sepan cuál es.
Parameter( evallist({imagenumber,machine}) ),La función de lectura captura la salida e intenta (alrededor del 90 % de las veces con éxito) detectar que RunProgram ha finalizado. Debido a que RunProgram aún está terminando y aún ejecutándose cuando se ejecuta este código, el identificador se mueve a imgNumToRPDone en lugar de simplemente eliminarlo. Eso le permite terminar limpiamente. Parm, abajo, recibe el valor de Parámetro, arriba. Hay otro control de cordura (lanzamiento) a continuación.
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 función makePythonCode genera el archivo .py con la ecuación oneStep incorporada. El programa python incluye una función init() que analiza la línea de comandos para obtener información sobre el tamaño de la imagen y el nombre del archivo de salida. Los datos se pasan de dos maneras: parámetros de línea de comandos y ecuaciones integradas. Posiblemente toda la información podría haber sido integrada, pero si no hay razón (como compilar la ecuación de un solo paso en el código) para incluir parámetros en el código, solo necesitaría una copia del código para todos los RunPrograms si estaba completamente controlado por parms de línea de comando. El código de Python está escrito para que pueda hacer una parte del problema si quisiera que varias computadoras trabajaran en una sola imagen. En general, es más rápido hacer una imagen completa a la vez y no recombinar partes más tarde.
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 siguiente línea es el bit integrado... OneStep es una cadena de caracteres JSL de una declaración Python válida que establece xNew a partir del valor xOld y la ecuación de Newton, como el comentario ### dos líneas arriba.
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 / iymaxAquí está el código de Python que escribe un archivo binario de enteros y dobles. Los dobles siguen a los números enteros para facilitar que JMP tome todo 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()
]\";El programa Python se guarda en el directorio NFS (como lo conoce Windows) con un número único. El nombre de archivo devuelto no se utiliza; la versión de Linux del nombre de archivo se crea mediante la función de inicio y se pasa a través de Plink. (Estamos de vuelta en JSL ahora...)
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);
);
El ciclo a continuación generó los 10,000 cuadros.
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 fórmula y la conversión de OneStep en una declaración de Python: esto parece complicado porque sustituye los valores reales en los parámetros F1...F6 y E1...E6. Oculto al final está el cambio de ^ de JSL a ** de Python, que se realiza en texto, no en expresión, porque JSL no tiene forma de representar ** como una expresión. Los F1...E6 se manejan como expresiones. Sin (número de imagen / 200) es aproximadamente 10 ondas sinusoidales completas de -1 a 1. Las ondas sinusoidales tienen frecuencias diferentes debido a los diferentes divisores. Los factores se mueven suavemente de -2 a 2 y los exponentes se mueven suavemente alrededor de su objetivo 1...6 en +/- 1. Entonces, la ecuación de una imagen a la siguiente es bastante similar, pero varía un poco de manera impredecible.
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 pythonLa matriz asociativa imgNumToLabel se puede usar más tarde para etiquetar imágenes. Simplificar xViejo->X...
imgNumToLabel[imagenumber]=substitute(regex(oneStep,"xNew = xOld - \((.*?)\) / ","\1")," ** ","^"," * ","·","xOld","X");
write("\!n",imgNumToLabel[imagenumber]);
//"/mnt/nfsc/farm/binary" || char(imagenumber) || ".generated"La prueba de un archivo existente se usó durante el desarrollo para acelerar las cosas cuando el programa tenía que reiniciarse. Puede que no funcione del todo bien, no se usó para la ejecución final (¡tres días!)
if(fileexists("\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated"),
write("Exists ","\\VBOXSVR\nfsc\farm\binary" || char(imagenumber) || ".generated");
continue();
);
makePythonCode(imagenumber,oneStep);Esto parece un código abandonado fuera de lugar.
n = 0; // eventually this becomes the offset // NoInterleave offset is always 0
sentinels = 0; // NoInterleave "sentinel" might be all it returns
//write("Launching... ");Aquí es donde comienzan los 10.000 RunPrograms. 50 o 60 se ejecutarán al mismo tiempo utilizando los valores de condat, muy por encima.
launch( imagenumber );
write("Launched. ")
);El resumen... espere a que finalicen los últimos RunPrograms, despegar puede encontrar uno o dos.
while( N Items( imgNumToRP ) > 0,
Wait( .1 );
unstick();
);
write("\!ndone");Las etiquetas se guardan para el generador de imágenes, más tarde. JSL/Python anterior ya guardó datos binarios que aún necesitan conversión al formato .png.
savetextfile("$desktop/imgNumToLabel.txt",char(imgNumToLabel));El Rpdone tiene 10.000 RunPrograms terminados. Déjalo ir.
imgNumToRPdone = 0;// clean out // NoInterleave picture generation happens after all count and basin are created
El generador de imágenes se puede ejecutar ahora o más tarde...
include("pythonFarmDisplay.jsl");Ese es el final del código Beowulf. El generador de imágenes podría haberlo usado, pero no lo hizo. Quiero que el video use los mismos colores de un cuadro al siguiente. Los colores representan cuencas de atracción, el cero que atrae a cualquier punto de la cuenca debe tener el mismo color. Pero la ecuación cambia, un poco, en cada fotograma, por lo que los ceros se mueven. El generador de imágenes intenta hacer coincidir los ceros de la imagen anterior con la imagen siguiente. A veces saltan, pero es lo suficientemente bueno.
El JSL para hacer las fotos sigue. Primero haz una tabla de colores para los lavabos. Estas son la versión más clara posible del color; se oscurecen más tarde por el gran número de iteraciones.
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;
);Las etiquetas de imagen guardadas anteriormente necesitan recargarse, las dimensiones de la imagen deben coincidir. Este código debería ser el mismo que el anterior y debería haberse extraído a un archivo de parámetros común. ¡Muy tarde ahora!
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;
Se supone que los colores de las cuencas siguen imágenes/ecuaciones ligeramente diferentes. Aquí es donde viven los ceros de la imagen anterior.
prevzeros = []; // what came before the first one
El disco NFS era demasiado pequeño. Aquí hay 309 GB de información de píxeles sin comprimir.
//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 );
Los archivos binarios se desempaquetarán en dos matrices JSL: un recuento de cuántas iteraciones necesitó el píxel y un pequeño número entero que identifica la cuenca del píxel. Hay un poco más de datos binarios después de eso (los bytes de balance) que llevan el número de cuenca a una coordenada compleja; la coordenada compleja es lo que se empareja de una imagen a la siguiente.
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 encuentra los vecinos más cercanos; el complejo plano es un excelente vecindario.
kdt = KDTable( prevzeros[0, 1 :: 2] ); // 1::2 means the real and imag cols
newbasin = (1 :: N Rows( zeros )) + 10000;
toomany = N Rows( zeros ) > N Rows( prevzeros );
For( i = 1, i <= N Rows( zeros ), i += 1,
If( zeros[i, 3] != i,
Throw( "basin index assumption wrong" )
); // zeros are sorted by basin number
{indx, dist} = kdt << knearestrows( 2, zeros[i, 1 :: 2] );
If( N Cols( dist ) > 1 & dist[1] < dist[2] / 10, // assume 10X closer is a good match
newbasin[i] = indx[1];
kdt << removerows( indx[1] );//
, /*else*/ !toomany, // if there are toomany new zeros, don't assign them this way
newbasin[i] = indx[1];
kdt << removerows( indx[1] );
);
);<
/*
zerobasin newbasin
1 4 change 1s to 4s in basinArray and zeros
2 3 2s to 3s
3 10001 too many case! no basin assigned. becomes 7, below
4 2
5 5
6 6
7 1
*/
unassigned = Loc( newbasin > 10000 );
If( N Rows( unassigned ) > 0,La nueva foto tiene más ceros que la foto anterior, colóquelos en el último recipiente.
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
Los ceros se ordenan arriba; la fila 1 obtiene el color 1, etc.
prevzeros = zeros;
Las imágenes con conteos altos se descartan para que el video sea más 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;
El GrayArray es oscuridad escalada de los conteos.
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
El gray * redmap[basinarray] (etc) es un componente de color rgb escalado
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;
Las correcciones para el desbordamiento/desbordamiento de los valores de color 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;
Los marcadores + se agregan en las ubicaciones cero... crudo, pero funciona. ¡No dibuje en el borde, la indexación de matriz hace filas o columnas enteras para el índice de cero!
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];La etiqueta se agrega en la parte inferior. Haz un pequeño mapa de bits de la etiqueta, luego cópialo en la imagen más grande.
if(1,
//label = char(imagenumber)||" zeros="||char(N Rows( zeros ))||" "||label;
labelpic = textbox(label,setwrap(3000),<<SETFONTSIZE></SETFONTSIZE>
La imagen está hecha aquí. Luego, decida si se mantendrá en el video nombrándolo.
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>¡Eso es todo! Aquí está el vídeo. Decidí hacer solo la versión científica, con marcos oscuros parpadeantes donde suceden y una leyenda ocupada debajo. Una versión bonita podría ser más suave e incluso más hipnótica. Es un video 4K, si tiene el hardware y la velocidad de conexión; de lo contrario, YouTube lo reducirá para usted.
Quedan dos horas de 2020 mientras publico esto. Tengo muchas ganas de 2021.
