text:
#!/usr/bin/env julia

# HOW TO RUN THIS SCRIPT:
# 1: check that you have a recent(ish, >1.0) version of Julia installed
# 2: within Julia, make sure you have the required packages (below, in the "Using..." statement) installed (I think they all come by default though)
# 3: make sure you have FFmpeg installed
# 4: setup an empty sub-directory called "tmpPlots"
# 5: go to line 693 and uncomment the `makeAnim()` function call
# 6: save this file
# 7: run this file (double click it and tell it to run with Julia, open a terminal and type `julia threeBodyProb.jl`, or start a julia session and type `include("threeBodyProb.jl")` )
# 8: profit

# following the above steps should generate a random three-body problem video in 720p without music
# to add music, specify initial conditions, or otherwise tinker...read the code/README and the comments in the code!
# this has only been tested on Linux (Mint) and Windows 10, so you might run into small additional issues on macOS?

using Plots, Random, Printf, Plots.Measures, Dates, CSV, DataFrames, DataFramesMeta

function initCondGen() #get random initial conditions for mass/radius, position, and velocity
    function getMass(nBodies) #generate random masses that better reflect actual stellar populations, although not currently using because it's boring
        mList=zeros(nBodies)
        N=(0.5^(-1.3)-150^(-1.3))/1.3 #crude approximation of IMF integral assuming alpha = 2.3, stellar mass range of 0.5:150 solar masses
        rescale=1e6
        max=floor(Int,N*rescale)
        for i=1:nBodies
            intTarget=rand(0:max,1)[1]/rescale
            m=(0.5^(-1.3)-intTarget*1.3)^(-1/1.3) #just algebra from above
            mList[i]=round(m,digits=2)
        end
        return mList
    end
    m=rand(1:1500,3)./10 #3 random masses between 0.1 and 150 solar masses, uniform distribution
    #m=getMass(3) #get mass from IMF -- this way is kind of boring...so not using it, but left here in case I change my mind?
    rad=m.^0.8 #3 radii based on masses in solar units
    m=m.*2e30 #convert to SI kg
    rad=rad.*7e8 #convert to SI m
    pos1=rand(-25:25,2) #random initial coordinates x & y for first body, AU
    function genPos2(pos1)
        accept2=false
        while accept2==false
            pos2=rand(-25:25,2) #random initial coordinates for second body, AU
            dist21=sqrt((pos1[1]-pos2[1])^2+(pos1[2]-pos2[2])^2)
            if (dist21*1.5e11)>(rad[1]+rad[2]) #they aren't touching
                accept2=true
                return pos2
            end
        end
    end
    pos2=genPos2(pos1)
    function genPos3(pos1,pos2)
        accept3=false
        while accept3==false
            pos3=rand(-25:25,2) #random initial coordinates for third body, AU
            dist31=sqrt((pos1[1]-pos3[1])^2+(pos1[2]-pos3[2])^2)
            dist32=sqrt((pos2[1]-pos3[1])^2+(pos2[2]-pos3[2])^2)
            if (dist31*1.5e11)>(rad[1]+rad[3]) && (dist32*1.5e11)>(rad[2]+rad[3]) #3rd isn't touching either
                accept3=true
                return pos3
            end
        end
    end
    pos3=genPos3(pos1,pos2)
    pos=[pos1[1],pos1[2],pos2[1],pos2[2],pos3[1],pos3[2]].*1.5e11 #convert accepted positions to SI, m
    v=rand(-7e3:7e3,6) #random x & y velocities with mag between -10 & 10 km/s, totally arbitrary...
    #r=[x1,y1,x2,y2,x3,y3,v1x,v1y,v2x,v2y,v3x,v3y] -- format if you want to specify your own initial conditions
    r=[pos[1],pos[2],pos[3],pos[4],pos[5],pos[6],v[1],v[2],v[3],v[4],v[5],v[6]]
    open("initCond.txt","w") do f #save initial conditions to file in folder where script is run
        write(f,"m1=$(@sprintf("%.1f",(m[1]/2e30))) m2=$(@sprintf("%.1f",(m[2]/2e30))) m3=$(@sprintf("%.1f",(m[3]/2e30))) (solar masses)\nv1x=$(v[1]/1e3) v1y=$(v[2]/1e3) v2x=$(v[3]/1e3) v2y=$(v[4]/1e3) v3x=$(v[5]/1e3) v3y=$(v[6]/1e3) (km/s)\nx1=$(pos1[1]) y1=$(pos1[2]) x2=$(pos2[1]) y2=$(pos2[2]) x3=$(pos3[1]) y3=$(pos3[2]) (AU from center)")
    end
    return r, rad, m
end

function dR(r,m;energyBool=0) #function we will use RK4 on to approximate solution
    G=6.67408313131313e-11# Nm^2/kg^2
    M1,M2,M3=m[1],m[2],m[3] #kg
    x1,x2,x3=r[1],r[3],r[5] #m
    y1,y2,y3=r[2],r[4],r[6] #m

    c1,c2,c3=G*M1,G*M2,G*M3 #Nm^2/kg
    r1_2=sqrt((x1-x2)^2+(y1-y2)^2) #distance from 1->2, m
    r1_3=sqrt((x1-x3)^2+(y1-y3)^2) #distance from 1->3, m
    r2_3=sqrt((x2-x3)^2+(y2-y3)^2) #distance from 2->3, m

    v1X,v2X,v3X=r[7],r[9],r[11] #these are our change in position after dt (dr/dt*dt=dr)
    v1Y,v2Y,v3Y=r[8],r[10],r[12] #m after * dt

    #get change in velocity from accelerations (d^2r/dt^2*dt=dv/dt*dt=dv)
    aX1=-(c2*(x1-x2)/(r1_2^3))-(c3*(x1-x3)/(r1_3^3)) #d^2x/dt^2 for 1, m/s after * dt
    aX2=-(c1*(x2-x1)/(r1_2^3))-(c3*(x2-x3)/(r2_3^3)) #d^2x/dt^2 for 2, m/s
    aX3=-(c1*(x3-x1)/(r1_3^3))-(c2*(x3-x2)/(r2_3^3)) #d^2x/dt^2 for 3, m/s
    aY1=-(c2*(y1-y2)/(r1_2^3))-(c3*(y1-y3)/(r1_3^3)) #d^2y/dt^2 for 1, m/s
    aY2=-(c1*(y2-y1)/(r1_2^3))-(c3*(y2-y3)/(r2_3^3)) #d^2y/dt^2 for 2, m/s
    aY3=-(c1*(y3-y1)/(r1_3^3))-(c2*(y3-y2)/(r2_3^3)) #d^2y/dt^2 for 3, m/s

    global energy #keep track of energy loss from RK4 error
    U=-G*M1*M2/r1_2-G*M2*M3/r2_3-G*M1*M3/r1_3 #grav potential
    K=0.5*M1*(v1X^2+v1Y^2)+0.5*M2*(v2X^2+v2Y^2)+0.5*M3*(v3X^2+v3Y^2) #kinetic
    if energyBool==1
        push!(energy,K+U) #total system energy
    end

    return [v1X,v1Y,v2X,v2Y,v3X,v3Y,aX1,aY1,aX2,aY2,aX3,aY3]
end

function gen3Body(stopCond=[10,100],numSteps=10000,initCond=nothing) #default stop conditions of 10 yrs and 100 AU sep
    tStop=stopCond[1]*365*24*3600 #convert to SI s
    sepStop=stopCond[2]*1.5e11 #convert to SI m
    stop=false
    currentT=0
    t=range(0,stop=tStop,length=(numSteps+1)) #+1 because I don't want 0 to count
    stepSize=tStop/numSteps
    x1=zeros(length(t))
    y1=zeros(length(t))
    x2=zeros(length(t))
    y2=zeros(length(t))
    x3=zeros(length(t))
    y3=zeros(length(t))
    if(initCond !== nothing)
        r,m,rad=initCond

        # println("r: $r")
        # println("n: $m")
        # println("rad: $rad")

        # open("initCond.txt","w") do f #save initial conditions to file in folder where script is run
        #     write(f,"preset: m1=$(@sprintf("%.1f",(m[1]/2e30))) m2=$(@sprintf("%.1f",(m[2]/2e30))) m3=$(@sprintf("%.1f",(m[3]/2e30))) (solar masses)\nv1x=$(r[7]/1e3) v1y=$(r[8]/1e3) v2x=$(r[9]/1e3) v2y=$(r[10]/1e3) v3x=$(r[11]/1e3) v3y=$(r[12]/1e3) (km/s)\nx1=$(r[1]/1.5e11) y1=$(r[2]/1.5e11) x2=$(r[3]/1.5e11) y2=$(r[4]/1.5e11) x3=$(r[5]/1.5e11) y3=$(r[6]/1.5e11) (AU from center)")
        # end
    else
        r,rad,m=initCondGen()

        # println("r: $r")
        # println("n: $m")
        # println("rad: $rad")
    end
    initV=copy(r[7:end])
    min12=rad[1]+rad[2]
    min13=rad[1]+rad[3]
    min23=rad[2]+rad[3]
    i=1
    stopT=maximum(t)
    collisionBool=false; collisionInds=[0,0]
    #implement RK4 to model solutions to differential equations
    while stop==false
        if currentT==stopT || currentT>stopT #in case of rounding error or something
            stop=true
        elseif i>(numSteps+1) #inf loop failsafe
            stop=true
            println("error: shouldn't have gotten here")
        else
            x1[i]=r[1] #store current positions
            y1[i]=r[2]
            x2[i]=r[3]
            y2[i]=r[4]
            x3[i]=r[5]
            y3[i]=r[6]

            k1=stepSize*dR(r,m,energyBool=1)
            k2=stepSize*dR(r.+0.5.*k1,m)
            k3=stepSize*dR(r.+0.5.*k2,m)
            k4=stepSize*dR(r.+k3,m)
            r+=(k1.+2.0*k2.+2.0.*k3.+k4)./6 #RK4 update to positions, velocities

            #check separation after each dt step
            sep12=sqrt((x1[i]-x2[i])^2+(y1[i]-y2[i])^2)
            sep13=sqrt((x1[i]-x3[i])^2+(y1[i]-y3[i])^2)
            sep23=sqrt((x3[i]-x2[i])^2+(y3[i]-y2[i])^2)

            if sep12<min12 || sep13<min13 || sep23<min23 || sep12>sepStop || sep13>sepStop || sep23>sepStop
                if sep12<min12 || sep13<min13 || sep23<min23
                    collisionBool=true
                    if sep12<min12
                        collisionInds=[1,2]
                    elseif sep13<min13
                        collisionInds=[1,3]
                    elseif sep23<min23
                        collisionInds=[2,3]
                    end
                else
                    collisionBool=false
                end
                stop=true #stop if collision happens or body is ejected
                t=range(0,stop=currentT,length=i) #t should match pos vectors
                x1=x1[1:i] #don't want trailing zeros
                y1=y1[1:i]
                x2=x2[1:i]
                y2=y2[1:i]
                x3=x3[1:i]
                y3=y3[1:i]
            end
            i+=1
            currentT+=stepSize #next step
        end
    end
    return [x1,y1,x2,y2,x3,y3], t, m, rad, collisionBool, collisionInds, initV
end

function getInteresting3Body(minTime=0, initCond=nothing) #in years, defaults to 0
    #sometimes (most of the time) random conditions result in a really short animation where things
    #just crash into each other/fly away, so this function throws away those
    yearSec=365*24*3600
    interesting=false
    i=1
    while interesting==false
        global energy=[] #re-initialize empty energy array
        plotData,t,m,rad,collisionBool,collisionInds,initV=gen3Body([60,200],600000, initCond)
        if initCond !== nothing  #if we are setting the valuas manually, don't log
            println(maximum(t)/yearSec) #tell me how many years we are simulating
            return plotData,t,m,rad,collisionBool,collisionInds
            interesting=true
        elseif (maximum(t)/yearSec)>minTime #only return if simulation runs for longer than minTime
            println(maximum(t)/yearSec) #tell me how many years we are simulating
            open("cron_log.txt","a") do f #for cron logging, a flag = append
                write(f,"$(maximum(t)/yearSec)\n")
            end
            open("3BodyStats.txt","a") do f #for stats logging
                initPos=[plotData[1][1],plotData[2][1],plotData[3][1],plotData[4][1],plotData[5][1],plotData[6][1]]./1.5e11 #AU
                write(f,"$(today()),$(maximum(t)/yearSec),$(m[1]/2e30),$(m[2]/2e30),$(m[3]/2e30),$(rad[1]/7e8),$(rad[2]/7e8),$(rad[3]/7e8),$collisionBool,$(collisionInds[1]),$(collisionInds[2]),$(initPos[1]),$(initPos[2]),$(initPos[3]),$(initPos[4]),$(initPos[5]),$(initPos[6]),$(initV[1]/1e3),$(initV[2]/1e3),$(initV[3]/1e3),$(initV[4]/1e3),$(initV[5]/1e3),$(initV[6]/1e3),$i\n")
            end
            return plotData,t,m,rad,collisionBool,collisionInds
            interesting=true
        elseif i>1999 #computationally expensive so don't want to go forever
            interesting=true #render it anyways I guess because sometimes it's fun?
            println("did not find interesting solution in number of tries allotted, running anyways")
            println(maximum(t)/yearSec) #how many years simulation runs for
            open("cron_log.txt","a") do f #for cron logging
                write(f,"found a solution with t = $(maximum(t)/yearSec) in $i iterations\n")
            end
            open("3BodyStats.txt","a") do f #log stats, see "3BodyAnalysis.ipynb"
                initPos=[plotData[1][1],plotData[2][1],plotData[3][1],plotData[4][1],plotData[5][1],plotData[6][1]]./1.5e11 #AU
                write(f,"$(today()),$(maximum(t)/yearSec),$(m[1]/2e30),$(m[2]/2e30),$(m[3]/2e30),$(rad[1]/7e8),$(rad[2]/7e8),$(rad[3]/7e8),$collisionBool,$(collisionInds[1]),$(collisionInds[2]),$(initPos[1]),$(initPos[2]),$(initPos[3]),$(initPos[4]),$(initPos[5]),$(initPos[6]),$(initV[1]/1e3),$(initV[2]/1e3),$(initV[3]/1e3),$(initV[4]/1e3),$(initV[5]/1e3),$(initV[6]/1e3),$i\n")
            end
            return plotData,t,m,rad,collisionBool,collisionInds
        end
        i+=1
    end
end

function detectOrbiting(d1_2,d1_3,d2_3,m,x,y) #determines if 2 bodies are orbiting, so we should use their center of mass for frame calculation
    if d1_2/d2_3 > 2 && d1_3/d2_3 > 2 #objects 2 and 3 are orbiting?
        orbiting=23
        cmX=(m[2]*x[2]+m[3]*x[3])/(m[2]+m[3]) #get centers of mass to use in limit calculations to prevent oscillations
        cmY=(m[2]*y[2]+m[3]*y[3])/(m[2]+m[3])
        xNew=[x[1],cmX]
        yNew=[y[1],cmY]
        return orbiting,xNew,yNew
    elseif d2_3/d1_2 > 2 && d1_3/d1_2 > 2 #objects 2 and 1 are orbiting?
        orbiting=21
        cmX=(m[2]*x[2]+m[1]*x[1])/(m[2]+m[1]) #get centers of mass
        cmY=(m[2]*y[2]+m[1]*y[1])/(m[2]+m[1])
        xNew=[x[3],cmX]
        yNew=[y[3],cmY]
        return orbiting,xNew,yNew
    elseif d1_2/d1_3 > 2 && d2_3/d1_3 > 2 #objects 1 and 3 are orbiting?
        orbiting=13
        cmX=(m[1]*x[1]+m[3]*x[3])/(m[1]+m[3]) #get centers of mass
        cmY=(m[1]*y[1]+m[3]*y[3])/(m[1]+m[3])
        xNew=[x[2],cmX]
        yNew=[y[2],cmY]
        return orbiting,xNew,yNew
    else #no pairs orbiting
        return 0,x,y
    end
end

function getLims(xNew,yNew,padding,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD) #computes limits given padding and offsets
    cX,cY=sum(xNew)/length(xNew),sum(yNew)/length(yNew) #next thing to change..?
    dx=maximum(xNew)-minimum(xNew); dy=maximum(yNew)-minimum(yNew)
    dF = dx<dy ? dy : dx #do we use dx or dy for the frame?
    xlims=[(cX+ΔCx)-padding-dF/2+ΔL,(cX+ΔCx)+padding+dF/2+ΔR]
    ylims=[(cY+ΔCy)-padding-dF/2+ΔD,(cY+ΔCy)+padding+dF/2+ΔU]
    return xlims,ylims
end

function getΔC(target,start,pos,extraDx,extraDy,x,y,padding,tol=0.0001,maxIter=100000) #find center shifts brute force
    targCx,targCy,targxlims,targylims = target #these are the "old" limits we want offset to
    cx,cy = start #from new limits
    ΔCx,ΔCy = cx-targCx,cy-targCy #initial "guess"
    diffxList = [0.,0.,0.]; diffyList = [0.,0.,0.]
    xtargList = [0.,0.,0.]; ytargList = [0.,0.,0.]
    ΔL,ΔR = extraDx; ΔU,ΔD = extraDy
    xlims,ylims = getLims(x,y,padding,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD)
    diff(r,rTarg)=abs(r-rTarg)
    for i = 1:length(pos)
        Bx,By = pos[i]
        rx,ry = relative(xlims,ylims,Bx,By) #where is it in terms of the frame width?
        rxTarg,ryTarg = relative(targxlims,targylims,Bx,By) #where was it? (where it should be)
        xtargList[i] = rxTarg; ytargList[i] = ryTarg
        diffxList[i] = diff(rx,rxTarg); diffyList[i] = diff(ry,ryTarg)
    end
    diffx,xInd = findmax(diffxList); diffy,yInd = findmax(diffyList) #find the one that shifted the most
    Bx = pos[xInd][1]; By = pos[yInd][2]
    rxTarg = xtargList[xInd]; ryTarg = ytargList[yInd]
    if diffx<tol && diffy<tol #not that big a difference, just return the initial guess
        return xlims,ylims,ΔCx,ΔCy
    else
        dx = targxlims[2]-targxlims[1]; dy = targylims[2]-targylims[1]
        function getDir(targxlims,targylims,ΔCx,ΔCy,tol,dx,dy) #figure out which way we need to shift
            acceptXDir = false; acceptYDir = false
            signdx = diffx<tol ? 0 : -1; signdy = diffy<tol ? 0 : -1
            counter = 1
            while acceptXDir == false || acceptYDir == false
                if acceptXDir == false
                    signdx=signdx^counter #flip the direction if what we tried before didn't work
                end
                if acceptYDir == false
                    signdy=signdy^counter
                end
                guessX = ΔCx+signdx*dx*tol/10; guessY = ΔCy+signdy*dy*tol/10
                xlims,ylims = getLims(x,y,padding,guessX,guessY,ΔL,ΔR,ΔU,ΔD)
                rx,ry = relative(xlims,ylims,Bx,By)
                newDiffx=diff(rx,rxTarg); newDiffy=diff(ry,ryTarg)
                acceptXDir = newDiffx<=diffx; acceptYDir = newDiffy<=diffy #are we going in the right direction?
                if counter>2
                    println("PROBLEM: changing center sign has no effect")
                    break
                end
                counter+=1
            end
            return signdx,signdy,guessX,guessY
        end
        signdx,signdy,guessX,guessY = getDir(targxlims,targylims,ΔCx,ΔCy,tol,dx,dy)
        counter = 2
        stopX = false; stopY = false
        while diffx>tol || diffy>tol #keep guessing in the right direction until we're within the tolerance
            #this could actaully probably just be calculated? whatever "if it ain't broke don't fix it" and
            #this is not high performance code
            if stopX == false
                guessX = ΔCx+signdx*dx*tol/10*counter
            end
            if stopY == false
                guessY = ΔCy+signdy*dy*tol/10*counter
            end
            xlims,ylims = getLims(x,y,padding,guessX,guessY,ΔL,ΔR,ΔU,ΔD)
            rx,ry = relative(xlims,ylims,Bx,By)
            diffx=diff(rx,rxTarg); diffy=diff(ry,ryTarg)
            stopX = diffx<tol; stopY = diffy<tol #compare new guesses and see if we should stop guessing
            counter+=1
            if counter == maxIter
                println("PROBLEM: did not converge in $maxIter iterations")
                break
            end
        end
        return xlims,ylims,guessX,guessY
    end
end

relative(xlims,ylims,x,y)=(x-xlims[1])/(xlims[2]-xlims[1]),(y-ylims[1])/(ylims[2]-ylims[1]) #returns position in terms of frame widths
center(xy) = sum(xy)/length(xy)

function comparePos(orbitOld,orbiting,m,x,y,padding,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD) #try to prevent jumps when switching modes
    orbitStr = orbiting!=0 ? string(orbiting) : string(orbitOld) #was the old one orbiting or is this one orbiting?
    i1,i2=parse(Int64,string(orbitStr[1])),parse(Int64,string(orbitStr[2])) #indices of two orbiting bodies
    inds = [1,2,3]; otherInd = 0 #trying to make generalization to n-bodies easier
    for i = 1:length(inds)
        if inds[i] != i1 && inds[i] != i2
            otherInd = i
        end
    end
    cmX=(m[i1]*x[i1]+m[i2]*x[i2])/(m[i1]+m[i2]) #get centers of mass
    cmY=(m[i1]*y[i1]+m[i2]*y[i2])/(m[i1]+m[i2])
    xNew=[x[otherInd],cmX] #new is a bit of a misnomer because sometimes it's old
    yNew=[y[otherInd],cmY] #point is we use this one for orbiting, the other for not
    xlimsOrbit,ylimsOrbit=getLims(xNew,yNew,padding,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD)
    xlimsNorm,ylimsNorm=getLims(x,y,padding,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD)
    function getOld(orbitOld,xlimsOrbit,ylimsOrbit,xlimsNorm,ylimsNorm,x,y,xNew,yNew) #which one was the old one?
        if orbitOld != 0 #thing was orbiting
            oldxlims,oldylims = xlimsOrbit,ylimsOrbit
            oldCx,oldCy = center(xNew),center(yNew)
            return oldxlims,oldylims,oldCx,oldCy
        else
            oldxlims,oldylims = xlimsNorm,ylimsNorm
            oldCx,oldCy = center(x),center(y)
            return oldxlims,oldylims,oldCx,oldCy
        end
    end
    oldxlims,oldylims,oldCx,oldCy = getOld(orbitOld,xlimsOrbit,ylimsOrbit,xlimsNorm,ylimsNorm,x,y,xNew,yNew)
    if orbiting != 0 #transitioning to orbiting, frame instantaneously wants to shrink
        cx = center(xNew); cy = center(yNew)
        ΔL = oldxlims[1]-xlimsOrbit[1]; ΔR = oldxlims[2]-xlimsOrbit[2]
        ΔU = oldylims[2]-ylimsOrbit[2]; ΔD = oldylims[1]-ylimsOrbit[1]
        extraDx = [ΔL,ΔR]; extraDy = [ΔU,ΔD]
        xlims,ylims,ΔCx,ΔCy = getΔC([oldCx,oldCy,oldxlims,oldylims],[cx,cy],[[x[1],y[1]],[x[2],y[2]],[x[3],y[3]]],extraDx,extraDy,xNew,yNew,padding)
        return xlims,ylims,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD
    else #transitioning from orbiting, frame instantaneously wants to expand
        cx = center(x); cy = center(y)
        ΔL = oldxlims[1]-xlimsNorm[1]; ΔR = oldxlims[2]-xlimsNorm[2]
        ΔU = oldylims[2]-ylimsNorm[2]; ΔD = oldylims[1]-ylimsNorm[1]
        extraDx = [ΔL,ΔR]; extraDy = [ΔU,ΔD]
        xlims,ylims,ΔCx,ΔCy = getΔC([oldCx,oldCy,oldxlims,oldylims],[cx,cy],[[x[1],y[1]],[x[2],y[2]],[x[3],y[3]]],extraDx,extraDy,x,y,padding)
        return xlims,ylims,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD
    end
end

function computeLimits(pos,posFuture,padding,m,orbitOld,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD) #determines plot limits at each frame, padding in units of pos
    x=[pos[1],pos[3],pos[5]]
    y=[pos[2],pos[4],pos[6]]
    d1_2=sqrt((x[1]-x[2])^2 + (y[1]-y[2])^2)
    d1_3=sqrt((x[1]-x[3])^2 + (y[1]-y[3])^2)
    d2_3=sqrt((x[2]-x[3])^2 + (y[2]-y[3])^2)
    orbiting,xNew,yNew = detectOrbiting(d1_2,d1_3,d2_3,m,x,y) #are they orbiting?
    if orbiting != orbitOld #are we switching modes?
        xlims,ylims,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD = comparePos(orbitOld,orbiting,m,x,y,padding,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD)
    else #"slowly" make the offsets go to zero, produces a nice smooth camera motion as we adjust
        relax = 0.95
        ΔCx*=relax;ΔCy*=relax;ΔL*=relax;ΔR*=relax;ΔU*=relax;ΔD*=relax
        xlims,ylims = getLims(xNew,yNew,padding,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD)
    end
    cNew = [(xlims[2]-xlims[1])/2+xlims[1],(ylims[2]-ylims[1])/2+ylims[1]] #dx+min(x); dy+min(y)
    return xlims,ylims,cNew,orbiting,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD
end

function getColors(m,c) #places colors of objects according to mass/size
    #c=[:biggest,:medium,:smallest] (order of input colors)
    maxM=maximum(m)
    minM=minimum(m)
    colors=[:blue,:blue,:blue] #testing
    if m[1]==maxM
        colors[1]=c[1]
        if m[2]==minM
            colors[2]=c[3]
            colors[3]=c[2]
        else
            colors[3]=c[3]
            colors[2]=c[2]
        end
    elseif m[2]==maxM
        colors[2]=c[1]
        if m[1]==minM
            colors[1]=c[3]
            colors[3]=c[2]
        else
            colors[3]=c[3]
            colors[1]=c[2]
        end
    else
        colors[3]=c[1]
        if m[1]==minM
            colors[1]=c[3]
            colors[2]=c[2]
        else
            colors[2]=c[3]
            colors[1]=c[2]
        end
    end
    return colors
end

function makeCircleVals(r,center=[0,0]) #makes circle values for the stars to plot
    xOffset=center[1]
    yOffset=center[2]
    xVals=[r*cos(i)+xOffset for i=0:(pi/64):(2*pi)]
    yVals=[r*sin(i)+yOffset for i=0:(pi/64):(2*pi)]
    return xVals,yVals
end

function saveFrames(plotData,t,m,rad,collisionBool,collisionInds)
    c=[:DodgerBlue,:Gold,:Tomato] #most massive to least massive, also roughly corresponds to temp
    colors=getColors(m,c)
    #adding fake stars
    numStars=2500
    starsX=zeros(numStars)
    starsY=zeros(numStars)
    for i=1:numStars
        num=rand(-200:200,2) #box size is 70 AU but we need some extra padding for movement
        starsX[i]=num[1]
        starsY[i]=num[2]
    end

    function getRatioRight(ratio,dx,dy) #makes sure the frame matches the ratio we want (for Twitter, square)
        if (dx/dy)!=ratio
            if dx>(ratio*dy)
                dy=dx/ratio
            else
                dx=dy*ratio
            end
        end
        return dx,dy
    end

    function relative(p::Plots.Subplot, rx, ry) #so I can plot in relative to parent
       xlims=Plots.xlims(p)
       ylims=Plots.ylims(p)
       return xlims[1]+rx*(xlims[2]-xlims[1]), ylims[1]+ry*(ylims[2]-ylims[1])
    end

    frameNum=1 #initialize frame counter
    stop=length(t) #what index the end is
    if collisionBool==true
        stop=length(t)-600 #-600 because we want to do the last 600 frames in slo-mo
    end
    #initialize a bunch of other things we'll need
    listInd=0
    limList=[]
    ratio=1
    offsetX = [0.,0.]; offsetY = [0.,0.]
    orbitOld = 0
    center = [0.,0.]; vel = [0.,0.]
    ΔCx = 0.;ΔCy = 0.;ΔL = 0.;ΔR = 0.;ΔU = 0.;ΔD = 0.
    println("energy loss = $((energy[end]-energy[1])/energy[1]*100) %") #anecdotally this is usually very small, but occasionally gets as high as ~1%
    for i=1:333:stop #this makes animation scale ~1 sec/year with other conditions
        GR.inline("png") #added to eneable cron/jobber compatibility, also this makes frames generate WAY faster? Prior to adding this when run from cron/jobber frames would stop generating at 408 for some reason.
        gr(legendfontcolor = plot_color(:white)) #legendfontcolor=:white plot arg broken right now (at least in this backend)
        print("$(@sprintf("%.2f",i/length(t)*100)) % complete\r") #output percent tracker
        pos=[plotData[1][i],plotData[2][i],plotData[3][i],plotData[4][i],plotData[5][i],plotData[6][i]] #current pos
        future = i+500<stop ? i+500 : i #make sure we don't go past end of data
        posFuture=[plotData[1][future],plotData[2][future],plotData[3][future],plotData[4][future],plotData[5][future],plotData[6][future]] #future pos
        limx,limy,center,orbitOld,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD=computeLimits(pos./1.5e11,posFuture./1.5e11,15,m,orbitOld,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD) #compute limits in AU, 15 AU padding
        dx,dy=(limx[2]-limx[1]),(limy[2]-limy[1])
        dx,dy=getRatioRight(ratio,dx,dy) #check ratio
        if listInd>1
            oldLimx,oldLimy=limList[listInd][1],limList[listInd][2]
            oldDx,oldDy=oldLimx[2]-oldLimx[1],oldLimy[2]-oldLimy[1]
            maxContraction=0.98; maxExpansion=1.02
            if dx/oldDx<maxContraction #frame shrunk more than x%
                limx[1]=center[1]-oldDx*maxContraction/2
                limx[2]=center[1]+oldDx*maxContraction/2
                limy[1]=center[2]-oldDx*maxContraction/2
                limy[2]=center[2]+oldDx*maxContraction/2
            elseif dx/oldDx>maxExpansion #grew more than x%
                limx[1]=center[1]-oldDx*maxExpansion/2
                limx[2]=center[1]+oldDx*maxExpansion/2
                limy[1]=center[2]-oldDx*maxExpansion/2
                limy[2]=center[2]+oldDx*maxExpansion/2
            elseif dy/oldDy<maxContraction #shrunk more than y%
                limx[1]=center[1]-oldDy*maxContraction/2
                limx[2]=center[1]+oldDy*maxContraction/2
                limy[1]=center[2]-oldDy*maxContraction/2
                limy[2]=center[2]+oldDy*maxContraction/2
            elseif dy/oldDy>maxExpansion #grew more than y%
                limx[1]=center[1]-oldDy*maxExpansion/2
                limx[2]=center[1]+oldDy*maxExpansion/2
                limy[1]=center[2]-oldDy*maxExpansion/2
                limy[2]=center[2]+oldDy*maxExpansion/2
            end
        end
        listInd+=1
        dx,dy=(limx[2]-limx[1]),(limy[2]-limy[1])
        dx,dy=getRatioRight(ratio,dx,dy) #check again
        limx = [center[1]-dx/2,center[1]+dx/2]; limy = [center[2]-dy/2,center[2]+dy/2]
        push!(limList,[limx,limy]) #record limits for later use, push! is bad and we should just preallocate this but whatever
        p=plot(plotData[1][1:33:i]./1.5e11,plotData[2][1:33:i]./1.5e11,label="",linewidth=2,linecolor=colors[1],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #plot orbits up to i
        p=plot!(plotData[3][1:33:i]./1.5e11,plotData[4][1:33:i]./1.5e11,label="",linewidth=2,linecolor=colors[2],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #linealpha argument causes lines to decay
        p=plot!(plotData[5][1:33:i]./1.5e11,plotData[6][1:33:i]./1.5e11,label="",linewidth=2,linecolor=colors[3],linealpha=max.((1:33:i) .+ 10000 .- i,2500)/10000) #example: alpha=max.((1:i) .+ 100 .- i,0) causes only last 100 to be visible
        p=scatter!(starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
        star1=makeCircleVals(rad[1],[plotData[1][i],plotData[2][i]]) #generate circles with appropriate sizes for each star
        star2=makeCircleVals(rad[2],[plotData[3][i],plotData[4][i]]) #at current positions
        star3=makeCircleVals(rad[3],[plotData[5][i],plotData[6][i]])
        p=plot!(star1[1]./1.5e11,star1[2]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true) #plot star circles with labels
        p=plot!(star2[1]./1.5e11,star2[2]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
        p=plot!(star3[1]./1.5e11,star3[2]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
        p=plot!(background_color=:black,background_color_legend=:transparent,foreground_color_legend=:transparent,
            background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white,legendfontfamily="Courier") #formatting for plot frame
        p=plot!(xlabel="x: AU",ylabel="y: AU",title="Random Three-Body Problem\nt:      years after start",
            legend=:best,xaxis=("x: AU",(limx[1],limx[2]),font(9,"Courier")),yaxis=("y: AU",(limy[1],limy[2]),font(9,"Courier")),tickfontcolor=:white,
            grid=false,titlefont=font(14,"Courier"),size=(720,721),legendfontsize=8,legendtitle="Mass (in solar masses)",legendtitlefontsize=8,legendtitlefont="Courier") #add in axes/title/legend with formatting

        tX,tY=relative(p[1],0.3,1.044)#static coords for time relative to parent
        p = annotate!(tX,tY,Plots.text((@sprintf("%0.2f",t[i]/365/24/3600)),"Courier",14,"black"))
        png(p,@sprintf("tmpPlots/frame_%06d.png",frameNum))
        frameNum+=1
        closeall() #close plots
    end

    if collisionBool==true #this condition makes 2 seconds of slo-mo right before the collision
        println("making collision cam")
        for i=1:10:600
            GR.inline("png") #added to eneable cron/jobber compatibility, also this makes frames generate WAY faster? Prior to adding this when run from cron/jobber frames would stop generating at 408 for some reason.
            gr(legendfontcolor = plot_color(:white)) #legendfontcolor=:white plot arg broken right now (at least in this backend)
            print("$(@sprintf("%.2f",i/600*100)) % complete\r") #output percent tracker
            pos=[plotData[1][end-(600-i)],plotData[2][end-(600-i)],plotData[3][end-(600-i)],plotData[4][end-(600-i)],plotData[5][end-(600-i)],plotData[6][end-(600-i)]] #current pos
            posFuture=pos #don't need future position at end
            limx,limy,center,orbitOld,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD=computeLimits(pos./1.5e11,posFuture./1.5e11,15,m,orbitOld,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD) #convert to AU, 10 AU padding
            p=plot(plotData[1][1:33:end-(600-i)]./1.5e11,plotData[2][1:33:end-(600-i)]./1.5e11,label="",linecolor=colors[1],linewidth=2,linealpha=max.((1:33:(i+length(t)-600)) .+ 10000 .- (i+length(t)-600),2500)/10000) #plot orbits up to i
            p=plot!(plotData[3][1:33:end-(600-i)]./1.5e11,plotData[4][1:33:end-(600-i)]./1.5e11,label="",linecolor=colors[2],linewidth=2,linealpha=max.((1:33:(i+length(t)-600)) .+ 10000 .- (i+length(t)-600),2500)/10000) #linealpha argument causes lines to decay
            p=plot!(plotData[5][1:33:end-(600-i)]./1.5e11,plotData[6][1:33:end-(600-i)]./1.5e11,label="",linecolor=colors[3],linewidth=2,linealpha=max.((1:33:(i+length(t)-600)) .+ 10000 .- (i+length(t)-600),2500)/10000) #example: alpha=max.((1:i) .+ 100 .- i,0) causes only last 100 to be visible
            p=scatter!(starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
            star1=makeCircleVals(rad[1],[plotData[1][end-(600-i)],plotData[2][end-(600-i)]]) #generate circles with appropriate sizes for each star
            star2=makeCircleVals(rad[2],[plotData[3][end-(600-i)],plotData[4][end-(600-i)]]) #at current positions
            star3=makeCircleVals(rad[3],[plotData[5][end-(600-i)],plotData[6][end-(600-i)]])
            p=plot!(star1[1]./1.5e11,star1[2]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true) #plot star circles with labels
            p=plot!(star2[1]./1.5e11,star2[2]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
            p=plot!(star3[1]./1.5e11,star3[2]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
            p=plot!(background_color=:black,background_color_legend=:transparent,foreground_color_legend=:transparent,
                background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white,legendfontfamily="Courier") #formatting for plot frame
            p=plot!(xlabel="x: AU",ylabel="y: AU",title="Random Three-Body Problem\nt: $(@sprintf("%0.2f",t[end-(600-i)]/365/24/3600)) years after start",
                legend=:best,xaxis=("x: AU",(limx[1],limx[2]),font(9,"Courier")),yaxis=("y: AU",(limy[1],limy[2]),font(9,"Courier")),
                grid=false,titlefont=font(14,"Courier"),size=(720,721),legendfontsize=8,legendtitle="Mass (in solar masses)",legendtitlefontsize=8,legendtitlefont="Courier") #add in axes/title/legend with formatting
            #collision cam zoom in
            i1,i2=collisionInds #these are the ones that are colliding, we use them to set the frame limits
            X=[plotData[1][end-(600-i)],plotData[3][end-(600-i)],plotData[5][end-(600-i)]]./1.5e11; Y=[plotData[2][end-(600-i)],plotData[4][end-(600-i)],plotData[6][end-(600-i)]]./1.5e11
            minX=(min(X[i1],X[i2])-1); maxX=(max(X[i1],X[i2])); minY=(min(Y[i1],Y[i2])-1); maxY=(max(Y[i1],Y[i2]))
            dx=maxX-minX; dy=maxY-minY
            dF=dx<dy ? dy : dx #use dy for frame if dx smaller, else dx

            #draw zoom box
            cornersX=[minX,minX+dF+1]; cornersY=[minY,minY+dF+1]
            p=plot!([cornersX[1],cornersX[2]],[cornersY[1],cornersY[1]],c=:white,label="") #side 1
            p=plot!([cornersX[2],cornersX[2]],[cornersY[1],cornersY[2]],c=:white,label="") #side 2
            p=plot!([cornersX[1],cornersX[2]],[cornersY[2],cornersY[2]],c=:white,label="") #side 3
            p=plot!([cornersX[1],cornersX[1]],[cornersY[1],cornersY[2]],c=:white,label="") #side 4
            offset = 0.0145 #for some reason the x corners don't quite match...
            s1x,s1y = relative(p[1],1/8-offset,7/8-0.25); s2x,s2y = relative(p[1],1/8+0.25-offset,7/8)
            subCornersX=[s1x,s2x]; subCornersY=[s1y,s2y] #physical coordinates, box in top left
            p=plot!([subCornersX[1],cornersX[1]],[subCornersY[1],cornersY[1]],c=:white,label = "") #corner 1 -> corner 1
            p=plot!([subCornersX[2],cornersX[2]],[subCornersY[2],cornersY[2]],c=:white,label = "") #corner 2 -> corner 2
            p=plot!([subCornersX[1],cornersX[1]],[subCornersY[2],cornersY[2]],c=:white,label = "") #corner 3 -> corner 3
            p=plot!([subCornersX[2],cornersX[2]],[subCornersY[1],cornersY[1]],c=:white,label = "") #corner 4 -> corner 4
            #draw box before plot so plot labels are on top
            p=plot!(title="COLLISION CAM\n(slo-mo x 33)",titlefontcolor=:orange,inset=(1,bbox(1/8,1/8,0.25,0.25)),
                    xlims=(minX,minX+dF+1),ylims=(minY,minY+dF+1),legend=:false,left_margin=0mm,right_margin=0mm,top_margin=0mm,bottom_margin=0mm,
                    foreground_color_border=:white,foreground_color_axis=:white,foreground_color_text=:white,grid=:false,
                    aspect_ratio=:equal,fontfamily=:Courier,subplot=2,framestyle=:box,titlefontsize=10,tickfontsize=6)
            p=plot!(p[2],plotData[1][1:10:end-(600-i)]./1.5e11,plotData[2][1:10:end-(600-i)]./1.5e11,label="",linecolor=colors[1],linewidth=2,linealpha=max.((1:10:(i+length(t)-600)) .+ 10000 .- (i+length(t)-600),2500)/10000)
            p=plot!(p[2],plotData[3][1:10:end-(600-i)]./1.5e11,plotData[4][1:10:end-(600-i)]./1.5e11,label="",linecolor=colors[2],linewidth=2,linealpha=max.((1:10:(i+length(t)-600)) .+ 10000 .- (i+length(t)-600),2500)/10000)
            p=plot!(p[2],plotData[5][1:10:end-(600-i)]./1.5e11,plotData[6][1:10:end-(600-i)]./1.5e11,label="",linecolor=colors[3],linewidth=2,linealpha=max.((1:10:(i+length(t)-600)) .+ 10000 .- (i+length(t)-600),2500)/10000)
            p=plot!(p[2],star1[1]./1.5e11,star1[2]./1.5e11,color=colors[1],fill=true)
            p=plot!(p[2],star2[1]./1.5e11,star2[2]./1.5e11,color=colors[2],fill=true)
            p=plot!(p[2],star3[1]./1.5e11,star3[2]./1.5e11,color=colors[3],fill=true)
            p=scatter!(p[2],starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
            #save frame
            png(p,@sprintf("tmpPlots/frame_%06d.png",frameNum))
            frameNum+=1
            closeall() #close plots
        end
        println("making freeze frame ending")
        for i=1:15 #make 0.5 s freeze frame ending
            GR.inline("png") #added to eneable cron/jobber compatibility, also this makes frames generate WAY faster? Prior to adding this when run from cron/jobber frames would stop generating at 408 for some reason.
            gr(legendfontcolor = plot_color(:white)) #legendfontcolor=:white plot arg broken right now (at least in this backend)
            print("$(@sprintf("%.2f",i/15*100)) % complete\r") #output percent tracker
            pos=[plotData[1][end],plotData[2][end],plotData[3][end],plotData[4][end],plotData[5][end],plotData[6][end]] #current pos
            posFuture=pos #don't need future position at end
            limx,limy,center,orbitOld,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD=computeLimits(pos./1.5e11,posFuture./1.5e11,15,m,orbitOld,ΔCx,ΔCy,ΔL,ΔR,ΔU,ΔD) #convert to AU, 10 AU padding
            p=plot(plotData[1][1:33:end]./1.5e11,plotData[2][1:33:end]./1.5e11,label="",linecolor=colors[1],linewidth=2,linealpha=max.((1:33:(length(t))) .+ 10000 .- (length(t)),2500)/10000) #plot orbits up to i
            p=plot!(plotData[3][1:33:end]./1.5e11,plotData[4][1:33:end]./1.5e11,label="",linecolor=colors[2],linewidth=2,linealpha=max.((1:33:(length(t))) .+ 10000 .- (length(t)),2500)/10000) #linealpha argument causes lines to decay
            p=plot!(plotData[5][1:33:end]./1.5e11,plotData[6][1:33:end]./1.5e11,label="",linecolor=colors[3],linewidth=2,linealpha=max.((1:33:(length(t))) .+ 10000 .- (length(t)),2500)/10000) #example: alpha=max.((1:i) .+ 100 .- i,0) causes only last 100 to be visible
            p=scatter!(starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
            star1=makeCircleVals(rad[1],[plotData[1][end],plotData[2][end]]) #generate circles with appropriate sizes for each star
            star2=makeCircleVals(rad[2],[plotData[3][end],plotData[4][end]]) #at current positions
            star3=makeCircleVals(rad[3],[plotData[5][end],plotData[6][end]])
            p=plot!(star1[1]./1.5e11,star1[2]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true) #plot star circles with labels
            p=plot!(star2[1]./1.5e11,star2[2]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
            p=plot!(star3[1]./1.5e11,star3[2]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
            p=plot!(background_color=:black,background_color_legend=:transparent,foreground_color_legend=:transparent,
                background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white,legendfontfamily="Courier") #formatting for plot frame
            p=plot!(xlabel="x: AU",ylabel="y: AU",title="Random Three-Body Problem\nt: $(@sprintf("%0.2f",t[end]/365/24/3600)) years after start",
                legend=:best,xaxis=("x: AU",(limx[1],limx[2]),font(9,"Courier")),yaxis=("y: AU",(limy[1],limy[2]),font(9,"Courier")),
                grid=false,titlefont=font(14,"Courier"),size=(720,721),legendfontsize=8,legendtitle="Mass (in solar masses)",legendtitlefontsize=8,legendtitlefont="Courier") #add in axes/title/legend with formatting
            #collision cam zoom in
            i1,i2=collisionInds #these are the ones that are colliding, we use them to set the frame limits
            X=[plotData[1][end],plotData[3][end],plotData[5][end]]./1.5e11; Y=[plotData[2][end],plotData[4][end],plotData[6][end]]./1.5e11
            minX=(min(X[i1],X[i2])-1); maxX=(max(X[i1],X[i2])); minY=(min(Y[i1],Y[i2])-1); maxY=(max(Y[i1],Y[i2]))
            dx=maxX-minX; dy=maxY-minY
            dF=dx<dy ? dy : dx #use dy for frame if dx smaller, else dx

            p=plot!(title="COLLISION CAM\n(slo-mo x 33)",titlefontcolor=:orange,inset=(1,bbox(1/8,1/8,0.25,0.25)),
                    xlims=(minX,minX+dF+1),ylims=(minY,minY+dF+1),legend=:false,left_margin=0mm,right_margin=0mm,top_margin=0mm,bottom_margin=0mm,
                    foreground_color_border=:white,foreground_color_axis=:white,foreground_color_text=:white,grid=:false,
                    aspect_ratio=:equal,fontfamily=:Courier,subplot=2,framestyle=:box,titlefontsize=10,tickfontsize=6)
            p=plot!(p[2],plotData[1][1:10:end]./1.5e11,plotData[2][1:10:end]./1.5e11,label="",linecolor=colors[1],linewidth=2,linealpha=max.((1:10:(length(t))) .+ 10000 .- (length(t)),2500)/10000)
            p=plot!(p[2],plotData[3][1:10:end]./1.5e11,plotData[4][1:10:end]./1.5e11,label="",linecolor=colors[2],linewidth=2,linealpha=max.((1:10:(length(t))) .+ 10000 .- (length(t)),2500)/10000)
            p=plot!(p[2],plotData[5][1:10:end]./1.5e11,plotData[6][1:10:end]./1.5e11,label="",linecolor=colors[3],linewidth=2,linealpha=max.((1:10:(length(t))) .+ 10000 .- (length(t)),2500)/10000)
            p=plot!(p[2],star1[1]./1.5e11,star1[2]./1.5e11,color=colors[1],fill=true)
            p=plot!(p[2],star2[1]./1.5e11,star2[2]./1.5e11,color=colors[2],fill=true)
            p=plot!(p[2],star3[1]./1.5e11,star3[2]./1.5e11,color=colors[3],fill=true)
            p=scatter!(p[2],starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
            #draw zoom box
            cornersX=[minX,minX+dF+1]; cornersY=[minY,minY+dF+1]
            p=plot!([cornersX[1],cornersX[2]],[cornersY[1],cornersY[1]],c=:white,label="") #side 1
            p=plot!([cornersX[2],cornersX[2]],[cornersY[1],cornersY[2]],c=:white,label="") #side 2
            p=plot!([cornersX[1],cornersX[2]],[cornersY[2],cornersY[2]],c=:white,label="") #side 3
            p=plot!([cornersX[1],cornersX[1]],[cornersY[1],cornersY[2]],c=:white,label="") #side 4
            offset = 0.0145 #for some reason the x corners don't quite match...
            s1x,s1y = relative(p[1],1/8-offset,7/8-0.25); s2x,s2y = relative(p[1],1/8+0.25-offset,7/8)
            subCornersX=[s1x,s2x]; subCornersY=[s1y,s2y] #physical coordinates, box in top left
            p=plot!([subCornersX[1],cornersX[1]],[subCornersY[1],cornersY[1]],c=:white,label = "") #corner 1 -> corner 1
            p=plot!([subCornersX[2],cornersX[2]],[subCornersY[2],cornersY[2]],c=:white,label = "") #corner 2 -> corner 2
            p=plot!([subCornersX[1],cornersX[1]],[subCornersY[2],cornersY[2]],c=:white,label = "") #corner 3 -> corner 3
            p=plot!([subCornersX[2],cornersX[2]],[subCornersY[1],cornersY[1]],c=:white,label = "") #corner 4 -> corner 4
            #save frame
            png(p,@sprintf("tmpPlots/frame_%06d.png",frameNum))
            frameNum+=1
            closeall() #close plots
        end
    end
end

function doOne()
    println("sit tight -- finding an interesting solution")

    #v1x=-1.673 v1y=-4.038 v2x=-4.276 v2y=6.851 v3x=2.399 v3y=0.942 (km/s)
    #x1=-7 y1=-11 x2=8 y2=-13 x3=19 y3=4 (AU from center)
    r=[-7,-11,8,-13,19,4,-1.673 * 1e3,-4.038 * 1e3,-4.276 * 1e3,6.851 * 1e3,2.399 * 1e3, 0.942 * 1e3]
    #m1=84.2 m2=26.9 m3=58.0 (solar masses)
    m=[84.2, 26.9, 58.0]
    rad=m.^0.8 #3 radii based on masses in solar units
    m=m.*2e30 #convert to SI kg
    rad=rad.*7e8 #convert to SI m

    initCond = [r,m,rad]
    initCond = nothing

    plotData,t,m,rad,collisionBool,collisionInds=getInteresting3Body(15, initCond) #find an interesting solution at least 15 years
    if collisionBool == true
        println("collision! inds = $collisionInds")
    else
        println("no collision")
    end

    #saveFrames(plotData,t,m,rad,collisionBool,collisionInds)

end

function mapCollisions()
    yearSec=365*24*3600

    xres = 800
    yres = 800
    xmin = -80
    xmax = 80
    ymin = -80
    ymax = 80

    data = rand(xres,yres)

    dfdata = DataFrame(xPos = Number[], yPos = Number[], timeToCrash = Number[])

    #Optional, preload dataframe from where we last left off
    #old_dfdata = DataFrame(CSV.File("tmpPlots/tmpSeparatedData_300y.csv"))
	old_dfdata = DataFrame(xPos = [], yPos= [])

    xvals = collect(1:xres) /  xres * (xmax - xmin) .+ xmin
    yvals = collect(1:yres) / yres * (ymax - ymin) .+ ymin

    for i = 1:xres
        x = xvals[i]
        for j = 1:yres
            y = yvals[j]
            #x1=-7 y1=-11 x2=8 y2=-13 x3=19 y3=4 (AU from center) #v1x=-1.673 v1y=-4.038 v2x=-4.276 v2y=6.851 v3x=2.399 v3y=0.942 (km/s)
            r=[-7*1.5e11,-11*1.5e11,8*1.5e11,-13*1.5e11,19*1.5e11,4*1.5e11,
                -1.673 * 1e3,-4.038 * 1e3,-4.276 * 1e3,6.851 * 1e3,2.399 * 1e3, 0.942 * 1e3]
            #m1=84.2 m2=26.9 m3=58.0 (solar masses)
            m=[84.2, 26.9, 58.0]

            rad=m.^0.8 #3 radii based on masses in solar units
            m=m.*2e30 #convert to SI kg
            rad=rad.*7e8 #convert to SI m

            #r[1] = (i / xres * (xmax - xmin) + xmin) * 1.5e11
            #r[2] = (j / yres * (ymax - ymin) + ymin) * 1.5e11

            r[1] = x * 1.5e11
            r[2] = y * 1.5e11

            println("starting coordinates: $x, $y")
            #println(r[1])
            #println(r[2])

            #Check if data for this point already exists (loading previous progress)
            alreadyexisting = @subset(old_dfdata, :xPos .== x, :yPos .== y)
            if nrow(alreadyexisting) > 0
                println("already existing data loaded")
                time = alreadyexisting[1, :timeToCrash]
                if time == -30
                    collisionBool = false
                else
                    collisionBool = true
                end
            else
                initCond = [r,m,rad]
                plotData,t,m,rad,collisionBool,collisionInds=getInteresting3Body(15, initCond) #find an interesting solution at least 15 years
                time = maximum(t)/yearSec
            end

            if collisionBool == true
                println("***********************")
                println("collision!")
                println("***********************")
                push!(dfdata, (x, y, time))
                data[i, j] = time
            else
                println("no collision, time $time")
                push!(dfdata, (x, y, -30))
                data[i, j] = -30
            end

        end
        CSV.write("tmpPlots/tmpSeparatedData.csv", dfdata)
    end

    CSV.write("tmpPlots/separatedData.csv", dfdata)

    df = DataFrame(data, :auto)
    open("tmpPlots/row-col-vals.txt","w") do f #save initial conditions to file in folder where script is run
        write(f,"rownames=$yvals \ncolvals=$xvals")
    end
    CSV.write("tmpPlots/heatmapdata.csv", df)

    gr()
    p=heatmap(xvals,
        yvals, data,
        c=cgrad([:blue, :white,:red, :yellow]),
        xlabel="x values", ylabel="y values",
        title="My title")

    png(p,@sprintf("tmpPlots/heatmap.png"))

end

function main() #pulls everything together, speeds things up to put everything in a function + gets rid of bad global syntax
    #doOne()

    mapCollisions()

end

#this is a function that will generate the animation for you without having to use the command line, works on Linux and Windows (run as administrator), untested on macOS
function makeAnim(clean=true)
    run(`ffmpeg -framerate 30 -i "tmpPlots/frame_%06d.png" -c:v libx264 -preset slow -coder 1 -movflags +faststart -g 15 -crf 18 -pix_fmt yuv420p -profile:v high -y -bf 2 -vf "scale=720:720,setdar=1/1" "threeBody.mp4"`)
    if clean==true
        println("cleaning up png files")
        foreach(rm,[string("tmpPlots/",x) for x in filter(endswith(".png"),readdir("tmpPlots"))])
    end
end

#generate frames!
main()

#generate the animation!
#makeAnim() #commented out because I compile the frames in the shell script (see 3BodyShell.sh)

#you'll want to uncomment the line above to tell the script to just make the animation
#or you can load the script into a julia terminal with include("threeBodyProb.jl")
#and then call main() and makeAnim() to your heart's content

#OLD STUFF....

#threeBodyFile="3Body_fps30.mp4"

#crf is compression value (17 or 18 "visually lossless"), pix_fmt is for twitter specific vid req, -b:v specifies target bitrate, -vcodec specifies codec (h264 in this case) -y says overwrite existing file
#run( `ffmpeg -framerate 30 -i $plotLoadPath"%06d.png" -vcodec libx264 -pix_fmt yuv420p -profile:v high -b:v 2048K -y -vf "scale=720:720,setdar=1/1" $threeBodyFile` ) #-vf scale=720:72 -crf 25
#run( `ffmpeg -framerate 30 -i $plotLoadPath"%06d.png" -c:v libx264 -preset slow -coder 1 -movflags +faststart -g 15 -crf 18 -pix_fmt yuv420p -profile:v high -y -bf 2 -vf "scale=720:720,setdar=1/1" $threeBodyFile` ) #all this bullshit to hopefully satisfy twitter requirements

#NOTE: moved ffmpeg commands to shell script

#old (simpler) way of generating animation
#uncomment and use this way if you just want a simple animation saved and don't
#care about performance/specific formatting of video.

# threeBodyAnim=@animate for i=1:length(t)
#     gr(legendfontcolor = plot_color(:white)) #plot arg broken right now in Julia
#     print("$(@sprintf("%.2f",i/length(t)*100)) % complete\r") #output percent tracker
#     pos=[plotData[1][i],plotData[2][i],plotData[3][i],plotData[4][i],plotData[5][i],plotData[6][i]] #current pos
#     limx,limy,center=getLims(pos./1.5e11,5) #convert to AU, 5 AU padding
#     plot(plotData[1][1:i]./1.5e11,plotData[2][1:i]./1.5e11,label="",linecolor=colors[1])
#     plot!(plotData[3][1:i]./1.5e11,plotData[4][1:i]./1.5e11,label="",linecolor=colors[2])
#     plot!(plotData[5][1:i]./1.5e11,plotData[6][1:i]./1.5e11,label="",linecolor=colors[3])
#     scatter!(starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
#     star1=makeCircleVals(rad[1],[plotData[1][i],plotData[2][i]])
#     star2=makeCircleVals(rad[2],[plotData[3][i],plotData[4][i]])
#     star3=makeCircleVals(rad[3],[plotData[5][i],plotData[6][i]])
#     plot!(star1[1]./1.5e11,star1[2]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true)
#     plot!(star2[1]./1.5e11,star2[2]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
#     plot!(star3[1]./1.5e11,star3[2]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
#     plot!(background_color=:black,background_color_legend=:transparent,background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white) #legendfontcolor=:white
#     plot!(xlabel="x: AU",ylabel="y: AU",title="Random Three Body Problem\nt: $(@sprintf("%0.2f",t[i]/365/24/3600)) yrs after start",
#         legend=:best,xaxis=("x: AU",(limx[1],limx[2]),font(12,"Courier")),yaxis=("y: AU",(limy[1],limy[2]),font(12,"Courier")),
#         grid=false,titlefont=font(24,"Courier"),size=(720,720),legendfontsize=12,legendtitle="Mass (in solar masses)",legendtitlefontsize=14)
#     end every 25

#mp4(threeBodyAnim,"3Body_fps30.mp4",fps=30)
#OR
#gif(threeBodyAnim,"3Body_fps30.gif",fps=30)


text:
# Taken from https://github.com/tk3369/YeoJohnsonTrans.jl/blob/master/src/YeoJohnsonTrans.jl

struct PowerTransform{T}
    λ::T
end

function yj(x, λ)
    y = float.(x)
    for (i, x) in enumerate(x)
        if x >= 0
            y[i] = λ ≈ 0 ? log(x + 1) : ((x + 1)^λ - 1)/λ
        else
            y[i] = λ ≈ 2 ? -log(-x + 1) : -((-x + 1)^(2 - λ) - 1) / (2 - λ)
        end
    end
    y
end

function (r::PowerTransform)(x)
    yj(x, r.λ)
end

mutable struct StepPowerTransform{S,F} <: AbstractStep
    selections::S
    params::Union{Nothing, Dict{Symbol,PowerTransform}}
    trained::Bool
    skip::Bool
    prehook::F
end

function log_likelihood(x, λ)
    N = length(x)
    y = yj(float.(x), λ)
    σ² = var(y, corrected = false)
    c = sum(sign.(x) .* log.(abs.(x) .+ one(eltype(x))))
    -N / 2.0 * log(σ²) + (λ - 1) * c
end

function lambda(x; interval = (-2.0, 3.0), optim_args...)
    i1, i2 = interval
    res = optimize(λ -> -log_likelihood(x, λ), i1, i2; optim_args...)
    round(0.1minimizer(res))/0.1
end

function transform!(s::StepPowerTransform, df)
    for (k,f) in s.params
        df[k] = f(df[k])
    end
    df
end
function transform(s::StepPowerTransform, df)
    transform!(s, copy(df))
end

function fit!(s::StepPowerTransform, df)
    s.params = Dict(col => PowerTransform(lambda(s.prehook(df[col]))) for col in getselectionkeys(df, s.selections))
    s.trained = true
end

function step_powertransform!(r::Recipe, s...; skip::Bool=false, prehook=identity)
    push!(r.steps,
          StepPowerTransform([s...],
                             nothing,
                             false,
                             skip,
                             prehook))
end

ext: jl


text:
include("DM1_1.jl")
include("DM2_1.jl")

function ReactiveGRASP(matrix,cost, n, m, ite, coupe,temps)
#    println("^^^^^^^^ Debut R-GRASP ^^^^^^^^")

    #initialisation
    (p,nb_iteration,z_cumul,zBest,zWorst) = intialiser(matrix,cost, n, m, ite, coupe)
    evol_p=Array{Float64}[]
#    println("valeurs initiales: ")
#    println("p: ", p)
#    println("nb_iteration: ", nb_iteration)
#    println("z_cumul: ", z_cumul)
#    println("zBest: ", zBest , "  zWorst: ", zWorst)

    t=time()
    while (time()-t <= temps)
        cpt = 1
        push!(evol_p,p)
        while (cpt <= ite)
            cpt = cpt+1

            prob=rand(Float64)

            alpha_choisit = choix_alpha(p,prob)
#            println("alpha_choisit: ", alpha_choisit, "  p: ", p)
            (SOL,z) = GRASP(cost, matrix, n, m, p[alpha_choisit])

            nb_iteration[alpha_choisit] += 1
            z_cumul[alpha_choisit] += z

            if z > zBest
                zBest = z
            end

            if z < zWorst
                zWorst = z
            end

#            println(" ============= \n ite: ",cpt)
        end
        println("zBest: ", zBest, "   zWorst: ", zWorst, "    nb_iteration : ", nb_iteration , "    p: ", p)
        recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)

    end
    println(evol_p)
end

#=======================================#

function intialiser(matrix,cost, n, m, ite, coupe)

    p = Float64[]
    for i in 1:coupe
        push!(p,i/coupe)
    end

    nb_iteration = zeros(Int64, length(p))

    z_cumul = zeros(Int64, length(p))

    (SOL,z) =Glouton(cost, matrix, n, m)

    zBest = z

    zWorst = z

    return(p,nb_iteration,z_cumul,zBest,zWorst)
end

#=======================================#

function recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
    q=zeros(Float64,length(p))

    somme_q=0
#    println("q: ", q)
#    println("nb_iteration: ", nb_iteration)

    for i in 1:length(p)
#        println("q: ", q, "   i:", i)
#        if nb_iteration[i]!=0
#            println("z_cumul[i]: ", z_cumul[i],"   nb_iteration[i]: ", nb_iteration[i] )
            moyenne=(z_cumul[i]/nb_iteration[i])
#            println("moyenne: ", moyenne)
#        else
#            moyenne=0
#        end
#        if moyenne>0
            q[i]=(moyenne-zWorst)/(zBest-zWorst)
#            println("q[i]: ", q[i])

#        else
#            q[i]= zWorst/(zBest-zWorst)
#        end
        somme_q+=q[i]

    end
#    println("q: ", q, "      somme_q", somme_q)

    for i in 1:length(q)
        p[i]=q[i]/somme_q
    end
#    println("p: ", p)

end

#=======================================#

function choix_alpha(p,prob)
    for i in 1:length(p)
        if prob < p[i]
            return(i)
        elseif i==length(p)
            return(i)
        end
    end
end


text:
function call(::Types.Unknown, args::Vector, storage, conn::IO)
    write(conn, "This command doesn't exist, you can check /help\n")
end


text:
#!/usr/bin/env julia

using Weave
using ArgParse

ap = ArgParseSettings("Weave Julia documents using Weave.jl",
                version = string(Pkg.installed("Weave")),
                add_version = true)

@add_arg_table ap begin
    "source"
        nargs = '+'
        help = "source document(s)"
        required = true
    "--doctype"
        default = :auto
        help = "output format"
    "--plotlib"
        arg_type = String
        default = "Gadfly"
        help = "output format"
    "--informat"
        default = :auto
        help = "output format"
    "--out_path"
        arg_type = String
        default = ":doc"
        help = "output directory"
    "--fig_path"
        arg_type = String
        default = "figures"
        help = "figure output directory"
    "--fig_ext"
        default = nothing
        help = "figure file format"
end

args = ArgParse.parse_args(ap)
source = args["source"]
delete!(args, "source")
args_col = []

args["plotlib"] == "nothing" && (args["plotlib"] = nothing)
#Check for special values of out_path

#args["informat"] == ":auto" && (args["informat"] = :auto)
#args["doctype"] == ":auto" && (args["informat"] = :auto)

if args["out_path"] == ":doc"
    args["out_path"] = :doc
elseif args["out_path"] == ":pwd"
    args["out_path"] = :pwd
end

for (key, val) in args
    push!(args_col, (parse(key), val))
end

for s=source
    weave(s; args_col...)
end


text:
using LinearAlgebra

export  ABtDecomposition, toCartesian, rank, size, *, getindex, ACA

struct ABtDecomposition{T} <: AbstractArray{T,2}
    A::Array{T,2}
    Bt::Array{T,2}
end

# required in order to inherit from AbstractArray
function ABtDecomposition(A::Array{T,2},Bt::Array{U,2}) where{T<:Number,U<:Number}
    if size(A,2) != size(Bt,1)
        throw(DomainError("size(A,2) != size(Bt,1). The input arrays should be compatible for the direct matrix-matrix multiplication."))
    end
    cT = promote_type(T,U)
    return ABtDecomposition{cT}(
        T == cT ? A  : convert(Array{cT,2},A),
        U == cT ? Bt : convert(Array{cT,2},Bt)
    )
end

# get cartesian index
function toCartesian(i,M,N)
    return ((i-1)%M + 1, i÷N + 1)
end

LinearAlgebra.:rank(A::ABtDecomposition) = size(A.A,2)

Base.size(ABt::ABtDecomposition) = (size(ABt.A,1),size(ABt.Bt,2))
Base.size(ABt::ABtDecomposition,i) = begin
    if i == 1
        return size(ABt.A,1)
    elseif i == 2
        return size(ABt.Bt,2)
    else
        throw(ArgumentError("Invalid dimension. Should be 1 or 2."))
    end
end

Base.:*(A::ABtDecomposition{T},B::Array{S}) where {T,S} = A.A*(A.Bt*B)
Base.:*(A::Array{T},B::ABtDecomposition{S}) where {T,S} = (A*B.A)*B.Bt
Base.:*(A::ABtDecomposition,B::ABtDecomposition) = ABtDecomposition(A.A,(A.Bt*B.A)*B.Bt)

Base.getindex(A::ABtDecomposition,i::Int) = sum(
    A.A[toCartesian(i,size(A.A,1),size(A.A,2))[1],:] .*
    A.Bt[:,toCartesian(i,size(A.Bt,1),size(A.Bt,2))[2]]
)
Base.getindex(A::ABtDecomposition,i::Int,j::Int) = sum(A.A[i,:] .* A.Bt[:,j])

function ACA(I, J, fct::Function, tol::T) where {T<:Real}
    #= Adaptive Cross Approximation algorithm for the
        low-rank approximation of full matrices. =#
    # functions returning the rows or the columns of the matrix
    row(i) = fct(I[i],J)
    col(j) = fct(I,J[j])

    Nr = length(I)
    Nc = length(J)
    # pivot lists
    Ir = collect(1:Nr)
    Ic = collect(1:Nc)

    # first row
    k  = 1
    Bt = reshape(row(Ir[k]),(1,Nc))
    deleteat!(Ir,k)

    k = argmax(abs.(Bt[Ic]))
    c = Ic[k]
    β = Bt[c]
    deleteat!(Ic,k)

    A = col(c)/β
    A = reshape(A,(Nr,1))

    aa = real(dot(A,A)); bb = real(dot(Bt,Bt)); res = aa*bb
    ε = sqrt(aa)*sqrt(bb)/sqrt(res)

    anp1 = A[:,1:1]

    flag = true
    n    = 1
    while ε > tol
        if n*(Nr+Nc) > Nr*Nc
            flag = false
            break
        end

        k = argmax(abs.(anp1[Ir]))
        r = Ir[k]
        deleteat!(Ir,k)

        bnp1 = reshape(row(r),(1,Nc))
        u    = A[r:r,:]
        bnp1 .-= u*Bt

        k = argmax(abs.(bnp1[Ic]))
        c = Ic[k]
        β = bnp1[c]
        deleteat!(Ic,k)

        anp1 = reshape(col(c),(Nr,1))
        v    = Bt[:,c:c]
        anp1 .= (anp1-A*v)/β

        aa = real(dot(anp1,anp1)); bb = real(dot(bnp1,bnp1))
        u = A'*anp1
        v = Bt*bnp1'
        ab = real(dot(u,v))

        res += 2*ab+aa*bb
        ε = sqrt(aa)*sqrt(bb)/sqrt(res)

        A  = hcat(A, anp1)
        Bt = vcat(Bt,bnp1)

        n += 1
    end
    return ABtDecomposition(A, Bt), flag
end

max_stars_repo_name: marcbakry/HierarchicalMatrices

max_issues_repo_name: marcbakry/HierarchicalMatrices

max_forks_repo_name: marcbakry/HierarchicalMatrices


text:
using Xpress, Dates

@enum DefinitionReadMode Ignore MainControl MainAttribute

FILE = joinpath(ENV["XPRESSDIR"], "include", "xprs.h")

if !isfile(FILE)
    error("File $FILE does not exist.")
end

function process(FILE)

lines = readlines(FILE)

out = """
#=

File automatically generated with script:
Xpress.jl/scripts/build_param_control_dicts.jl

Last build: $(now())

$("Optimizer version: $(Xpress.getversion())")

Banner from lib:
$(Xpress.getbanner())

Banner from header (xprs.h):
$(lines[7])
$(lines[8])

=#

"""

controls_start = " * control parameters for XPRSprob"
string_control = "/* String control parameters */"
double_control = "/* Double control parameters */"
integer_control = "/* Integer control parameters */"

attributes_start = " * attributes for XPRSprob "
string_attribute = "/* String attributes */"
double_attribute = "/* Double attributes */"
integer_attribute = "/* Integer attributes */"

the_end = " * control parameters for XPRSmipsolpool"

current_definiton = Ignore
current_type = Nothing

dict = Dict(
    i => Any[] for i in [(j,k)
        for j in [Int32, String, Float64], k in [MainAttribute, MainControl]]
)

for line in lines
    # @show line
    if startswith(line, the_end)
        break
    end
    if startswith(line, controls_start)
        current_definiton = MainControl
        continue
    elseif startswith(line, attributes_start)
        current_definiton = MainAttribute
        continue
    elseif startswith(line, string_control) || startswith(line, string_attribute)
        current_type = String
        continue
    elseif startswith(line, double_control) || startswith(line, double_attribute)
        current_type = Float64
        continue
    elseif startswith(line, integer_control) || startswith(line, integer_attribute)
        current_type = Int32
        continue
    elseif startswith(line, "#define") && current_type != Nothing && current_definiton != Ignore
        data = split(line)
        name = data[2]
        # @show data[3]
        value = tryparse(Int32, data[3])
        if value === nothing
            println("failed to parse $line")
            continue
        end
        push!(dict[current_type, current_definiton], (name, value))
        continue
    else
        continue
    end
end

is_first = true

order = [
    (String, MainControl),
    (Float64, MainControl),
    (Int32, MainControl),
    (String, MainAttribute),
    (Float64, MainAttribute),
    (Int32, MainAttribute),
]

for (tp, mode) in order
    vec = dict[(tp, mode)]
    if is_first
        is_first = false
    else
        out *= ")\n\n"
    end
    if tp == String && mode == MainAttribute
        out *= "const STRING_ATTRIBUTES = Dict{String, Int32}(\n"
    elseif tp == Int32 && mode == MainAttribute
        out *= "const INTEGER_ATTRIBUTES = Dict{String, Int32}(\n"
    elseif tp == Float64 && mode == MainAttribute
        out *= "const DOUBLE_ATTRIBUTES = Dict{String, Int32}(\n"
    elseif tp == String && mode == MainControl
        out *= "const STRING_CONTROLS = Dict{String, Int32}(\n"
    elseif tp == Int32 && mode == MainControl
        out *= "const INTEGER_CONTROLS = Dict{String, Int32}(\n"
    elseif tp == Float64 && mode == MainControl
        out *= "const DOUBLE_CONTROLS = Dict{String, Int32}(\n"
    end
    for (name, value) in vec
        out *= "    \"$(replace(string(name), "XPRS_"=>""))\" => $value,\n"
    end
end
out *= ")\n\n"

out *= """
const STRING_CONTROLS_VALUES = values(STRING_CONTROLS)

const DOUBLE_CONTROLS_VALUES = values(DOUBLE_CONTROLS)

const INTEGER_CONTROLS_VALUES = values(INTEGER_CONTROLS)

const STRING_ATTRBUTES_VALUES = values(STRING_ATTRIBUTES)

const DOUBLE_ATTRBUTES_VALUES = values(DOUBLE_ATTRIBUTES)

const INTEGER_ATTRBUTES_VALUES = values(INTEGER_ATTRIBUTES)
"""

open(joinpath(@__DIR__, "..", "src", "attributes_controls.jl"), "w") do f
    write(f, out)
end

return nothing
end

process(FILE);


text:
# # Getting started
# ExplainableAI.jl can be used on any classifier.
# In this first example, we will look at attributions on a LeNet5 model that was pretrained on MNIST.
#
# ### Loading the model
#md # !!! note
#md #
#md #     Outside of these docs, you should be able to load the model using
#md #     ```julia
#md #     using BSON: @load
#md #     @load "model.bson" model
#md #     ```

using ExplainableAI
using Flux
using BSON

model = BSON.load("../model.bson", @__MODULE__)[:model]

#md # !!! warning "Strip softmax"
#md #
#md #     For models with softmax activations on the output, it is necessary to call
#md #     ```julia
#md #     model = strip_softmax(model)
#md #     ```
#md #     before analyzing.

# ### Loading MNIST
# We use MLDatasets to load a single image from the MNIST dataset:
using MLDatasets
using ImageCore
using ImageIO
using ImageShow

index = 10
x, _ = MNIST(Float32, :test)[10]

convert2image(MNIST, x)

# By convention in Flux.jl, this input needs to be resized to WHCN format by adding a color channel and batch dimensions.
input = reshape(x, 28, 28, 1, :);

#md # !!! warning "Input format"
#md #
#md #     For any attribution of a model, ExplainableAI.jl assumes the batch dimension to be come last in the input.
#md #
#md #     For the purpose of heatmapping, the input is assumed to be in WHCN order
#md #     (width, height, channels, batch), which is Flux.jl's convention.

# ## Calling the analyzer
# We can now select an analyzer of our choice
# and call [`analyze`](@ref) to get an `Explanation`:
analyzer = LRP(model)
expl = analyze(input, analyzer);

# This `Explanation` bundles the following data:
# * `expl.attribution`: the analyzer's attribution
# * `expl.output`: the model output
# * `expl.neuron_selection`: the neuron index of used for the attribution
# * `expl.analyzer`: a symbol corresponding the used analyzer, e.g. `:LRP`

# Finally, we can visualize the `Explanation` through heatmapping:
heatmap(expl)

# Or get the same result by combining both analysis and heatmapping into one step:
heatmap(input, analyzer)

# ## Neuron selection
# By passing an additional index to our call to `analyze`, we can compute the attribution
# with respect to a specific output neuron.
# Let's see why the output wasn't interpreted as a 4 (output neuron at index 5)
heatmap(input, analyzer, 5)

# This heatmap shows us that the "upper loop" of the hand-drawn 9 has negative relevance
# with respect to the output neuron corresponding to digit 4!

#md # !!! note
#md #
#md #     The ouput neuron can also be specified when calling `analyze`:
#md #     ```julia
#md #     expl = analyze(img, analyzer, 5)
#md #     ```

# ## Input batches
# ExplainableAI also supports input batches:
batchsize = 100
xs, _ = MNIST(Float32, :test)[1:batchsize]
batch = reshape(xs, 28, 28, 1, :) # reshape to WHCN format
expl_batch = analyze(batch, analyzer);

# This will once again return a single `Explanation` `expl_batch` for the entire batch.
# Calling `heatmap` on `expl_batch` will detect the batch dimension and return a vector of heatmaps.
#
# Let's check if the digit at `index = 10` still matches.
hs = heatmap(expl_batch)
hs[index]

# JuliaImages' `mosaic` can be used to return a tiled view of all heatmaps:
mosaic(hs; nrow=10)

# We can also evaluate a batch with respect to a specific output neuron, e.g. for the digit zero at index `1`:
mosaic(heatmap(batch, analyzer, 1); nrow=10)

# ## Automatic heatmap presets
# Currently, the following analyzers are implemented:
#
# ```
# ├── Gradient
# ├── InputTimesGradient
# ├── SmoothGrad
# ├── IntegratedGradients
# └── LRP
#     ├── LRPZero
#     ├── LRPEpsilon
#     └── LRPGamma
# ```
#
# Let's try `InputTimesGradient`
analyzer = InputTimesGradient(model)
heatmap(input, analyzer)

# and `Gradient`
analyzer = Gradient(model)
heatmap(input, analyzer)

# As you can see, the function `heatmap` automatically applies common presets for each method.
#
# Since `InputTimesGradient` and LRP both compute attributions, their presets are similar.
# Gradient methods however are typically shown in grayscale.

# ## Custom heatmap settings
# We can partially or fully override presets by passing keyword arguments to `heatmap`:
using ColorSchemes
heatmap(expl; cs=ColorSchemes.jet)
#
heatmap(expl; reduce=:sum, rangescale=:extrema, cs=ColorSchemes.inferno)

# This also works with batches
mosaic(heatmap(expl_batch; rangescale=:extrema, cs=ColorSchemes.inferno); nrow=10)

# For the full list of keyword arguments, refer to the [`heatmap`](@ref) documentation.


text:
# Command: Get Timeouts
"""
	timeouts(session::Session)
Get the timeouts according to [webdriver2](https://w3c.github.io/webdriver/#get-timeouts).
Currently not supported by Selenium.
"""
function timeouts(session::Session)
    @unpack addr, id = session
    response =
        HTTP.get("$addr/session/$id/timeouts", [("Content-Type" => "application/json; charset=utf-8")])
    @assert response.status == 200
    JSON3.read(response.body).status
end


text:
trig1arg = [
(:acos, :(-1./sqrt(1-abs2(x))), (-1,1)),
(:acosd, :(-(180/pi)./sqrt(1-abs2(x))), (-1,1)),
(:acosh, :(1./sqrt(abs2(x)-1)), (1,Inf)),
(:acot, :(-1./(1+abs2(x))), (-Inf,Inf)),
(:acotd, :(-(180/pi)./(1+abs2(x))), (-Inf,Inf)),
(:acoth, :(1./(1-abs2(x))), (-Inf,-1,1,Inf)),
(:acsc, :(-1./sqrt(x.*x.*(x-1).*(x+1))), (-Inf,-1,1,Inf)),
(:acscd, :(-(180/pi)./sqrt(x.*x.*(x-1).*(x+1))), (-Inf,-1,1,Inf)),
(:acsch, :(-1./sqrt(x.^4+x.^2)), (-Inf,Inf)),
(:asec, :(1./sqrt(x.^4-x.^2)), (-Inf,-1,1,Inf)),
(:asecd, :((180/pi)./sqrt(x.^4-x.^2)), (-Inf,-1,1,Inf)),
(:asech, :(-1./sqrt(x.^2-x.^4)), (0,1)),
(:asin, :(1./sqrt(1-abs2(x))), (-1,1)),
(:asind, :((180/pi)./sqrt(1-abs2(x))), (-1,1)),
(:asinh, :(1./sqrt(1+abs2(x))), (-Inf,Inf)),
(:atan, :(1./(1+abs2(x))), (-Inf,Inf)),
(:atand, :((180/pi)./(1+abs2(x))), (-Inf,Inf)),
(:atanh, :(1./(1-abs2(x))), (-1,1)),
(:cos, :(-sin(x)), (-Inf,Inf)),
# cos_kernel: Not exported
(:cosc, :(-2y./x-sinc(x)*(pi^2)), (-Inf,Inf)),
(:cosd, :(-sind(x)*pi/180), (-Inf,Inf)),
(:cosh, :(sinh(x)), (-Inf,Inf)),
(:cospi, :(-sinpi(x).*pi), (-Inf,Inf)),
(:cot, :(-abs2(csc(x))), (-Inf,Inf)),
(:cotd, :(-abs2(cscd(x))*pi/180), (-Inf,Inf)),
(:coth, :(-abs2(csch(x))), (-Inf,Inf)),
(:csc, :(-y.*cot(x)), (-Inf,Inf)),
(:cscd, :(-y.*cotd(x)*pi/180), (-Inf,Inf)),
(:csch, :(-y.*coth(x)), (-Inf,Inf)),
# deg2rad_ext: Not exported
# mulpi_ext: Not exported
(:sec, :(y.*tan(x)), (-Inf,Inf)),
(:secd, :(y.*tand(x)*pi/180), (-Inf,Inf)),
(:sech, :(-y.*tanh(x)), (-Inf,Inf)),
(:sin, :(cos(x)), (-Inf,Inf)),
# sin_kernel: Not exported
(:sinc, :(cosc(x)), (-Inf,Inf)), # sin(πx)/(πx)
(:sind, :(cosd(x)*pi/180), (-Inf,Inf)),
(:sinh, :(cosh(x)), (-Inf,Inf)),
(:sinpi, :(cospi(x).*pi), (-Inf,Inf)),
(:tan, :(1+abs2(y)), (-Inf,Inf)),
(:tand, :((1+abs2(y))*pi/180), (-Inf,Inf)),
(:tanh, :(1-abs2(y)), (-Inf,Inf)),
]

for (f,g,r) in trig1arg
    @eval @primitive $f(x),dy,y  (dy.*($g))
    addtest1(f,r)
end


text:
@testset "Datasets Kernel Tests" begin
    df1 = createdummydatasetone()
    dataset1 = Persa.Dataset(df1)

    df2 = createdummydatasettwo()
    dataset2 = Persa.Dataset(df2, 10, 10)

    @testset "Dummy Tests" begin
        @test dataset1.users == 7
        @test dataset1.items == 6
        @test Persa.users(dataset1) == 7
        @test Persa.items(dataset1) == 6

        @test dataset2.users == 10
        @test dataset2.items == 10
        @test Persa.users(dataset2) == 10
        @test Persa.items(dataset2) == 10

        @test Persa.size(dataset1) == (7, 6)
        @test Persa.size(dataset2) == (10, 10)

        @test Persa.length(dataset1) == 35
        @test Persa.length(dataset2) == 35
    end

    @testset "Index Tests" begin
        @testset "Cartesian Index Tests" begin
            for i = 1:size(df1)[1]
                user = df1[!, :user][i]
                item = df1[!, :item][i]
                rating = df1[!, :rating][i]
                @test Persa.value(dataset1[user, item]) == rating
            end

            for i = 1:size(df2)[1]
                user = df2[!, :user][i]
                item = df2[!, :item][i]
                rating = df2[!, :rating][i]
                @test Persa.value(dataset2[user, item]) == rating
            end

            for i = 1:length(dataset1)
                (user, item, rating) = dataset1[i]
                @test dataset1[user, item] == rating
            end

            for i = 1:length(dataset2)
                (user, item, rating) = dataset2[i]
                @test dataset2[user, item] == rating
            end
        end

        @testset "Column Index Tests" begin
            for user in 1:Persa.users(dataset1)
                for (_, item, rating) in dataset1[user, :]
                    @test dataset1[user, item] == rating
                end
            end

            for item in 1:Persa.items(dataset1)
                for (user, _, rating) in dataset1[:, item]
                    @test dataset1[user, item] == rating
                end
            end

            for user in 1:Persa.users(dataset2)
                for (_, item, rating) in dataset2[user, :]
                    @test dataset2[user, item] == rating
                end
            end

            for item in 1:Persa.items(dataset2)
                for (user, _, rating) in dataset2[:, item]
                    @test dataset2[user, item] == rating
                end
            end
        end

        @testset "Limits Tests" begin
            (users, items) = Persa.size(dataset1)

            @test_throws ArgumentError dataset1[users + 1, items]
            @test_throws ArgumentError dataset1[users, items + 1]
            @test_throws ArgumentError dataset1[users + 1, items + 1]
            @test_throws ArgumentError dataset1[(users * items) + 1]
        end
    end

    @testset "Column Index Tests" begin
        for user=1:Persa.users(dataset1)
            @test length(dataset1[user, :]) == length(dataset1[user, 1:end])
        end

        for item=1:Persa.items(dataset1)
            @test length(dataset1[:, item]) == length(dataset1[1:end, item])
        end

        for user=1:Persa.users(dataset2)
            @test length(dataset2[user, :]) == length(dataset2[user, 1:end])
        end

        for item=1:Persa.items(dataset2)
            @test length(dataset2[:, item]) == length(dataset2[1:end, item])
        end

        start = 2
        last = Persa.users(dataset1) - 1

        for item=1:Persa.items(dataset1)
            elements = dataset1[start:last, item]
            for element in elements
                @test element[1] >= start && element[1] <= last
            end
        end

        start = 2
        last = Persa.items(dataset1) - 1

        for user=1:Persa.users(dataset1)
            elements = dataset1[user, start:last]
            for element in elements
                @test element[2] >= start && element[2] <= last
            end
        end

        start = 2
        last = Persa.users(dataset2) - 1

        for item=1:Persa.items(dataset2)
            elements = dataset2[start:last, item]
            for element in elements
                @test element[1] >= start && element[1] <= last
            end
        end

        start = 2
        last = Persa.items(dataset2) - 1

        for user=1:Persa.users(dataset2)
            elements = dataset2[user, start:last]
            for element in elements
                @test element[2] >= start && element[2] <= last
            end
        end

        @test length(dataset1[1:end, 1:end]) == length(dataset1)
        @test length(dataset2[1:end, 1:end]) == length(dataset2)

        @test dataset1[1] == dataset1[[1]][1]
        @test dataset1[1] in dataset1[[1]]
        @test dataset1[1] in dataset1[[1, 2]]
        @test length(dataset1[[1,2,3]]) == 3

        @test dataset1[1] == dataset1[1:1][1]
        @test dataset1[1] in dataset1[1:1]
        @test dataset1[1] in dataset1[1:2]
        @test length(dataset1[1:3]) == 3

        @test dataset1[1] == dataset1[:][1]
        @test dataset1[1] in dataset1[:]
        @test dataset1[1] in dataset1[:]
        @test length(dataset1) == length(dataset1[:])
        @test length(dataset1) == length(dataset1[:, :])
    end

    @testset "Matrix Conversion Tests" begin
        matrix = Array(dataset1)

        for user in 1:Persa.users(dataset1)
            for item in 1:Persa.items(dataset1)
                @test Persa.value(dataset1[user, item]) === matrix[user, item]
            end
        end

        matrix = Array(dataset2)

        for user in 1:Persa.users(dataset2)
            for item in 1:Persa.items(dataset2)
                @test Persa.value(dataset2[user, item]) === matrix[user, item]
            end
        end
    end

    @testset "Iterator Tests" begin
        for (user, item, rating) in dataset1
            @test rating == dataset1[user, item]
        end

        for (user, item, rating) in dataset2
            @test dataset2[user, item] == rating
        end
    end

    @testset "Clone Tests" begin
        cloned = Persa.sample(dataset1, [1:length(dataset1)...])

        @test length(dataset1) == length(cloned)
        @test size(dataset1) == size(cloned)
        @test Persa.users(dataset1) == Persa.users(cloned)
        @test Persa.items(dataset1) == Persa.items(cloned)

        cloned = Persa.sample(dataset2, [1:length(dataset2)...])

        @test length(dataset2) == length(cloned)
        @test size(dataset2) == size(cloned)
        @test Persa.users(dataset2) == Persa.users(cloned)
        @test Persa.items(dataset2) == Persa.items(cloned)
    end
end


text:
module Systems

using DocStringExtensions

export NavierStokes, PointForce, SpatialGauss, Gaussian!, Gaussian, vorticity, velocity, streamfunction, nl

using ..Fields
using ..TimeMarching
using ..RigidBodyMotions
using ..Bodies
#import ViscousFlow: r₁
#import ViscousFlow: plan_intfact

using LinearAlgebra
using SparseArrays

"""
    assign_velocity!(V::VectorData,X::VectorData,
                     xc::Real,yc::Real,α::Real,
                     mlist::Vector{RigidBodyMotion},t::Real)

Assign the components of rigid body velocity for every body (in inertial coordinate system)
at time `t` in the overall data structure `V`, using coordinates described by `X` (also in inertial
coordinate system), based on array of supplied motion `mlist` for each body.
"""
function RigidBodyMotions.assign_velocity!(V::VectorData{N},X::VectorData{N},
                                           bl::BodyList,tlist::Vector{RigidTransform},
                                           mlist::Vector{RigidBodyMotion},t::Real) where {N}
    N == numpts(bl) || error("Inconsistent size of data structures")
    for i in 1:length(bl)
        ui = view(V.u,bl,i)
        vi = view(V.v,bl,i)
        xi = view(X.u,bl,i)
        yi = view(X.v,bl,i)
        Ti = tlist[i]
        assign_velocity!(ui,vi,xi,yi,Ti.trans[1],Ti.trans[2],Ti.α,mlist[i],t)
    end
end

include("systems/rigidbodies.jl")
include("systems/navier_stokes.jl")
include("systems/histories.jl")
include("systems/utils.jl")

end


text:
function map_socs!(v::Vector{Float64}, cones::ConicSets, a::AuxiliaryData)
    cont = 0
    for (idx, sdp) in enumerate(cones.sdpcone)
        cont += div(sdp.sq_side * (sdp.sq_side + 1), 2)
    end
    sizehint!(a.soc_v, length(cones.socone))
    sizehint!(a.soc_s, length(cones.socone))
    for (idx, soc) in enumerate(cones.socone)
        len = soc.len
        push!(a.soc_s, view(v, cont + 1))
        push!(a.soc_v, view(v, cont + 2:cont + len))
        cont += len
    end
    return nothing
end

function ivech!(out::AbstractMatrix{T}, v::AbstractVector{T}) where T
    n = sympackeddim(length(v))
    n1, n2 = size(out)
    @assert n == n1 == n2
    c = 0
    for j in 1:n, i in 1:j
        c += 1
        out[i,j] = v[c]
    end
    return out
end

function ivech(v::AbstractVector{T}) where T
    n = sympackeddim(length(v))
    out = zeros(n, n)
    ivech!(out, v)

    return out
end

ivec(X) = Matrix(LinearAlgebra.Symmetric(ivech(X),:U))


text:
@testset "ConditionalMvNormal" begin

    xlength = 3
    zlength = 2
    batchsize = 10

    σvector(x::AbstractVector) = ones(Float32,xlength) .* 3
    σvector(x::AbstractMatrix) = ones(Float32,xlength,size(x,2)) .* 3
    σscalar(x::AbstractVector) = 2
    σscalar(x::AbstractMatrix) = ones(Float32,size(x,2)) .* 2

    x = rand(Float32, xlength) |> gpu
    z = rand(Float32, zlength) |> gpu
    X = rand(Float32, xlength, batchsize) |> gpu
    Z = rand(Float32, zlength, batchsize) |> gpu

    cases = [
        ("vector μ / vector σ",
         SplitLayer(zlength, [xlength,xlength], [identity,abs]), Vector, 4),
        ("vector μ / scalar σ",
         SplitLayer(zlength, [xlength,1], [identity,abs]), Real, 4),
        ("vector μ / fixed vector σ",
         SplitLayer(Dense(zlength,xlength), σvector), Vector, 2),
        ("vector μ / fixed scalar σ",
         SplitLayer(Dense(zlength,xlength), σscalar), Real, 2),
        ("vector μ / unit σ",
         Dense(zlength,xlength), Real, 2),
        ("vector μ / shared, trainable vector σ",
         SplitLayer(Dense(zlength,xlength), ones(Float32,xlength)), Vector, 3),
        ("vector μ / shared, trainable scalar σ",
         SplitLayer(Dense(zlength,xlength), 1f0), Real, 3)
    ]

    disttypes(::Type{<:Vector}) = (TuringDiagMvNormal,ConditionalDists.BatchDiagMvNormal)
    disttypes(::Type{<:Real}) = (TuringScalMvNormal,ConditionalDists.BatchScalMvNormal)

    for (name,mapping,T,nrps) in cases
        @testset "$name" begin
            p = ConditionalMvNormal(mapping) |> gpu
            (Texample,Tbatch) = disttypes(T)

            res = condition(p,z)
            μ = mean(res)
            σ2 = var(res)
            @test res isa Texample
            @test size(μ) == (xlength,)
            @test size(σ2) == (xlength,)

            loss() = logpdf(p,x,z)
            ps = Flux.params(p)
            @test length(ps) == nrps
            @test loss() isa Float32
            @test_nowarn Flux.gradient(loss, ps)

            f() = sum(rand(p,z))
            gs = Flux.gradient(f,ps)
            for p in ps
                g = gs[p]
                @test all(isfinite.(g)) && all(g .!= 0)
            end

            # batch tests
            res = condition(p,Z)
            μ = mean(res)
            σ2 = var(res)
            @test res isa Tbatch
            @test size(μ) == (xlength,batchsize)
            @test size(σ2) == (xlength,batchsize)

            loss() = sum(logpdf(p,X,Z))
            @test loss() isa Float32
            @test_nowarn Flux.gradient(loss, ps)

            f() = sum(rand(p,Z))
            gs = Flux.gradient(f,ps)
            for p in ps
                g = gs[p]
                @test all(isfinite.(g)) && all(g .!= 0)
            end
        end
    end
end


text:
MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{<:_C_SETS}) = true

function MOI.is_valid(opt::Optimizer, c::CI{SAF,T}) where {T<:_C_SETS}
    return 1 <= c.value <= num_constraints(opt.form, CCI{T})
end

function MOI.add_constraint(opt::Optimizer, f::SAF, s::T) where {T<:_C_SETS}
    if !iszero(f.constant)
        throw(MOI.ScalarFunctionConstantNotZero{Float64,SAF,T}(f.constant))
    end
    add_constraint(opt.form, AffineConstraint(f, s))
    return CI{SAF,T}(num_constraints(opt.form, CCI{T}))
end

function MOI.get(
    opt::Optimizer,
    ::MOI.NumberOfConstraints{SAF,S},
) where {S<:_C_SETS}
    return num_constraints(opt.form, CCI{S})
end

function MOI.get(
    opt::Optimizer,
    ::MOI.ListOfConstraintIndices{SAF,S},
) where {S<:_C_SETS}
    return [CI{SAF,S}(i) for i in 1:num_constraints(opt.form, CCI{S})]
end

function MOI.get(opt::Optimizer, ::MOI.ListOfConstraints)
    return error("TODO: Implement!")
end

MOI.supports(::Optimizer, ::MOI.ConstraintPrimal, ::CI{SAF,<:_C_SETS}) = true
function MOI.get(
    opt::Optimizer,
    ::MOI.ConstraintPrimal,
    ci::CI{SAF,S},
) where {S<:_C_SETS}
    MOI.throw_if_not_valid(opt, ci)
    if opt.result === nothing
        error("Model does not have a solution available.")
    end
    cci = convert(CCI{S}, ci)
    ac = get_constraint(opt.form, cci)
    val = 0.0
    for (coeff, cvi) in zip(ac.f.coeffs, ac.f.indices)
        val += coeff * opt.result.best_solution[index(cvi)]
    end
    return val
end


text:
module MadNLPTests

import LinearAlgebra: norm
import JuMP: Model, @variable, @constraint, @objective, @NLconstraint , @NLobjective, optimize!,
    MOI, termination_status, LowerBoundRef, UpperBoundRef, value, dual
import Test: @test, @testset

function solcmp(x,sol;atol=1e-4,rtol=1e-4)
    aerr = norm(x-sol,Inf)
    rerr = aerr/norm(sol,Inf)
    return (aerr < atol || rerr < rtol)
end

function test_madnlp(name,optimizer_constructor::Function,exclude)
    @testset "$name" begin
        for f in [infeasible,unbounded,lootsma,eigmina]
            !(string(f) in exclude) && f(optimizer_constructor)
        end
    end
end

function infeasible(optimizer_constructor::Function)
    @testset "infeasible" begin
        m=Model(optimizer_constructor)
        @variable(m,x>=1)
        @constraint(m,x==0.)
        @objective(m,Min,x^2)
        optimize!(m)
        @test termination_status(m) == MOI.LOCALLY_INFEASIBLE
    end
end

function unbounded(optimizer_constructor::Function)
    @testset "unbounded" begin
        m=Model(optimizer_constructor)
        @variable(m,x,start=1)
        @objective(m,Max,x^2)
        optimize!(m)
        @test termination_status(m) == MOI.INFEASIBLE_OR_UNBOUNDED
    end
end

function lootsma(optimizer_constructor::Function)
    @testset "lootsma" begin
        m=Model(optimizer_constructor)
        @variable(m, par == 6.)
        @variable(m,0 <= x[i=1:3] <= 5, start = 0.)
        l=[
            @NLconstraint(m,-sqrt(x[1]) - sqrt(x[2]) + sqrt(x[3]) >= 0.)
            @NLconstraint(m,sqrt(x[1]) + sqrt(x[2]) + sqrt(x[3]) >= 4.)
        ]
        @NLobjective(m,Min,x[1]^3 + 11. *x[1] - par*sqrt(x[1])  +x[3] )

        optimize!(m)

        @test solcmp(value.(x),[0.07415998565403112,2.9848713863700236,4.0000304145340415])
        @test solcmp(dual.(l),[2.000024518601535,2.0000305441119535])
        @test solcmp(dual.(LowerBoundRef.(x)),[0.,0.,0.])
        @test solcmp(dual.(UpperBoundRef.(x)),[0.,0.,0.])

        @test termination_status(m) == MOI.LOCALLY_SOLVED
    end
end

function eigmina(optimizer_constructor::Function)
    @testset "eigmina" begin
        m=Model(optimizer_constructor)
        @variable(m,-1 <= x[1:101] <= 1,start = .1)
        @constraint(m, x[1]*x[1] + x[2]*x[2] + x[3]*x[3] + x[4]*x[4] + x[5]*x[5] + x[6]*x[6] +
            x[7]*x[7] + x[8]*x[8] + x[9]*x[9] + x[10]*x[10] + x[11]*x[11] + x[12]*x[12] +
            x[13]*x[13] + x[14]*x[14] + x[15]*x[15] + x[16]*x[16] + x[17]*x[17] + x[18]*x[18] +
            x[19]*x[19] + x[20]*x[20] + x[21]*x[21] + x[22]*x[22] + x[23]*x[23] + x[24]*x[24] +
            x[25]*x[25] + x[26]*x[26] + x[27]*x[27] + x[28]*x[28] + x[29]*x[29] + x[30]*x[30] +
            x[31]*x[31] + x[32]*x[32] + x[33]*x[33] + x[34]*x[34] + x[35]*x[35] + x[36]*x[36] +
            x[37]*x[37] + x[38]*x[38] + x[39]*x[39] + x[40]*x[40] + x[41]*x[41] + x[42]*x[42] +
            x[43]*x[43] + x[44]*x[44] + x[45]*x[45] + x[46]*x[46] + x[47]*x[47] + x[48]*x[48] +
            x[49]*x[49] + x[50]*x[50] + x[51]*x[51] + x[52]*x[52] + x[53]*x[53] + x[54]*x[54] +
            x[55]*x[55] + x[56]*x[56] + x[57]*x[57] + x[58]*x[58] + x[59]*x[59] + x[60]*x[60] +
            x[61]*x[61] + x[62]*x[62] + x[63]*x[63] + x[64]*x[64] + x[65]*x[65] + x[66]*x[66] +
            x[67]*x[67] + x[68]*x[68] + x[69]*x[69] + x[70]*x[70] + x[71]*x[71] + x[72]*x[72] +
            x[73]*x[73] + x[74]*x[74] + x[75]*x[75] + x[76]*x[76] + x[77]*x[77] + x[78]*x[78] +
            x[79]*x[79] + x[80]*x[80] + x[81]*x[81] + x[82]*x[82] + x[83]*x[83] + x[84]*x[84] +
            x[85]*x[85] + x[86]*x[86] + x[87]*x[87] + x[88]*x[88] + x[89]*x[89] + x[90]*x[90] +
            x[91]*x[91] + x[92]*x[92] + x[93]*x[93] + x[94]*x[94] + x[95]*x[95] + x[96]*x[96] +
            x[97]*x[97] + x[98]*x[98] + x[99]*x[99] + x[100]*x[100] == 1)
        @constraint(m, x[1]*x[101] - x[1] == 0)
        @constraint(m, x[2]*x[101] - 2*x[2] == 0)
        @constraint(m, x[3]*x[101] - 3*x[3] == 0)
        @constraint(m, x[4]*x[101] - 4*x[4] == 0)
        @constraint(m, x[5]*x[101] - 5*x[5] == 0)
        @constraint(m, x[6]*x[101] - 6*x[6] == 0)
        @constraint(m, x[7]*x[101] - 7*x[7] == 0)
        @constraint(m, x[8]*x[101] - 8*x[8] == 0)
        @constraint(m, x[9]*x[101] - 9*x[9] == 0)
        @constraint(m, x[10]*x[101] - 10*x[10] == 0)
        @constraint(m, x[11]*x[101] - 11*x[11] == 0)
        @constraint(m, x[12]*x[101] - 12*x[12] == 0)
        @constraint(m, x[13]*x[101] - 13*x[13] == 0)
        @constraint(m, x[14]*x[101] - 14*x[14] == 0)
        @constraint(m, x[15]*x[101] - 15*x[15] == 0)
        @constraint(m, x[16]*x[101] - 16*x[16] == 0)
        @constraint(m, x[17]*x[101] - 17*x[17] == 0)
        @constraint(m, x[18]*x[101] - 18*x[18] == 0)
        @constraint(m, x[19]*x[101] - 19*x[19] == 0)
        @constraint(m, x[20]*x[101] - 20*x[20] == 0)
        @constraint(m, x[21]*x[101] - 21*x[21] == 0)
        @constraint(m, x[22]*x[101] - 22*x[22] == 0)
        @constraint(m, x[23]*x[101] - 23*x[23] == 0)
        @constraint(m, x[24]*x[101] - 24*x[24] == 0)
        @constraint(m, x[25]*x[101] - 25*x[25] == 0)
        @constraint(m, x[26]*x[101] - 26*x[26] == 0)
        @constraint(m, x[27]*x[101] - 27*x[27] == 0)
        @constraint(m, x[28]*x[101] - 28*x[28] == 0)
        @constraint(m, x[29]*x[101] - 29*x[29] == 0)
        @constraint(m, x[30]*x[101] - 30*x[30] == 0)
        @constraint(m, x[31]*x[101] - 31*x[31] == 0)
        @constraint(m, x[32]*x[101] - 32*x[32] == 0)
        @constraint(m, x[33]*x[101] - 33*x[33] == 0)
        @constraint(m, x[34]*x[101] - 34*x[34] == 0)
        @constraint(m, x[35]*x[101] - 35*x[35] == 0)
        @constraint(m, x[36]*x[101] - 36*x[36] == 0)
        @constraint(m, x[37]*x[101] - 37*x[37] == 0)
        @constraint(m, x[38]*x[101] - 38*x[38] == 0)
        @constraint(m, x[39]*x[101] - 39*x[39] == 0)
        @constraint(m, x[40]*x[101] - 40*x[40] == 0)
        @constraint(m, x[41]*x[101] - 41*x[41] == 0)
        @constraint(m, x[42]*x[101] - 42*x[42] == 0)
        @constraint(m, x[43]*x[101] - 43*x[43] == 0)
        @constraint(m, x[44]*x[101] - 44*x[44] == 0)
        @constraint(m, x[45]*x[101] - 45*x[45] == 0)
        @constraint(m, x[46]*x[101] - 46*x[46] == 0)
        @constraint(m, x[47]*x[101] - 47*x[47] == 0)
        @constraint(m, x[48]*x[101] - 48*x[48] == 0)
        @constraint(m, x[49]*x[101] - 49*x[49] == 0)
        @constraint(m, x[50]*x[101] - 50*x[50] == 0)
        @constraint(m, x[51]*x[101] - 51*x[51] == 0)
        @constraint(m, x[52]*x[101] - 52*x[52] == 0)
        @constraint(m, x[53]*x[101] - 53*x[53] == 0)
        @constraint(m, x[54]*x[101] - 54*x[54] == 0)
        @constraint(m, x[55]*x[101] - 55*x[55] == 0)
        @constraint(m, x[56]*x[101] - 56*x[56] == 0)
        @constraint(m, x[57]*x[101] - 57*x[57] == 0)
        @constraint(m, x[58]*x[101] - 58*x[58] == 0)
        @constraint(m, x[59]*x[101] - 59*x[59] == 0)
        @constraint(m, x[60]*x[101] - 60*x[60] == 0)
        @constraint(m, x[61]*x[101] - 61*x[61] == 0)
        @constraint(m, x[62]*x[101] - 62*x[62] == 0)
        @constraint(m, x[63]*x[101] - 63*x[63] == 0)
        @constraint(m, x[64]*x[101] - 64*x[64] == 0)
        @constraint(m, x[65]*x[101] - 65*x[65] == 0)
        @constraint(m, x[66]*x[101] - 66*x[66] == 0)
        @constraint(m, x[67]*x[101] - 67*x[67] == 0)
        @constraint(m, x[68]*x[101] - 68*x[68] == 0)
        @constraint(m, x[69]*x[101] - 69*x[69] == 0)
        @constraint(m, x[70]*x[101] - 70*x[70] == 0)
        @constraint(m, x[71]*x[101] - 71*x[71] == 0)
        @constraint(m, x[72]*x[101] - 72*x[72] == 0)
        @constraint(m, x[73]*x[101] - 73*x[73] == 0)
        @constraint(m, x[74]*x[101] - 74*x[74] == 0)
        @constraint(m, x[75]*x[101] - 75*x[75] == 0)
        @constraint(m, x[76]*x[101] - 76*x[76] == 0)
        @constraint(m, x[77]*x[101] - 77*x[77] == 0)
        @constraint(m, x[78]*x[101] - 78*x[78] == 0)
        @constraint(m, x[79]*x[101] - 79*x[79] == 0)
        @constraint(m, x[80]*x[101] - 80*x[80] == 0)
        @constraint(m, x[81]*x[101] - 81*x[81] == 0)
        @constraint(m, x[82]*x[101] - 82*x[82] == 0)
        @constraint(m, x[83]*x[101] - 83*x[83] == 0)
        @constraint(m, x[84]*x[101] - 84*x[84] == 0)
        @constraint(m, x[85]*x[101] - 85*x[85] == 0)
        @constraint(m, x[86]*x[101] - 86*x[86] == 0)
        @constraint(m, x[87]*x[101] - 87*x[87] == 0)
        @constraint(m, x[88]*x[101] - 88*x[88] == 0)
        @constraint(m, x[89]*x[101] - 89*x[89] == 0)
        @constraint(m, x[90]*x[101] - 90*x[90] == 0)
        @constraint(m, x[91]*x[101] - 91*x[91] == 0)
        @constraint(m, x[92]*x[101] - 92*x[92] == 0)
        @constraint(m, x[93]*x[101] - 93*x[93] == 0)
        @constraint(m, x[94]*x[101] - 94*x[94] == 0)
        @constraint(m, x[95]*x[101] - 95*x[95] == 0)
        @constraint(m, x[96]*x[101] - 96*x[96] == 0)
        @constraint(m, x[97]*x[101] - 97*x[97] == 0)
        @constraint(m, x[98]*x[101] - 98*x[98] == 0)
        @constraint(m, x[99]*x[101] - 99*x[99] == 0)
        @constraint(m, x[100]*x[101] - 100*x[100] == 0)
        @objective(m, Min, x[101])
        optimize!(m)

        @test termination_status(m) == MOI.LOCALLY_SOLVED
    end
end

export test_madnlp, solcmp

end # module


text:
_ndims(::Type{CartesianIndex{N}}) where N = N
_ndims(::Type{Int64}) = 1

"""
	generate_spatial_bounding_box(img, select_function)

Generate the spatial bounding box of foreground in the image with start-end positions.
Users can define arbitrary function to select expected foreground from the whole image or specified channels.
And it can also add margin to every dim of the bounding box.
"""
function generate_spatial_bounding_box(img, select_function)
	idx = findall(select_function, img)
	N = _ndims(eltype(idx))
	start_bounds = map(d -> minimum(getindex.(idx, d)), 1:N)
	end_bounds = map(d -> maximum(getindex.(idx, d)), 1:N)
	return start_bounds, end_bounds
end

function generate_spatial_bounding_box(img, select_function, channel_indices)
	idx = findall(select_function, img[channel_indices...])
	N = _ndims(eltype(idx))
	start_bounds = map(d -> minimum(getindex.(idx, d)), 1:N)
	end_bounds = map(d -> maximum(getindex.(idx, d)), 1:N)
	return start_bounds, end_bounds
end


text:
# NOTE: based on test/pkg.jl::capture_stdout, but doesn't discard exceptions
macro grab_output(ex)
    quote
        let fname = tempname()
            try
                ret = nothing
                open(fname, "w") do fout
                    redirect_stdout(fout) do
                        ret = $(esc(ex))
                    end
                end
                ret, read(fname, String)
            finally
                rm(fname, force=true)
            end
        end
    end
end

# variant on @test_throws that checks the CuError error code
macro test_throws_cuerror(code, ex)
    # generate a test only returning CuError if it is the correct one
    test = quote
        try
            $(esc(ex))
        catch err
            isa(err, CuError) || rethrow()
            err.code == $code || error("Wrong CuError: ", err.code, " instead of ", $code)
            rethrow()
        end
    end

    # now re-use @test_throws (which ties into @testset, etc)
    quote
        @test_throws CuError $test
    end
end

mutable struct NoThrowTestSet <: Test.AbstractTestSet
    results::Vector
    NoThrowTestSet(desc) = new([])
end
Test.record(ts::NoThrowTestSet, t::Test.Result) = (push!(ts.results, t); t)
Test.finish(ts::NoThrowTestSet) = ts.results
fails = @testset NoThrowTestSet begin
    # OK
    @test_throws_cuerror CUDAdrv.ERROR_UNKNOWN throw(CuError(CUDAdrv.ERROR_UNKNOWN))
    # Fail, wrong CuError
    @test_throws_cuerror CUDAdrv.ERROR_UNKNOWN throw(CuError(CUDAdrv.ERROR_INVALID_VALUE))
    # Fail, wrong Exception
    @test_throws_cuerror CUDAdrv.ERROR_UNKNOWN error()
end
@test isa(fails[1], Test.Pass)
@test isa(fails[2], Test.Fail)
@test isa(fails[3], Test.Fail)


text:
using NeXLMatrixCorrection
using Test

@testset "aspure" begin
    krp = aspure(KRatio(
        [n"O K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        mat"SiO2",
        uv(0.605, 0.01),
    ))
    @test isapprox(value(krp), 0.1302, atol=0.0001)

    krp = aspure(KRatio(
        [n"Si K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        mat"SiO2",
        uv(0.139188, 0.01),
    ))
    @test isapprox(value(krp), 0.0529, atol=0.0001)
    @test isapprox(σ(krp), 0.0038, atol=0.0001)

    krp = aspure(KRatio(
        [n"Ca K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        mat"CaF2",
        uv(0.193347, 0.001),
    ))
    @test isapprox(value(krp), 0.0962, atol=0.0001)
    @test isapprox(σ(krp), 0.0005, atol=0.0001)

    krps = aspure(KRatios(
        [n"O K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        mat"SiO2",
        [ 0.605 , 0.547978 ],
    ))
    @test isapprox(value(krps[1]), 0.1302, atol = 0.0001)
    @test isapprox(value(krps[2]), 0.1179, atol = 0.0001)
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
module MyPackage

using LinearAlgebra: norm

include("myfile.jl")

end


text:
@testset "Splitting Test" begin

    net_file = "$(@__DIR__)/../examples/networks/R2_R2.nnet"

    net = read_nnet(net_file, last_layer_activation = Id())
    A = [1. 0.; 0. 1.; -1. 0.; 0. -1.;]
    b = [1.1, 1.1, 0.1, 0.1]
    X = Hyperrectangle(low=[-1,-1], high=[1,1])
    # The true output set is a union of two segments : (0,0)->(0,1) U (0,0)->(1,0)
    Y = HPolytope(A, b) # Output constraints. Hyperrectangle: -0.1, -0.1, 1.1, 1.1
    problem = Problem(net, X, Y)

    @testset "Neurify" begin
        solver = Neurify()
        result = solve(solver, problem)
        @test result.status == :holds
    end

    @testset "ReluVal" begin
        solver = ReluVal()
        result = solve(solver, problem)
        @test result.status == :holds
    end

end