cancel
Showing results for 
Show  only  | Search instead for 
Did you mean: 
Browse apps to extend the software in the new JMP Marketplace
Choose Language Hide Translation Bar
Craige_Hales
Super User
Multiple precision arithmetic with Python

Playing with JMP + Python + gmpy2 (python library for multiple precision arithmetic) + mpfr (multiple precision real) + mpc (multiple precision complex). Using the same Beowulf cluster descibed in Beowulf, Newton, and Mr Hanson . It is enormously slow, even with 8 machines/~44 CPUs available. In the video you can see a window monitoring 100 threads in a 10x10 grid; each thread is computing the pixels for a small patch of a zoomed in Mandelbrot. With normal (fast) arithmetic, getting past a zoom of 1e-13 or so is impossible because there are only about 15 significant digits. In the video, at the end, the zoom is around  1e-21. It is getting really slow because the iteration count is 200,000 for the hardest pixels, and probably close to that for most of them. Thanks to online clock for the real time clock in the bottom right.

Here's the JSL I'm playing with; it is set up for my environment and can't possibly "just work" but it might be useful for an example of

  • many RunPrograms working in parallel
  • using a cluster of Linux machines running Python
  • using the Python Multiple Precision library (you might not need that, apparently Python natively supports really big integers.)
  • putting thumbnails in a data table (better places to see that than here maybe)

The multiple precision numbers are represented as strings in JMP. The Python program has a tiny calculator JMP uses when zooming in; JMP tells the calculator the 0.0 to 1.0 coordinates of a selection rectangle in the graph box (no axes visible, but still a graph box) and Python  does the mpfr calculation and returns the long string of resulting digits.

 


// 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

condat = { // load balence with audacity running on thinkb and using 100...good enough...
// the 3,3,4,4,4,5,6,8 balance is not just about machine speed, but about the
// poor algorithm for asigning machines that uses up line 1, then 2, ... then 8
// then sets back to line 1. for 100 that means 
// 3+3+4+4+4+5+6+8+ 3+3+4+4+4+5+6+8 +3+3+4+4+4+5+3
// 9 9 12 12 12 15(b) 17(c) 19(d)
// then 'think' machines are a little faster than the 9010/7010 machines
// and twice as fast as the 790 machines. (not 8/3, but maybe 17/9 or 19/9)
{"4e:c3:8a:3b:e7:d4:bd:0d:85:62:36:40:17:15:b4:a3", "dell790a.local", 3, 0, 22},
{"00:c2:37:09:dc:d7:cf:fb:33:e3:14:e3:d6:fb:39:be", "dell790b.local", 3, 0, 14},
{"54:c8:1a:6a:50:e3:51:02:ed:1b:f8:fa:bd:f4:e4:17", "dell7010a.local", 4, 0, 30},
{"66:0d:ca:3c:1f:11:eb:45:d4:92:d0:0b:ae:01:b8:6b", "dell7010b.local", 4, 0, 30},
{"94:84:2f:a5:56:2d:db:16:b5:80:a9:0d:87:07:ee:53", "dell9010.local", 5, 0, 30},
{"b2:b7:b3:c3:83:63:11:e7:56:e7:c9:c2:e6:fc:f0:c9", "thinkb.local", 6, 0, 30},
{"2b:ca:fc:45:f0:64:cd:ac:70:ce:93:73:e8:84:2d:75", "thinkc.local", 6, 0, 30},
{"62:b7:67:23:8d:86:fe:80:1a:62:ec:d8:41:1b:b4:1b", "thinkd.local", 6, 0, 30}
//{"9a:e1:a7:05:ed:f7:d0:4c:1d:c6:cb:d2:db:f3:02:c4", "d790main.local", 0, 0, 30}, // host for vm
//{"14:26:bc:41:48:ef:43:31:61:29:e5:a6:3f:d5:43:07", "closet.local", 0, 0, 14} // tv capture
};
machinename = {};
for(i=1,i<=nitems(condat),i+=1,machinename[i]=regex(condat[i][server],"[^.]+"));

good = 0;
For( iCondat = 1, iCondat <= N Items( condat ), iCondat += 1,
    condat[iCondat][running] = 0;
    ai = Get Addr Info( condat[iCondat][server] );
    For( j = 3, j <= N Items( ai ), j += 1,
        If( "PF_INET" == ai[j][1], // needs [ ] wrapped around "PF_INET6" addrs
            Try(
                ip = ai[j][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] ) )
            )
        )
    );
    Wait( .01 );
);
If( good != N Items( condat ),
    Throw( "did not find all condat" );
    Wait( .1 );
);
Write( "\!n", good, " of ", N Items( condat ), " names resolved." );

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 -haltr /mnt/nfsc/hi"; // fstab: closet.local:/mnt/nfs /mnt/nfsc/ nfs defaults 0 0
    //command = " \!"echo `hostname`>>/mnt/nfsc/hi;cat /mnt/nfsc/hi"; // 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 = Run Program(
            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" );
);

pythonCodeFileName = "DeleteMeFrom_" || Regex( Include File List()[1], "([^\!\\!\]*?\.jsl)", "\1" ) || "_.py";
pythonCodeFilePathWin = "\\VBOXSVR\nfsc\farm\";
pythonCodeFilePathLinux = "/mnt/nfsc/farm/";

xdimension = 1920/2;
ydimension = 1080/2;
squaresize = 10; // squaresize^2 processes will run (30 works, slow to launch)
subImageWidth = xdimension / squaresize; // needed to unpack .pbd
if(subImageWidth != floor(subImageWidth),throw("subImageWidth "||char(subImageWidth)));
subImageHeight = ydimension / squaresize; // maybe not needed
if(subImageHeight != floor(subImageHeight),throw("subImageHeight "||char(subImageHeight)));


edges = [=>]; // holds the x/y lef/rig/top/bot mpfr strings for labeling graph xlef
progressbarlength = 50; // in the monitor window
resultWindowScale = 1; // scale preview
digitPrecision = 25; // prec
drillDepth = 800; // deep
colorpower = 1.3;

savetextfile(pythonCodeFilePathWin||pythonCodeFileName,"\[import threading, time, sys # , tty, termios
#import sys, getopt, time
import math
import numpy as np
import gmpy2
from gmpy2 import mpfr
from gmpy2 import mpc
from PIL import Image


eventWorkAvailable = threading.Event()
eventWorkAvailable.clear()

eventWorkCaptured = threading.Event()
eventWorkCaptured.clear()

# all parm names must be exactly 4 chars...simple parser below...

parms = { # all nums kept as strings because the parser is very simple
    'xcen': "0.0", # left edge 
    'xrad': "0.0", # radius, xmax = xmin + 2 * xrad
    'ycen': "0.0", # top edge
    'yrad': "0.0",
    'xpix': "0", # bitmap size
    'ypix': "0", # bitmap size
    'deep': "0", # depth...max iterations
    'xlef': "", # edges of whole pic (not just my bit)
    'xrig': "",
    'ytop': "",
    'ybot': "",
#    'redp': "", # rgb exponents
#    'grnp': "",
#    'blup': "",
    # next 3 calculate a new xcen/ycen/xrad/yrad, in place, from current values
    'xzfr': "0.0", # for the zoom command, 0..1 of xmin..xmax
    'yzfr': "0.0", # for the zoom command, 0..1 of ymin..ymax
    'zfra': "0.0", # typ .5 to 2.0, the zoom in/out fraction
    'name': "/mnt/xxx/abc001.pbd", # output file name
    'ncpu': "100", # number of workers total, must be a square 1,4,16,...100(sweet spot?),...256 may be near a windows limit for puttys
    'icpu': "0", # worker id in square grid of ncpu, 0 based, left to right starting at top (give all workers same xcen, etc for whole pic)
    'prec': "100", # precision, in base 10 digits (increases with 'deep' and smaller xrad/yrad)
    'prog': "." # progressbar?  no comma on last one. ***READ ONLY*** progress bar status
    }

gmpy2.get_context().precision = int(3.32*int(parms['prec']))

command = "" # global for main worker thread to copy
drop = 0 # global to communicate to main worker thread "drop the picture in progress asap and wait for next command"

def listener(): # listen to stdin for parms and commands
    global parms
    global command
    global drop
    global digits

    gmpy2.get_context().precision = int(3.32*int(parms['prec'])) # roughly digits # The thread gets its own context
    formatstring = '{0:.' + str(int(parms['prec'])+2) + 'G}' # save a few extra, general format

    ndrop=0
    while 1:
        line = sys.stdin.readline().strip()
        if not line: # eof detected and turns into normal quit command
            line = "quit"

        cmd = line[:4] # commands/parms are all exactly 4 chars
        op = line[4:].strip() # strip space between cmd...operand and trailing newline

        if cmd in parms: # asynchronous set a parm
            # cmd should NOT be 'prog', but any other 'xcen'...'prec' parm to be set to a value
            parms[cmd] = op # set a parameter value
            if cmd == 'prec':
                gmpy2.get_context().precision = int(3.32*int(parms['prec'])) # just this thread for now
                formatstring = '{0:.' + str(int(parms['prec'])+2) + 'G}' # save a few extra, general format

        elif cmd == "read": # asynchronous retrieve a parm, including 'prog' progressbar value
            # a leading space...often the first character arrives early and
            # gets separated from the message. This gets trimed if present.
            print(" "+op,parms[op]) # retrieve a parameter value
            sys.stdout.flush()

        elif cmd == "drop": # asynchronous
            drop = 1 # drop the current picture in progress
            ndrop+=1 
#            parms["prog"]="drp"+str(ndrop)
#            
#            print(" prog",parms["prog"]) # announce dropped pic
#            sys.stdout.flush()
            
        elif cmd == "zoom": # asynchronous
            # tiny mpfr calculator for zooming in/out ***using current prec***
            xcen = mpfr(parms["xcen"])
            xrad = mpfr(parms["xrad"])
            xmin = xcen - xrad
            xmax = xcen + xrad
            ycen = mpfr(parms["ycen"])
            yrad = mpfr(parms["yrad"])
            ymin = ycen - yrad
            ymax = ycen + yrad
            xzfr = mpfr(parms["xzfr"])
            yzfr = mpfr(parms["yzfr"])
            zfra = mpfr(parms["zfra"])
            parms["xcen"] = formatstring.format(xmin + xzfr * (xmax - xmin))
            parms["ycen"] = formatstring.format(ymin + yzfr * (ymax - ymin))
            parms["xrad"] = formatstring.format(xrad * zfra)
            parms["yrad"] = formatstring.format(yrad * zfra)

        else: # synchronous: these are interpreted below
            # quit or make. other garbage is ignored, mostly, but will ***block*** until prev make finishes
            command = cmd 
            #operand = op
            eventWorkAvailable.set() # main wakes up and starts capturing details
            eventWorkCaptured.wait() # main signals the details are captured (***block*** happens here if make in progress)
            eventWorkCaptured.clear() # main will need to set this
            if cmd == "quit": # and quit is also interpreted here so it doesn't read again
                break



# start the listener thread to process commands on stdin                        

stdiothread = threading.Thread(target=listener,args=(()))
stdiothread.start()




def mandelbrot(z,maxiter,horizon,log_horizon,log_2):
    c = z
    for n in range(maxiter):
        az = abs(z)
        if az > horizon:
            return float(n) - float(gmpy2.log(gmpy2.log(az))/log_2 + log_horizon) # https://gist.github.com/jfpuget/16eb22a3c26b6275d9dd -- the smooth log part
        z = z*z + c
#    return 0
    return float(maxiter)

def mandelbrot_set(xmin,xstep,ymin,ystep,nx,ny,maxiter):
    global parms
    global drop
    horizon = mpfr(2.0 ** 40)
    log_2 = gmpy2.log(2)
    log_horizon = gmpy2.log(gmpy2.log(horizon))/log_2
    step10 = min(xstep,ystep)/10 # radius of small triangle of slopes about the point
#    power=1
#    dx0 = mpfr(step10 * math.sin(0))
#    dy0 = mpfr(step10 * math.cos(0))
#    dx1 = mpfr(step10 * math.sin(math.pi*120/180))
#    dy1 = mpfr(step10 * math.cos(math.pi*120/180))
#    dx2 = mpfr(step10 * math.sin(math.pi*240/180))
#    dy2 = mpfr(step10 * math.cos(math.pi*240/180))
# old code
#    redp = float(parms["redp"])
#    grnp = float(parms["grnp"])
#    blup = float(parms["blup"])
#    n3 = np.empty((ny,nx,3),np.uint8)
# new code
    n2 = np.empty((ny,nx),dtype=np.float64)

    x = xmin
    for ix in range(nx):
        #time.sleep(0.001)
        if drop == 1:
            break
        parms["prog"] = str(int(100*ix/nx))+"%"
        if ix % 3 == 0:
            print(" prog",parms["prog"]) # announce %
            sys.stdout.flush()
        
        
        y = ymin
        for jy in range(ny):
            if drop == 1:
                break
#            p0 = mandelbrot( mpc(x+dx0, y+dy0), maxiter, horizon, log_horizon, log_2)
#            p1 = mandelbrot( mpc(x+dx1, y+dy1), maxiter, horizon, log_horizon, log_2)
#            p2 = mandelbrot( mpc(x+dx2, y+dy2), maxiter, horizon, log_horizon, log_2)
#            close = (max(p0,p1,p2)-min(p0,p1,p2))/((p0+p1+p2)/3) # 1,2,3 -> 2/6   501,502,503-> 2/502  0,500,1000 -> 1000/500
#            if close > 1:
#                close = 1
#            close = close
#            d = math.atan2(p0,p1)/(math.pi/2) # -1 ... + 1
#            e = math.atan2(p1,p2)/(math.pi/2)
#            f = math.atan2(p2,p0)/(math.pi/2)
#            n3[jy,ix] = (255*close*(d+1)/2,255*close*(e+1)/2,255*close*(f+1)/2)
            
            p = mandelbrot( mpc(x, y), maxiter, horizon, log_horizon, log_2)/maxiter
# old code
#            try:
#                n3[jy,ix] = (255*p**redp,255*p**grnp,255*p**blup)
#            except:
#                print(p,type(p),redp,grnp,blup)
# new code
            n2[jy,ix] = p
            
            y = y + ystep
        x = x + xstep
    return n2

while 1: # process work...
    eventWorkAvailable.wait()
    eventWorkAvailable.clear() # tell the listener we are not available
    copycommand = command
    command = ""
    eventWorkCaptured.set() # tell the listener we've captured the details
    
    if copycommand=="quit":
        break
    
    if copycommand=="make":
        
        gmpy2.get_context().precision = int(3.32*int(parms['prec']))
        formatstring = '{0:.' + str(int(parms['prec'])+2) + 'G}'
        
        drop = 0 # if this becomes 1, drop this picture and return to waiting
#        parms["prog"]="start" # typically 1/8192 step size, this is a smaller, barely started, value
#        print(" prog",parms["prog"]) # announce pre-percent
#        sys.stdout.flush()
        
        # figure out my responsibility in the grid
        mynumber = int(parms["icpu"]) # 0..99
        totalnumber = int(parms["ncpu"]) # 100
        sqrtnumber = math.sqrt(totalnumber) # 10
        myxgridLeft = int(mynumber % sqrtnumber) # 0..9
        myygridTop = int(mynumber / sqrtnumber) # 0..9
        # make 0..1 fractions of my grid position
        myxgridLeft = myxgridLeft / sqrtnumber
        myygridTop = myygridTop / sqrtnumber
        
        # unpack the parms for the group's grid responsibility, apply my part

        xcen = mpfr(parms["xcen"]) # center and size for the group picture
        xrad = mpfr(parms["xrad"])
        xmin = xcen - xrad # for the group
        parms["xlef"] = formatstring.format(xmin) # record group picture edges
        xsize = 2 * xrad # for the group
        parms["xrig"] = formatstring.format(xmin+xsize) # record group picture edges
        xmin += myxgridLeft * xsize # for me
        xsize /= sqrtnumber # for me

        ycen = mpfr(parms["ycen"]) # center and size for the group picture
        yrad = mpfr(parms["yrad"])
        ymin = ycen - yrad # for the group
        parms["ytop"] = formatstring.format(ymin) # record group picture edges
        ysize = 2 * yrad # for the group
        parms["ybot"] = formatstring.format(ymin+ysize) # record group picture edges
        ymin += myygridTop * ysize # for me
        ysize /= sqrtnumber # for me

        # only need to say this once (every thread says it)
        print(  " xlef",parms["xlef"],
                "\nxrig",parms["xrig"],
                "\nytop",parms["ytop"],
                "\nybot",parms["ybot"],
                "\nprec",parms["prec"],
                "\ndeep",parms["deep"],
                "\nxzfr",parms["xzfr"],
                "\nyzfr",parms["yzfr"],
                "\nzfra",parms["zfra"],
                "\nxcen",parms["xcen"],
                "\nxrad",parms["xrad"],
                "\nycen",parms["ycen"],
                "\nyrad",parms["yrad"],
                "\nprog","go go go","\n") # announce edge parms
        sys.stdout.flush()
        
        # I'm only responsible for my fraction of the pixels
        width = int(int(parms["xpix"]) / sqrtnumber)
        height = int(int(parms["ypix"]) / sqrtnumber)
        
        # make the picture
        z = mandelbrot_set(xmin,xsize/width,ymin,ysize/height,width,height,int(parms["deep"]))
        if drop == 0:
#            parms["prog"]="save" # after 100%
#            print(" prog",parms["prog"]) # announce saving
#            sys.stdout.flush()

# old code returned a png file...
#            img = Image.fromarray(z)
#            img.save(parms["name"])
# new code now binary... ".pbd" extension "Python Binary Data"...I made that up.
            fileobj = open(parms["name"], mode='wb')
            z.tofile(fileobj)
            fileobj.close()
            parms["prog"]="okpic" 
        else:
            parms["prog"]="*BAD*" 

        print(" prog",parms["prog"]) # announce fail to make pic or good pic
        sys.stdout.flush()
        
# post-quit cleanup        
#print("wait for thread to join")
stdiothread.join() # wait until listener actually finishes
print("done bye")


]\");

dtThumb = New Table( "thumbnails",
    New Column( "thumb Picture", Expression, "None", Set Display Width( 88 ) ),
    New Column( "thumb Time", Format( "hr:m:s", 13, 0 ), Input Format( "hr:m:s", 0 ), Set Display Width( 74 ) ),
    New Column( "thumb Data", Expression, "None", Set Display Width( 68 ) ),
    New Column( "thumb CenterX", Character, Set Display Width( 82 ) ),
    New Column( "thumb RadiusX", Character, Set Display Width( 86 ) ),
    New Column( "thumb CenterY", Character, Set Display Width( 78 ) ),
    New Column( "thumb RadiusY", Character, Set Display Width( 80 ) ),
    New Column( "thumb Left", Character, Set Display Width( 68 ) ),
    New Column( "thumb Right", Character, Set Display Width( 67 ) ),
    New Column( "thumb Top", Character, Set Display Width( 73 ) ),
    New Column( "thumb Bottom", Character, Set Display Width( 84 ) ),
    New Column( "thumb Digits", Format( "Fixed Dec", 12, 0 ), Set Display Width( 93 ) ),
    New Column( "thumb Iteration", Format( "Fixed Dec", 12, 0 ), Set Display Width( 79 ) )
);

/*
1:aaa
2:bbb
3:dddd
4:eeee
5:fffff
6:gggggg
7:hhhhhh
8:iiiiii
here's a pattern to spatially separate the machines
defaghibfdeghiagfhdeibghfiadeghifbghi
use items from this list, in a ring, to balance the load (aaa is a weaker machine than gggggg)
and distribute the load spatially (don't use aa next to each other in the bitmap, they might be
both in a hard patch)

built by ababab is len 6, put ghi between each is len 24, put de every 6 is len 32, f every 6 is len 37
the prime 37 should work with a 16x16 square OK
*/
separatePlan = [3 4 5 1 6 7 8 2 5 3 4 6 7 8 1 6 5 7 3 4 8 2 6 7 5 8 1 3 4 6 7 8 5 2 6 7 8];
separateIndex=1;
findMachine = Function( {},
    {iCondat, n}, // find a machine that is not filled up
    separateIndex+=1;
    if(separateIndex>nitems(separatePlan),separateIndex=1);
    return(separatePlan[separateIndex]);
);

// build/replace/insert new image 
histcolor = heatcolor((0::1022)/1022,"spectral")||0;
//shortHigh = j(1,128,1);
//longHigh = j(1,256,1);
//shortLow = j(1,128,0);
//longLow = j(1,256,0);
//shortRiseFromMiddle = (128::255)/255;
//longRiseFromBottom = (0::255)/255;
//shortFallToMiddle = (255::128)/255;
//longFallToBottom = (255::0)/255;
//red   = longLow || longLow || longLow || longRiseFromBottom;
//red[1024]=0;
//green = longRiseFromBottom || shortHigh || longHigh || longHigh || shortFallToMiddle;  
//green[1024]=0;
//blue  = shortRiseFromMiddle || longHigh || longFallToBottom || shortLow || longRiseFromBottom;
//blue[1024]=.25;
//histcolor = rgbcolor( red,green,blue);

generateImage = Function( {},
    {minmaxrange, normalized, img, imgseg, hist, color}, 
    minmaxrange = Range( compositionMatrix );
    normalized = (compositionMatrix - minmaxrange[1]) / (minmaxrange[2] - minmaxrange[1]);

    hist = Sort Ascending( normalized[0] );
    hist = hist[1 :: N Rows( hist ) :: N Rows( hist ) / 1024];
    loc1 = Loc( hist, 1 );
    If( N Items( loc1 ),
        hist = hist[1 :: loc1[1]]
    ); // only one slot for max
    posterized = Round( (1024 / N Rows( hist )) * Loc Sorted( hist, normalized ^ colorpower ) ); // <<<<<<<<< trim power
    
    color = Shape( histcolor[posterized], N Rows( compositionMatrix ), N Cols( compositionMatrix ) );
    
    img = New Image( color );
    //newwindow("x",img)
    //    img = newimage("rgb",{normalized^redpower,normalized^greenpower,normalized^bluepower});
    Try(
        imgseg = resultGB[framebox( 1 )] << FindSeg( PictSeg( 1 ) );
        try(imgseg << remove); // delete old and replace with new...fail expected 1st time
        resultGB[framebox( 1 )] << Add Image( image( img ), bounds( top( 0 ), Left( 0 ), bottom( 1 ), Right( 1 ) ) );//
        imgseg = resultGB[framebox( 1 )] << FindSeg( PictSeg( 1 ) );
        imgseg << lock( 1 ); // so the mouse won't interact with the image and drag it around	
    ,//
        show(exception_msg)
    );
    return(img);
);

// launch a thread in the farm
launch = Function( {imagenumber,machine = findMachine()},
    {command},
    
    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.8 " || pythonCodeFilePathLinux || pythonCodeFileName;
    imgNumToRP[imagenumber] = Run Program(
        executable( "plink" ),
        options(
            localkey || /* -v */ " -hostkey " || condat[machine][hostkey] || " " || remoteuser || condat[machine][server] || " " || command || // 
            ""
        ),
        Parameter( Eval List( {imagenumber, machine} ) ),
        readfunction(
            Function( {this, parm},
                {t, txt, txttrim, imagenumber, machine},
                {imagenumber, machine} = parm;
                txt = "";
                While( this << canread,
                    t =  (this << read);
                    txt = txt || t;
                    if(endswith(txt,"\!n")|endswith(txt,"\!r"),break());
                    wait(.01);
                );
                /////////////////////
                if(contains(txt,"rror"), write(txt));
                if(contains(txt,"okpic") & imgNumToIsStoredInGraph[imagenumber] != generation,
                    // grab the picture from OKPIC
                    local({left,right,top,bot,filename,rc,imgseg,bindata,matdata,img},
                        //show(imagenumber);
                        left = Mod( (imagenumber - 1), squaresize );
                        col = left * subImageWidth;
                        right = left + 1;
                        top = Floor( (imagenumber - 1) / squaresize );
                        row = top * subImageHeight;
                        bot = top + 1;
                        left = left / squaresize;
                        right = right / squaresize;
                        top = top / squaresize;
                        bot = bot / squaresize;
                        filename=Eval Insert( pythonCodeFilePathWin || "mbrot_^generation^_^Char( imagenumber )^.pbd" );
                        try(//
                            //show(filename);
                            bindata = Load Text File( filename, blob() );
                            matdata = Blob To Matrix( bindata, "float", 8, "little", subImageWidth );
                            compositionMatrix[(row+1)::(row+subImageHeight), (col+1)::(col+subImageWidth)] = matdata[0,0];
                            img = generateImage();// <<<<<<<<<<<<<<< resultGB<<inval;//generateImage();
                            imgNumToIsStoredInGraph[imagenumber] = generation;
                            imagesReceived += 1;
                            if(imagesReceived==squaresize^2,
                            	show(edges);
                            	dtThumb<<addrows(1);
                            	dtThumb:thumbTime[nrows(dtThumb)] = TickSeconds()-StartTime;
                            	dtThumb:thumbPicture[nrows(dtThumb)] = img;
                            	dtThumb:thumbData[nrows(dtThumb)] = compositionMatrix;
                            	dtThumb:thumbIteration[nrows(dtThumb)] = edges["deep"];
                            	dtThumb:thumbDigits[nrows(dtThumb)] = edges["prec"];
                            	dtThumb:thumbCenterX[nrows(dtThumb)] = edges["xcen"];
                            	dtThumb:thumbLeft[nrows(dtThumb)] = edges["xlef"];
                            	dtThumb:thumbRadiusX[nrows(dtThumb)] = edges["xrad"];
                            	dtThumb:thumbRight[nrows(dtThumb)] = edges["xrig"];
                            	dtThumb:thumbBottom[nrows(dtThumb)] = edges["ybot"];
                            	dtThumb:thumbCenterY[nrows(dtThumb)] = edges["ycen"];
                            	dtThumb:thumbRadiusY[nrows(dtThumb)] = edges["yrad"];
                            	dtThumb:thumbTop[nrows(dtThumb)] = edges["ytop"];
                            )
                        ,//
                            show(exception_msg);//
                        );//
                    );
                );
                
                // check all lines for xlef/xrig/ytop/ybot reports and stash them away
                gotedge=0;
                patmatch(txt,("xlef"|"xrig"|"ytop"|"ybot"|"prec"|"deep"|"xcen"|"xrad"|"ycen"|"yrad")>>edgekey + patspan(" ") +
                    patbreak("\!n")>>edgevalue + pattest(edges[edgekey]=edgevalue;gotedge=1;1) + "\!n" + patfail(/*get all of them at once*/));
                if(gotedge,
                    titleabove<<settext(edges["ytop"]);
                    titlebelow<<settext(edges["ybot"]);
                    titleleft<<settext(edges["xlef"]);
                    titleright<<settext(edges["xrig"]);
                );
                txttrim = regex(txt,"([^\!n\!r]+)[\!n\!r]+$","\1"); // keep just the most recent (last) line
                //show("trim",chartoblob(txt),ismissing(txttrim),txttrim);
                if(!ismissing(txttrim) & length(txttrim) > 5 & try(!contains(imgNumToMonTextBox[imagenumber]<<gettext,"bye"),1),//
                    Try(// 
                        if(contains(txttrim,"ERROR"),
                            if(try(imgNumToMachine[imagenumber]<<isreadeof,1),
                                imgNumToMachine[imagenumber] = 0;
                                show(txttrim,imagenumber,machine,"relaunching");
                                launch(imagenumber, machine);	
                            );
                            txttrim="aaaaa"||char(machine);
                        );
                        txttrim = trim(Substr( txttrim, 6 ));//
                        if(contains(txttrim,"%"),
                            try(imgNumToMonTextIfBox[imagenumber]<<set(0));
                            siz=round(progressbarlength*num(txttrim));
                            try((imgNumToMonProgressBox[imagenumber]<<child)<<width(siz));
                            try(((imgNumToMonProgressBox[imagenumber]<<child)<<sib)<<width(progressbarlength-siz));
                            imgNumToMonProgressIfBox[imagenumber]<<set(1);
                        ,
                            if(txttrim != "." & txttrim != "go go go",//speedup?
                                //show(txttrim);
                                try(imgNumToMonTextBox[imagenumber] << settext( left(txttrim,10) ));//
                            );
                            try(imgNumToMonTextIfBox[imagenumber]<<set(1));
                            try(imgNumToMonProgressIfBox[imagenumber]<<set(0));
                        );
                        imgNumToResponse[imagenumber] = txttrim;
                       // imgNumToMonTextBox[imagenumber] << updatewindow;//
                    , //
                        Show( exception_msg,"bad", chartoblob(txt) );//
                    );//
                );//
////////////////////////////////////               
            );
        );
    );
    //Show( imagenumber );
    imgNumToRP[imagenumber] << Write( "read prog\!n" );
    return( machinename[machine] );
);

imgNumToRP = [=> ];  // holder for active RunProgram objects
imgNumToRPdone = [=> ]; // holder for RunProgram objects that need a moment to terminate
imgNumToMachine = [=> 0]; // which machine is this running on?
imgNumToResponse = [=> ""]; // f(x) labels for graphs
imgNumToMonTextBox = [=> ];
imgNumToMonTextIfBox = [=> ];
imgNumToMonProgressBox = [=> ];
imgNumToMonProgressIfBox = [=> ];
imgNumToIsStoredInGraph = [=>-1];
compositionMatrix = j(ydimension,xdimension,1);//imgNumToSegInGraph = [=>];

stopImgNum = Function( {imgNum},
    imgNumToRP[imgNum] << Write( "quit\!n" )
);
stopall = Function( {},
    {k, rp, i},
    dropall();
    For Each( {{k, rp}, i}, imgNumToRP, stopImgNum( k ) );
    print("all stopping messages sent");
);
dropall = Function( {},
    {k, rp, i},
    For Each( {{k, rp}, i}, imgNumToRP, rp << Write( "drop\!n" ) );
);

startdrag=function({x,y},
    if(ismissing(dragpos0[1]),
    	dragpos0[1]=x;dragpos0[2]=y;//
    	dragpos1[1]=x;dragpos1[2]=y;//
    //	show("cap",dragpos0);
    ,//
    	dragpos1[1]=x;dragpos1[2]=y;//
     //   show("-->",dragpos1);
    )
);

enddrag=function({x,y},
	dragpos1[1]=x;dragpos1[2]=y;//
	dropall(); // stop old work

/*
redp ^redpower^
grnp ^greenpower^
blup ^bluepower^
*/	
    commands =
"prec ^digitPrecision^
deep ^drillDepth^
xzfr ^(dragpos0[1]+dragpos1[1])/2^
yzfr ^(dragpos0[2]+dragpos1[2])/2^
zfra ^max(abs(dragpos0[1]-dragpos1[1]),abs(dragpos0[2]-dragpos1[2]))^
zoom
name ^pythonCodeFilePathLinux^mbrot_^generation^_^Char( k )^.pbd
make
";
imgNumToIsStoredInGraph = [=>-1];
//foreach({{k,v},i},imgNumToSegInGraph,v<<remove;);
compositionMatrix[0,0] = 1; //imgNumToSegInGraph = [=>];
imagesReceived = 0;
sendCommandsToAll_k_rp_i(commands);

dragpos0[1]=.;dragpos0[2]=.;
dragpos1[1]=.;dragpos1[2]=.;

	StartTime=tickseconds();
);

imagesReceived = 0;
StartTime=tickseconds();

dragpos0=[. .];
dragpos1=[.,.];
resultWindow = New Window( "result",
    V List Box(
        H Center Box( titleabove = Text Box( "above" ) ),
        H List Box(
            V Center Box( titleleft = Text Box( "left", <<Rotate Text( "left" ) ) ),
            resultGB = Graph Box(//
                framesize( xdimension * resultWindowScale, ydimension * resultWindowScale ),
                X Scale( 0, 1 ),
                suppressaxes,
                <<yaxis( Min( 1 ), Max( 0 ) ),
                <<AddGraphicsScript("front",
                	mousetrap(startdrag(x,y),enddrag(x,y));
                	
                	//nope. not here. generateImage(); // can this be here? it might be modifying the display tree during show...
                	
                	if( !ismissing(dragpos0[1]),
                        pencolor("red");
                        linestyle("dotted");
                        line({dragpos0[1],dragpos0[2]},{dragpos0[1],y});
                        line({dragpos0[1],y},{x,y});
                        line({x,y},{x,dragpos0[2]});
                        line({x,dragpos0[2]},{dragpos0[1],dragpos0[2]});
                        //rect(min(dragpos0[1],x),max(dragpos0[2],y),max(dragpos0[1],x),min(dragpos0[2],y),1)
                	);
                )
            ),
            V Center Box( titleright = Text Box( "right", <<Rotate Text( "left" ) ) )
        ),
        H Center Box( titlebelow = Text Box( "below" ) ),
        hlistbox(
        	textbox("color"),sliderbox(.1,10,colorpower,generateimage()),
        	spacerbox(size(20,1)),
        	globalbox(drillDepth),
        	spacerbox(size(20,1)),
        	globalbox(digitPrecision)
        	
        )
    )
);

monitorWindow = New Window( "monitor",
    V List Box(//
        H List Box(//
            Button Box( "stop all", stopall() ), //
           // Button Box( "drop all", dropall() ) //
        ), //
        lub = Lineup Box( N Col( squaresize ), 
            spacing( 3, 3 ),
        	<<backgroundcolor(rgbcolor(.94,.97,.97))
        ) //
    ) //
);

write("\!ncheckpoint 1\!n");wait(.1);

monbox = Function( {imgNum, macnam},
    Eval(
        Eval Expr(
            H List Box(
                Spacer Box( size( 1, 30 ), color( rgbcolor(1,1,1) ) ),// left edge
                Spacer Box( size( 5, 1 ), color( rgbcolor(1,1,1) ) ),// first part of top edge
                V List Box(
                    Spacer Box( size( progressbarlength, 1 ), color( rgbcolor(1,1,1) ) ), // most of top edge
                    spacerbox(size(1,5)),
                    imgNumToMonTextIfBox[imgNum] = If Box( 1, imgNumToMonTextBox[imgNum] = Text Box( Char( imgNum ) ) ),
                    imgNumToMonProgressIfBox[imgNum] = If Box(
                        0,
                        imgNumToMonProgressBox[imgNum] = H List Box(
                            Spacer Box( size( 0, 7 ), color( "green" ) ),
                            Spacer Box( size( progressbarlength, 7 ), color( rgbcolor(.9,.9,.9) ) )
                        )
                    ),
                    Text Box( macnam )
                ),
                <<backgroundcolor(rgbcolor(.95,.95,.95))
            )
        )
    )
);

write("\!ncheckpoint 2\!n");wait(.1);

wait(.1);monitorWindow<<movewindow(0,0);wait(.1);

write("\!ncheckpoint 3\!n");wait(.1);

For( i = 1, i <= squaresize^2, i += 1,
    lub << append( monbox( i, launch( i ) ) );
   //wait(0);
);

write("\!ncheckpoint 4\!n");wait(.1);

generation = 0;
sendCommandsToAll_k_rp_i = function({commands},
    if(contains(commands,"make"), generation+=1;show(generation));
    imgNumToResponse = [=> ""];
    For Each( {{k, rp},i}, imgNumToRP,
        rp << Write( regex(evalinsert(commands),"[\!r\!n]+","\!n",GLOBALREPLACE ));// change returns to newlines for linux
    );
    show("broadcast commands done");
);

// bring it up. the threads need some coaxing.
commands = 
"zfra ^z^
read zfra
";

write("\!ncheckpoint 5\!n");wait(.1);

For( z = 1, z < 10, z += 1,
    sendCommandsToAll_k_rp_i( Eval Insert( commands ) );
    Print( "Ready", z );
    start = tickseconds();
    While( N Items( imgNumToResponse ) < N Items( imgNumToRP ),
        if(tickseconds()-start > 10 & N Items( imgNumToResponse ) == 0, throw("nothing is happening"));
        write(evalinsert("\!nwaiting for responces ^N Items( imgNumToResponse )^, ^N Items( imgNumToRP )^")); 
        Wait( 1 );
    );
    needmore = 0;
    For Each( {{k, v}}, imgNumToResponse, If( trim(v) != char(z), needmore = 1 ) );
    if(needmore==0,break());
);

write("\!ncheckpoint 6\!n");wait(.1);

begin=tickseconds();
// make same image... a 10x10 processor grid on the beowulf cluster seems good.
// 1080      /2/2/2/3/3/3/5
// 1920  /2/2/2/2/2/2/2/3/5
// 2 3 4 5 6 8 10 12 15 20 24 30 40
/*
x  -1.94062988168240384898793330252067703470368414336
x  -1.94062988168240384898793330252067703470368414330
y   0.00100077622840272515587559395423166491124004188
y   0.00100077622840272515587559395423166491124004184
*/

/*
redp ^redpower^
grnp ^greenpower^
blup ^bluepower^
xcen -0.5396920368085179367059520359131710652013065 
ycen -0.6611804024210354957147719385629364551294705 

ycen  -0.049842152966738802163856592043818211312978962993082255019266085623 
      -0.049842152966738802163856592043818211312978962993082255019266085622
xcen   0.336447185016324548337106220487657679616547644208117764858297769859
       0.336447185016324548337106220487657679616547644208117764858297769858
*/
commands =
"read prog
prec ^digitPrecision^
deep ^drillDepth^
xcen -1.9406298816824038489879333025206770347036841433
xrad 3
ycen 0.0010007762284027251558755939542316649112400418
yrad 2
xpix ^xdimension^
ypix ^ydimension^
ncpu ^squaresize*squaresize^
icpu ^i-1^
name ^pythonCodeFilePathLinux^mbrot_^generation^_^Char( k )^.pbd
make
";
sendCommandsToAll_k_rp_i(commands);

print("polling...");// check for okpic and *BAD*
keepGoing = 1;
While( keepGoing,
    keepGoing = 0;
    For Each( {{k, rp},i}, imgNumToRP,
        if( (rp<<CanRead) | (rp<<CanWrite()) | !(rp<<IsReadEOF), 
            if( trim(imgNumToResponse[k]) != "okpic" & trim(imgNumToResponse[k]) != "*BAD*",
                keepGoing = 1;//
                commands = "read prog\!n";
                rp << Write( commands );
                wait(1);//
            );
        );
    );
);
show(tickseconds()-begin,edges);
print("JSL stopped");

 

Last Modified: Jan 1, 2022 6:59 PM