text:
using Compat: ComplexF32, ComplexF64, undef, popfirst!, IOBuffer

using FortranFiles
using Compat.Test
using Base.Iterators: product

include("test_macro_throws.jl")

using FortranFiles: RecordMarkerType, AccessMode, SequentialAccess, DirectAccess

# To test long records being split into subrecords, without needing
# test data files of several GB, define a custom RecordMarkerType with
# a short maximum subrecord length
const maxsubreclen = 2^15 - 9
const RECMRK4Bwssr = FortranFiles.WithSubrecords(maxsubreclen)
const RECMRK4Bwosr = FortranFiles.WithoutSubrecords{Int32}()
const nelem        = 73
const dareclen     = nelem*sizeof(ComplexF64)

# --- Definition of different types of record marker ---

struct RecordTypeTest
   rectyp :: Any
   desc   :: String
   target :: String
   tag    :: String
   fflags :: String
end

const rectyp_tests = [
   RecordTypeTest( DirectAccess, "fixed length (direct access)",
                   "gendatadirect", "DA", "" ),
   RecordTypeTest( RECMRK4Bwssr, "markers of 4-byte with short subrecords",
                   "gendataseq", "4Bshort", "-fmax-subrecord-length=$(maxsubreclen)" ),
   RecordTypeTest( RECMRK4B, "markers of 4-byte with default subrecords",
                   "gendataseq", "4Bdef", "" ),
   RecordTypeTest( RECMRK8B, "markers of 8-byte without subrecords",
                   "gendataseq", "8B", "-frecord-marker=8" ),
   RecordTypeTest( RECMRK4Bwosr, "markers of 4-byte without subrecords",
                   "gendataseq", "4Bdumb", "" )
   ]

# --- Definition of different byte order conversions ---

struct ByteOrderTest
   name   :: String
   tag    :: String
   fflags :: String
end

const byteorder_tests = [
   ByteOrderTest( "native",        "no",   "" ),
   ByteOrderTest( "little-endian", "LE", "-fconvert=little-endian" ),
   ByteOrderTest( "big-endian",    "BE", "-fconvert=big-endian" ),
   ]

# --- Generation of test data and Julia reading/writing code ---

function gendata(target, tag, fflags)
   cmd = `make -C codegen $target SUFFIX=_$tag XFLAGS=$fflags`
   run(cmd)
   cmd = `codegen/$(target)_$(tag).x data_$(tag).bin`
   run(cmd)
end

function genalldata(tests...)
   for test in product(tests...)
      tag = join( (t.tag for t in test), "_" )
      fflags = join( (t.fflags for t in test), " " )
      gendata(test[1].target, tag, fflags)
   end
   return nothing
end

genalldata(rectyp_tests, byteorder_tests)

include("codegen/jread.jl")
include("codegen/jfread.jl")
include("codegen/jskip.jl")
include("codegen/jwrite.jl")

# --- Auxiliary functions ---

function cmpfiles(a::String, b::String)
   cmd = `cmp $a $b`
   try
      run(cmd)
      return true
   catch
      return false
   end
end

readdata(f::FortranFile{DirectAccess})  = _readdata(f, nelem)
freaddata(f::FortranFile{DirectAccess}) = _freaddata(f, nelem)
skipdata(f::FortranFile{DirectAccess})  = _readdata(f, 1)

function _readdata(f, n::Integer)
   vari1  = read(f, rec=1,  Array{Int8}(undef, n))
   @test typeof(vari1)==Array{Int8,1}
   @test sizeof(vari1)==sizeof(Int8)*n
   vari2  = read(f, rec=2,  Array{Int16}(undef, n))
   @test typeof(vari2)==Array{Int16,1}
   @test sizeof(vari2)==sizeof(Int16)*n
   vari4  = read(f, rec=3,  Array{Int32}(undef, n))
   @test typeof(vari4)==Array{Int32,1}
   @test sizeof(vari4)==sizeof(Int32)*n
   vari8  = read(f, rec=4,  Array{Int64}(undef, n))
   @test typeof(vari8)==Array{Int64,1}
   @test sizeof(vari8)==sizeof(Int64)*n
   varr4  = read(f, rec=11, Array{Float32}(undef, n))
   @test typeof(varr4)==Array{Float32,1}
   @test sizeof(varr4)==sizeof(Float32)*n
   varr8  = read(f, rec=12, Array{Float64}(undef, n))
   @test typeof(varr8)==Array{Float64,1}
   @test sizeof(varr8)==sizeof(Float64)*n
   varc8  = read(f, rec=21, Array{ComplexF32}(undef, n))
   @test typeof(varc8)==Array{ComplexF32,1}
   @test sizeof(varc8)==sizeof(ComplexF32)*n
   varc16 = read(f, rec=22, Array{ComplexF64}(undef, n))
   @test typeof(varc16)==Array{ComplexF64,1}
   @test sizeof(varc16)==sizeof(ComplexF64)*n
   varstr = read(f, rec=30,  Array{FString{11}}(undef, n))
   @test typeof(varstr)==Array{FString{11},1}
   @test sizeof(varstr)==sizeof(FString{11})*n
   return ( vari1, vari2, vari4, vari8, varr4, varr8, varc8, varc16, varstr )
end

function _freaddata(f, n::Integer)
   @fread f rec=1 vari1::Array{Int8}(undef, n)
   @test typeof(vari1)==Array{Int8,1}
   @test sizeof(vari1)==sizeof(Int8)*n
   @fread f rec=2 vari2::Array{Int16}(undef, n)
   @test typeof(vari2)==Array{Int16,1}
   @test sizeof(vari2)==sizeof(Int16)*n
   @fread f rec=3 vari4::Array{Int32}(undef, n)
   @test typeof(vari4)==Array{Int32,1}
   @test sizeof(vari4)==sizeof(Int32)*n
   @fread f rec=4 vari8::Array{Int64}(undef, n)
   @test typeof(vari8)==Array{Int64,1}
   @test sizeof(vari8)==sizeof(Int64)*n
   @fread f rec=11 varr4::Array{Float32}(undef, n)
   @test typeof(varr4)==Array{Float32,1}
   @test sizeof(varr4)==sizeof(Float32)*n
   @fread f rec=12 varr8::Array{Float64}(undef, n)
   @test typeof(varr8)==Array{Float64,1}
   @test sizeof(varr8)==sizeof(Float64)*n
   @fread f rec=21 varc8::Array{ComplexF32}(undef, n)
   @test typeof(varc8)==Array{ComplexF32,1}
   @test sizeof(varc8)==sizeof(ComplexF32)*n
   @fread f rec=22 varc16::Array{ComplexF64}(undef, n)
   @test typeof(varc16)==Array{ComplexF64,1}
   @test sizeof(varc16)==sizeof(ComplexF64)*n
   @fread f rec=30 varstr::Array{FString{11}}(undef, n)
   @test typeof(varstr)==Array{FString{11},1}
   @test sizeof(varstr)==sizeof(FString{11})*n
   return ( vari1, vari2, vari4, vari8, varr4, varr8, varc8, varc16, varstr )
end

function writedata(f::FortranFile{DirectAccess}, data)
   vari1, vari2, vari4, vari8, varr4, varr8, varc8, varc16, varstr = data
   nwritten = write(f, rec=1,  vari1)
   @test nwritten==sizeof(vari1)
   nwritten = write(f, rec=2,  vari2)
   @test nwritten==sizeof(vari2)
   nwritten = write(f, rec=3,  vari4)
   @test nwritten==sizeof(vari4)
   nwritten = write(f, rec=4,  vari8)
   @test nwritten==sizeof(vari8)
   nwritten = write(f, rec=11, varr4)
   @test nwritten==sizeof(varr4)
   nwritten = write(f, rec=12, varr8)
   @test nwritten==sizeof(varr8)
   nwritten = write(f, rec=21, varc8)
   @test nwritten==sizeof(varc8)
   nwritten = write(f, rec=22, varc16)
   @test nwritten==sizeof(varc16)
   nwritten = write(f, rec=30, varstr)
   @test nwritten==sizeof(varstr)
   # going back
   nwritten = write(f, rec=16, varc16)
   @test nwritten==sizeof(varc16)
   # overwriting
   nwritten = write(f, rec=16, varstr)
   @test nwritten==sizeof(varstr)
end

# --- TESTING STARTS HERE ---

@testset "Exceptions" begin
   @testset "Constructor" begin
      open("/dev/null", "w") do io
         # sequential access file with fixed-length records
         @test_throws FortranFilesError FortranFile(io, recl=80)
         # direct access file without record length
         @test_throws FortranFilesError FortranFile(io, access="direct")
         # unknown access mode
         @test_throws FortranFilesError FortranFile(io, access="stream")
         # unknown byte-order conversion
         @test_throws FortranFilesError FortranFile(io, convert="host")
      end
   end

   @testset "Writing" begin
      open("/dev/null", "w") do io
         # direct access files
         fdir = FortranFile(io, access="direct", recl=80)
         # no record number
         @test_throws FortranFilesError write(fdir, zeros(10))
         # write too much
         @test_throws FortranFilesError write(fdir, rec=1, zeros(11))

         # sequential access files
         for recmrk in (RECMRK4B,RECMRK4Bwosr)
            fseq = FortranFile(io, marker=recmrk)
            # with record number
            @test_throws MethodError write(fseq, rec=1, zeros(10))
            # inhomogeneous array, see PR#4
            inhomA = Integer[1, big(2)]
            @test_throws FortranFilesError write(fseq, inhomA)
         end
      end
   end

   @testset "Reading" begin
      let A=ones(10), B=2A
         # direct access files
         data = UInt8[]
         io = IOBuffer(data, write=true)
         fdir = FortranFile(io, access="direct", recl=80)
         write(fdir, rec=1, A)
         write(fdir, rec=2, B)
         close(fdir)
         io = IOBuffer(data)
         fdir = FortranFile(io, access="direct", recl=80)
         # no record number
         @test_throws FortranFilesError AA = read(fdir, (Float64,10))
         # macro version
         @test_throws FortranFilesError @fread fdir AA::(Float64,10)
         # read too much
         @test_throws FortranFilesError AA = read(fdir, rec=1, (Float64,11))
         close(fdir)

         # sequential access files
         for recmrk in (RECMRK4B,RECMRK4Bwosr)
            data = UInt8[]
            io = IOBuffer(data, write=true)
            fseq = FortranFile(io, marker=recmrk)
            write(fseq, A)
            write(fseq, B)
            close(fseq)
            io = IOBuffer(data)
            fseq = FortranFile(io, marker=recmrk)
            # with record number
            @test_throws MethodError AA = read(fseq, rec=1, (Float64,10))
            # macro version
            @test_throws FortranFilesError @fread fseq rec=1 AA::(Float64,10)
            # read too much
            @test_throws FortranFilesError AA = read(fseq, (Float64,11))
            close(fseq)
            # garbage data: wrong trailing record marker
            data = UInt8[]
            io = IOBuffer(data, write=true)
            write(io, UInt32(80), A, UInt32(40))
            close(io)
            io = IOBuffer(data)
            fseq = FortranFile(io, marker=recmrk)
            @test_throws FortranFilesError AA = read(fseq, (Float64,10))
            close(fseq)
         end
      end
   end

   @testset "Macro syntax" begin
      @test_macro_throws ArgumentError eval(:(@fread f foo=1))
      @test_macro_throws ArgumentError eval(:(@fread f 10))
      @test_macro_throws ArgumentError eval(:(@fread f (Int,Int)))
   end
end

@testset "Strings" begin
   jstr = "Hello World!"
   N    = length(jstr)
   fstr = FString(N, jstr)
   @test typeof(fstr) == FString{N}
   @test trimstring(fstr) == jstr
   @test trimlen(fstr) == N
   @test trim(fstr) == fstr
   @test String(fstr) == jstr
   fstr80 = FString(80, jstr)
   @test typeof(fstr80) == FString{80}
   @test trimstring(fstr80) == jstr
   @test trimlen(fstr80) == N
   @test trim(fstr80) == fstr
   @test String(fstr80) == rpad(jstr, 80)
   fstr8 = FString(8, jstr)
   @test typeof(fstr8) == FString{8}
   @test trimstring(fstr8) == jstr[1:8]
   @test trimlen(fstr8) == 8
   @test trim(fstr8) == fstr8
   @test String(fstr8) == jstr[1:8]
   jstr = "Hällo Wörld!"
   @test_throws InexactError FString(80, jstr)
end

@testset "Testing records with $(rectest.desc), $(botest.name) byte order" for rectest in rectyp_tests, botest in byteorder_tests

   local infile, outfile, data, data2
   tag = "$(rectest.tag)_$(botest.tag)"
   infilename  = "data_$(tag).bin"
   outfilename = "chck_$(tag).bin"

   @testset "Opening files" begin
      if isa(rectest.rectyp, RecordMarkerType)
         infile  = FortranFile(infilename, "r", marker = rectest.rectyp, convert = botest.name)
         @test infile.acctyp.recmrktyp == rectest.rectyp
         outfile = FortranFile(outfilename, "w", marker = rectest.rectyp, convert = botest.name)
         @test outfile.acctyp.recmrktyp == rectest.rectyp
      else
         infile  = FortranFile(infilename, "r", access="direct", recl=dareclen, convert = botest.name)
         @test infile.acctyp.reclen == dareclen
         outfile = FortranFile(outfilename, "w", access="direct", recl=dareclen, convert = botest.name)
         @test infile.acctyp.reclen == dareclen
      end
   end

   @time @testset "Reading data" begin
      data = readdata(infile)
   end

   @time @testset "@fread'ing data " begin
      rewind(infile)
      data2 = freaddata(infile)
      @test data == data2
   end

   @time @testset "Reading data with skipping" begin
      rewind(infile)
      skipdata(infile)
      close(infile)
   end

   @time @testset "Writing data" begin
      writedata(outfile, data)
      close(outfile)
   end

   @testset "Verifying data" begin
      @test cmpfiles(infilename, outfilename)
   end

end;


text:
## Abstract Markov Chain
## =====================

AbstractChain{T<:Domain} = Vector{T}

"Get point samples out of a abstract Markov chain"
point_sample(chain::AbstractChain) = [rand(state) for state in chain]

ext: jl


text:
using Publish
using ParameterSchedulers
using Artifacts, LazyArtifacts

# override default theme
Publish.Themes.default() = artifact"flux-theme"

p = Publish.Project(ParameterSchedulers)

# serve documentation
serve(ParameterSchedulers)


text:
using AA228FinalProject
using POMDPGifs
using Random
using ParticleFilters
using POMDPPolicies
using ARDESPOT
using POMDPs

rng = MersenneTwister(2)
speed = 2.0
as = vec([RoombaAct(v, om) for v in (0.0, speed), om in (-1.0, 0.0, 1.0)])
m = RoombaPOMDP(sensor=Bumper(), mdp=RoombaMDP(config=1, aspace=as, contact_pen=-0.1));

default = FunctionPolicy(x->[speed, 0.0])
bounds = IndependentBounds(DefaultPolicyLB(default), 10.0, check_terminal=true)
solver = DESPOTSolver(K=20, T_max=1.0, bounds=bounds, rng=rng)
planner = solve(solver, m)

spf = SimpleParticleFilter(m, BumperResampler(5000), rng=rng)
filter = RoombaParticleFilter(spf, 2.0, 0.5);

makegif(m, planner, filter, filename="out.gif", rng=rng, max_steps=100, show_progress=true)


text:
function faceconnectivity2(t,f2t,dim,elemtype,porder)

nf = size(f2t,2);

#display([porder dim elemtype])
philocvl,~,~,plocfc,perm = localbasis(porder,dim,elemtype);
permind = permindex(plocfc,dim,elemtype);
face = getelemface(dim,elemtype);

npf = size(perm,1);
npe = size(philocvl,1);

facecon = zeros(Int,npf,2,nf);
if dim<=2
    for i = 1:nf
        e1 = f2t[1,i];
        l1 = f2t[2,i];
        e2 = f2t[3,i];
        l2 = f2t[4,i];

        facecon[:,1,i] = (e1-1)*npe .+ perm[:,l1];

        if e2>0 # face i is an interior face
            facecon[:,2,i] = (e2-1)*npe .+ perm[permind,l2];
        else
            facecon[:,2,i] = facecon[:,1,i];
        end
    end
else
    for i = 1:nf
        e1 = f2t[1,i];
        l1 = f2t[2,i];
        e2 = f2t[3,i];
        l2 = f2t[4,i];

        facecon[:,1,i] = (e1-1)*npe .+ perm[:,l1];
        if e2>0 # face i is an interior face
            f1 = t[face[:,l1],e1];
            f2 = t[face[:,l2],e2];
            if elemtype==0
                if (f1[1]==f2[1]) && (f1[2]==f2[3]) && (f1[3]==f2[2])
                    k = 1;
                elseif (f1[1]==f2[2]) && (f1[2]==f2[1]) && (f1[3]==f2[3])
                    k = 2;
                elseif (f1[1]==f2[3]) && (f1[2]==f2[2]) && (f1[3]==f2[1])
                    k = 3;
                else
                    error("Mesh connectivity is wrong");
                end
            else
                if (f1[1]==f2[1]) && (f1[2]==f2[4]) && (f1[3]==f2[3]) && (f1[4]==f2[2])
                    k = 1;
                elseif (f1[1]==f2[4]) && (f1[2]==f2[3]) && (f1[3]==f2[2]) && (f1[4]==f2[1])
                    k = 2;
                elseif (f1[1]==f2[3]) && (f1[2]==f2[2]) && (f1[3]==f2[1]) && (f1[4]==f2[4])
                    k = 3;
                elseif (f1[1]==f2[2]) && (f1[2]==f2[1]) && (f1[3]==f2[4]) && (f1[4]==f2[3])
                    k = 4;
                else
                    error("Mesh connectivity is wrong");
                end
            end
            facecon[:,2,i] = (e2-1)*npe .+ perm[permind[:,k],l2];
        else
            facecon[:,2,i] = facecon[:,1,i];
        end
    end
end

return facecon

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
isgzipped( filename::String ) = hasextension( filename, "gz" )

function hasextension( filename::String, ext::String )
   restr = "\.$ext\$"
   re = Base.match(Regex(restr), filename)
   return re == nothing ? false : true
end

plug_write{S <: AbstractString}( io, str::S; plug::Char='\t' ) = (write( io, str ); write( io, plug  ))
plug_write( io, str::Char; plug::Char='\t' ) = (write( io, str ); write( io, plug  ))

tab_write{S <: AbstractString}( io, str::S ) = plug_write( io, str, plug='\t' )
tab_write( io, str::Char ) = plug_write( io, str, plug='\t' )

end_write{S <: AbstractString}( io, str::S ) = plug_write( io, str, plug='\n' )
end_write( io, str::Char ) = plug_write( io, str, plug='\n' )

function tab_write( io, dat::Tuple )
   for i in 1:length(dat)
      write( io, string(dat[i]) )
      if i < length(dat)
         write( io, '\t' )
      end
   end
   write( io, '\n' )
end

default_colors( n::Int ) = [Scale.color_continuous().f(p) for p in linspace(0, 1, n)]
default_colors( n::Int, opacity::Float64 ) = [(t = convert(Gadfly.RGBA, c);
                                              return Gadfly.RGBA(t.r, t.g, t.b, opacity))
                                              for c in default_colors(n)]

const DEF_COLORS = default_colors( 100, 1.0 )

random_rgb() = (rand(0:255), rand(0:255), rand(0:255))
random_rgb_str() = join(random_rgb(), ",")

rgba_to_int( col::Gadfly.RGBA ) = Int(floor(col.r*255)),Int(floor(col.g*255)),Int(floor(col.b*255))

function weighted_rgb( dg::Float64, maxenergy::Float64 )
   if dg < 0.0
      ind = min(Int(floor((abs(dg) / abs(maxenergy)) * 100)), 100)
      col = DEF_COLORS[ind]
   else
      col = DEF_COLORS[1]
   end
   join(rgba_to_int(col), ",")
end

function randcode( arr=[collect('0':'9'); collect('A':'Z')], len=7 )
   shuffle!(arr)
   join(arr[1:len])
end

function writebed{T}( io, dup::DuplexInterval{T}, name::String; maxdistance::Int=100_000, maxenergy::Float64=-35.0 )
   const dist = dup.last.first - dup.first.last
   if dup.first.seqname == dup.last.seqname &&
      dist <= maxdistance
      tab_write( io, string(dup.first.seqname) )
      tab_write( io, string(dup.first.first-1) )
      tab_write( io, string(dup.last.last) )
      tab_write( io, "EvoDuplexes: " * name )
      tab_write( io, string(min(1000, energy(dup.duplex) * -250)) )
      tab_write( io, dup.first.strand == STRAND_POS ? '+' : '-' )
      tab_write( io, string(dup.first.first-1) )
      tab_write( io, string(dup.last.last) )
      tab_write( io, weighted_rgb(energy(dup.duplex), maxenergy) )
      tab_write( io, "2" )
      tab_write( io, string(length(dup.first.first:dup.first.last)) * "," * string(length(dup.last.first:dup.last.last)) )
      end_write( io, string(0) * "," * string(dup.last.first-dup.first.first) )
   end
end

function writebed{T}( io, v::Vector{DuplexInterval{T}}, name::String; maxdistance::Int=100_000, maxenergy::Float64=-35.0 )
   for i in v
      writebed( io, i, name, maxdistance=maxdistance, maxenergy=maxenergy )
   end
end


text:
#EXERCISE MADE BY ERIK MARTÍN GARZÓN & DAVID BOCES OBIS
# the packages we need
using Gtk, Graphics, Logging, Printf

include("affin_transformation.jl")
#File opening and reading matrix and vor together
A=readdlm("circle.txt", Float64)

#Reading and saving matrix size
s=size(A)
s1=s[1] #rows
s2=s[2] #columns

# the main window
win = GtkWindow("SO(3)")

function init_canvas(h,w)
    # create the drawing canvas
    c = @GtkCanvas(h,w)

    # create the initial drawing inside a try/catch loop
    @guarded draw(c) do widget
        #  draw the background with the canvas drawing context
        # the code for this comes later
        draw_the_canvas(c)
    end

    show(c)
    return c
end

# make the canvas
the_canvas = init_canvas(1500,750)

# an array to store the entry boxes
entry_list = []

# an array of labels that we use to display normalized inputs,
# and which also gets modified from the callback
normalized_labels = []

function find_by_name(list, name)
    for item in list
        if get_gtk_property(item, :name, String) == name
            return item
        end
    end
    @warn name, "not found in list"
    @warn "available names are"
    for item in list
        @warn get_gtk_property(item, :name, String)
    end
end

function bold_label(label_string)
    label = GtkLabel("")
    GAccessor.markup(label, """<b>""" * label_string * """</b>""")
    return label
end

# Now put everything into the window,
# including the canvas
function init_window(win, canvas)
    # make a box for the drawing canvas
    canvas_box = GtkBox(:v)
    push!(canvas_box, bold_label("EJERCICIO 3.1:                                                                                EJERCICIO 3.2:"))
    push!(canvas_box, bold_label("Proyecta los puntos del archivo al plano focal de la camara y dibuja el resultado.                                           Haz un dibujo 3d de toda la escena, con los dos sistemas ortonormales de coordenadas y
                                                                                los puntos de la circunferencia."))
    push!(canvas_box, canvas)

    # make a containing box that will stack the widgets and the canvas side by side
    global_box = GtkBox(:h)
    push!(global_box, GtkLabel("   ")) # a very basic separator
    push!(global_box, canvas_box)

    # put it all inside the window
    push!(win, global_box)
end

# the background drawing
function draw_the_canvas(canvas)
    h   = height(canvas)
    w   =  width(canvas)
    ctx =  getgc(canvas)

    # clear the canvas
    rectangle(ctx, 750, 0, w, h)
    set_source_rgb(ctx, 1, 1, 1)
    fill(ctx)

    # Paint red rectangle
    # rectangle(ctx, 3w/4, 0, w/4, h/4)
    # set_source_rgb(ctx, 1, 0, 0)
    # fill(ctx)

    # Paint blue rectangle
    # rectangle(ctx, 0, 3h/4, w/4, h/4)
    # set_source_rgb(ctx, 0, 0, 1)
    # fill(ctx)

    # Paint an example line
    # set_line_width(ctx, 5)
    # set_source_rgb(ctx, 1, 0, 1)
    # move_to(ctx, 10, 10)
    # line_to(ctx, 100, 50)
    # stroke(ctx)

    # ------------------------------------------------------------------ DATA
    originX = 850
    originY = 300

    #AXIS R3
    X = [1;0;0]
    Y = [0;1;0]
    Z = [0;0;1]

    #WE PROJECT THE ROTATED AXIS IN R3 TO R2
    X2 = to_2d(X)
    Y2 = to_2d(Y)
    Z2 = to_2d(Z)

    cam = [1;6;1]

    ay = 90
    az = -20
    focus = 3

    #Setting camera coords
    camCoord2d = to_2d(cam)*100

    #Camera axis rotations
    camXRY = rotate_phi_z(ay, Y, X)
    camYRY = rotate_phi_z(ay, Y, Y)
    camZRY = rotate_phi_z(ay, Y, Z)

    camXRZ = rotate_phi_z(az, camZRY, camXRY)
    camYRZ = rotate_phi_z(az, camZRY, camYRY)
    camZRZ = rotate_phi_z(az, camZRY, camZRY)

    camX2d = to_2d(camXRZ)
    camY2d = to_2d(camYRZ)
    camZ2d = to_2d(camZRZ)

    #Setting up variables we will need for calculations
    Ry = [cosd(ay) 0 sind(ay);0 1 0;-sind(ay) 0 cosd(ay)]
    Rz = [cosd(az) -sind(az) 0; sind(az) cosd(az) 0;0 0 1]
    C = Rz * Ry
    D = hcat(C, cam)
    B = vcat(D, [0 0 0 1])
    Binv = B^-1

    # for i=1:s1
    #     println(to_2d(focal_camera_points_in_space(B, focus, to_2d(focal_camera(Binv, focus, A[i,:])))[1:3,1]))
    #     println()
    #     println(focal_camera(Binv, focus, A[i,:]))
    # end

    # ----------------------------------------------------------- EXERCISE 3.1
    for i=1:s1
        circle(ctx, 375 + focal_camera(Binv, focus, A[i,:])[1]*100, 200 + focal_camera(Binv, focus, A[i,:])[2]*100, 1)
        set_source_rgb(ctx, 1, 0, 0)
        fill(ctx)
    end

    # ----------------------------------------------------------- EXERCISE 3.2
    # DRAW CANON BASIS
    #DRAW X AXIS CANON (RED)
    set_line_width(ctx, 2)
    set_source_rgb(ctx, 1, 0, 0)
    move_to(ctx, originX, originY)
    line_to(ctx, originX + X2[1]*100, originY - X2[2]*100)
    stroke(ctx)

    #circle at the end of the axis
    circle(ctx, originX + X2[1]*100, originY - X2[2]*100, 5)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)

    #DRAW Y AXIS CANON (GREEN)
    set_line_width(ctx, 2)
    set_source_rgb(ctx, 0, 1, 0)
    move_to(ctx, originX, originY)
    line_to(ctx, originX + Y2[1]*100, originY - Y2[2]*100)
    stroke(ctx)

    #circle at the end of the axis
    circle(ctx, originX + Y2[1]*100, originY - Y2[2]*100, 5)
    set_source_rgb(ctx, 0, 1, 0)
    fill(ctx)

    #DRAW Z AXIS CANON (BLUE)
    set_line_width(ctx, 2)
    set_source_rgb(ctx, 0, 0, 1)
    move_to(ctx, originX, originY)
    line_to(ctx, originX + Z2[1]*100, originY - Z2[2]*100)
    stroke(ctx)

    #circles at the end of the axis
    circle(ctx, originX + Z2[1]*100, originY - Z2[2]*100, 5)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)

    # DRAW CIRCLE POINTS IN SPACE
    for i=1:s1
        circle(ctx, originX + to_2d(A[i,:])[1]*100, originY - to_2d(A[i,:])[2]*100, 2)
        set_source_rgb(ctx, 0, 0, 1)
        fill(ctx)
    end

    # DRAW FOCAL CAMERA POINTS IN SPACE
    for i=1:s1
        circle(ctx, originX + to_2d(focal_camera_points_in_space(B, focus, focal_camera(Binv, focus, A[i,:])[1:2,:])[1:3,:])[1]*100,
         originY - to_2d(focal_camera_points_in_space(B, focus, focal_camera(Binv, focus, A[i,:])[1:2,:])[1:3,:])[2]*100, 1)
        set_source_rgb(ctx, 1, 0, 0)
        fill(ctx)
    end

    # WE DRAW THE AXIS' CAMERA
    #DRAW X AXIS CAMERA (RED)
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 1, 0, 0)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + camCoord2d[1] + camX2d[1]*50, originY - camCoord2d[2] - camX2d[2]*50)
    stroke(ctx)

    #circle at the end of the axis
    circle(ctx, originX + camCoord2d[1] + camX2d[1]*50, originY - camCoord2d[2] - camX2d[2]*50, 2)
    set_source_rgb(ctx, 1, 0, 0)
    fill(ctx)

    #DRAW Y AXIS CAMERA (GREEN)
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 1, 0)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + camCoord2d[1] + camY2d[1]*50, originY - camCoord2d[2] - camY2d[2]*50)
    stroke(ctx)

    #circle at the end of the axis
    circle(ctx, originX + camCoord2d[1] + camY2d[1]*50, originY - camCoord2d[2] - camY2d[2]*50, 2)
    set_source_rgb(ctx, 0, 1, 0)
    fill(ctx)

    #DRAW Z AXIS CANON (BLUE)
    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 0, 1)
    move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
    line_to(ctx, originX + camCoord2d[1] + camZ2d[1]*50, originY - camCoord2d[2] - camZ2d[2]*50)
    stroke(ctx)

    #circles at the end of the axis
    circle(ctx, originX + camCoord2d[1] + camZ2d[1]*50, originY - camCoord2d[2] - camZ2d[2]*50, 2)
    set_source_rgb(ctx, 0, 0, 1)
    fill(ctx)

    # LINES FORM EACH CIRCLE POINT TO ORIGIN CAMERA
    for i=1:s1
        set_line_width(ctx, 0.2)
        set_source_rgb(ctx, 0, 0.5, 1)
        move_to(ctx, originX + camCoord2d[1], originY - camCoord2d[2])
        line_to(ctx, originX + to_2d(A[i,:])[1]*100, originY - to_2d(A[i,:])[2]*100)
        stroke(ctx)
    end

    stroke(ctx)

end

# -------- initialize everything ---------

# prepare and show the initial widgets
init_window(win, the_canvas)
showall(win)


text:
using AutomotiveDrivingModels
using AutoDrivers.GaussianMixtureRegressionDrivers

########################################
#            CollisionRate             #
########################################

type CollisionRate <: TraceMetricExtractor
    running_sum::Float64
    n_obs::Int

    CollisionRate(running_sum::Float64=0.0, n_obs::Int=0) = new(running_sum, n_obs)
end
Base.Symbol(::CollisionRate) = :collisionrate
get_score(m::CollisionRate) = m.running_sum / m.n_obs
function reset!(metric::CollisionRate)
    metric.running_sum = 0.0
    metric.n_obs = 0
    metric
end
function extract_collision_rate(rec::SceneRecord, roadway::Roadway, egoid::Int)
    collided = 0.0
    for pastframe in 1 - length(rec) : 0
        vehicle_index = get_index_of_first_vehicle_with_id(rec, egoid, pastframe)
        collide = convert(Float64, get(IS_COLLIDING, rec, roadway, vehicle_index, pastframe))
        collided = min(collide + collided, 1.0) # only count if first collision
    end
    collided
end
function extract!(
    metric::CollisionRate,
    rec_orig::SceneRecord, # the records are exactly as long as the simulation (ie, contain everything)
    rec_sim::SceneRecord,
    roadway::Roadway,
    egoid::Int;
    real::Bool=false,
    )
    if real
        metric.running_sum += extract_collision_rate(rec_orig, roadway, egoid)
    else
        metric.running_sum += extract_collision_rate(rec_sim, roadway, egoid)
    end
    metric.n_obs += 1

    metric
end

########################################
#            OffRoadRate               #
########################################

type OffRoadRate <: TraceMetricExtractor
    running_sum::Float64
    n_obs::Int

    OffRoadRate(running_sum::Float64=0.0, n_obs::Int=0) = new(running_sum, n_obs)
end
Base.Symbol(::OffRoadRate) = :offroadrate
get_score(m::OffRoadRate) = m.running_sum / m.n_obs
function reset!(metric::OffRoadRate)
    metric.running_sum = 0.0
    metric.n_obs = 0
    metric
end
function extract_off_road_rate(rec::SceneRecord, roadway::Roadway, egoid::Int)
    num_off_road = 0.0
    for pastframe in 1 - length(rec) : 0
        vehicle_index = get_index_of_first_vehicle_with_id(rec, egoid, pastframe)
        off_road = (convert(Float64, get(MARKERDIST_LEFT, rec, roadway, vehicle_index, pastframe)) < -1.0) +
                    (convert(Float64, get(MARKERDIST_RIGHT, rec, roadway, vehicle_index, pastframe)) < -1.0)
        num_off_road += off_road
    end
    num_off_road
end
function extract!(
    metric::OffRoadRate,
    rec_orig::SceneRecord, # the records are exactly as long as the simulation (ie, contain everything)
    rec_sim::SceneRecord,
    roadway::Roadway,
    egoid::Int;
    real::Bool=false,
    )
    if real
        metric.running_sum += extract_off_road_rate(rec_orig, roadway, egoid)
    else
        metric.running_sum += extract_off_road_rate(rec_sim, roadway, egoid)
    end
    metric.n_obs += 1

    metric
end

########################################
#            LaneChangeRate            #
########################################

type LaneChangeRate <: TraceMetricExtractor
    running_sum::Float64
    n_obs::Int

    LaneChangeRate(running_sum::Float64=0.0, n_obs::Int=0) = new(running_sum, n_obs)
end
Base.Symbol(::LaneChangeRate) = :lanechangerate
get_score(m::LaneChangeRate) = m.running_sum / m.n_obs
function reset!(metric::LaneChangeRate)
    metric.running_sum = 0.0
    metric.n_obs = 0
    metric
end
function extract_lane_change_rate(rec::SceneRecord, roadway::Roadway, egoid::Int)
    num_lane_change = 0.0
    vehicle_index = get_index_of_first_vehicle_with_id(rec, egoid, 1 - length(rec))
    veh_ego = rec[vehicle_index, 1 - length(rec)]
    prev_lane = roadway[veh_ego.state.posF.roadind.tag] # Find initial lane
    for pastframe in 2 - length(rec) : 0
        vehicle_index = get_index_of_first_vehicle_with_id(rec, egoid, pastframe)
        veh_ego = rec[vehicle_index, pastframe]
        lane = roadway[veh_ego.state.posF.roadind.tag]
        num_lane_change += (lane != prev_lane)
        prev_lane = lane
    end
    num_lane_change
end
function extract!(
    metric::LaneChangeRate,
    rec_orig::SceneRecord, # the records are exactly as long as the simulation (ie, contain everything)
    rec_sim::SceneRecord,
    roadway::Roadway,
    egoid::Int;
    real::Bool=false,
    )
    if real
        metric.running_sum += extract_lane_change_rate(rec_orig, roadway, egoid)
    else
        metric.running_sum += extract_lane_change_rate(rec_sim, roadway, egoid)
    end
    metric.n_obs += 1

    metric
end

########################################
#            HardBrakeRate             #
########################################

type HardBrakeRate <: TraceMetricExtractor
    running_sum::Float64
    n_obs::Int

    HardBrakeRate(running_sum::Float64=0.0, n_obs::Int=0) = new(running_sum, n_obs)
end
Base.Symbol(::HardBrakeRate) = :hardbrakerate
get_score(m::HardBrakeRate) = m.running_sum / m.n_obs
function reset!(metric::HardBrakeRate)
    metric.running_sum = 0.0
    metric.n_obs = 0
    metric
end
function extract_hard_brake_rate(rec::SceneRecord, roadway::Roadway, egoid::Int)
    num_hard_brake = 0.0
    for pastframe in 2 - length(rec) : 0
        vehicle_index = get_index_of_first_vehicle_with_id(rec, egoid, pastframe)
        acc_ego = convert(Float64, get(ACC, rec, roadway, vehicle_index, pastframe))
        num_hard_brake += (acc_ego < -3.0)
    end
    num_hard_brake
end
function extract!(
    metric::HardBrakeRate,
    rec_orig::SceneRecord, # the records are exactly as long as the simulation (ie, contain everything)
    rec_sim::SceneRecord,
    roadway::Roadway,
    egoid::Int;
    real::Bool=false,
    )
    if real
        metric.running_sum += extract_hard_brake_rate(rec_orig, roadway, egoid)
    else
        metric.running_sum += extract_hard_brake_rate(rec_sim, roadway, egoid)
    end
    metric.n_obs += 1

    metric
end

########################################
#        RWSE Global Position          #
########################################

type RWSEPosG <: TraceMetricExtractor
    horizon::Float64 # [s]
    running_sum::Float64
    n_obs::Int
end
RWSEPosG(horizon::Float64) = RWSEPosG(horizon, 0.0, 0)

Base.Symbol(m::RWSEPosG) = Symbol(@sprintf("RWSE_posG_%d_%02d", floor(Int, m.horizon), floor(Int, 100*rem(m.horizon, 1.0))))
get_score(m::RWSEPosG) = sqrt(m.running_sum / m.n_obs)
function reset!(metric::RWSEPosG)
    metric.running_sum = 0.0
    metric.n_obs = 0
    metric
end
function extract!(
    metric::RWSEPosG,
    rec_orig::SceneRecord, # the records are exactly as long as the simulation (ie, contain everything)
    rec_sim::SceneRecord,
    roadway::Roadway,
    egoid::Int,
    )

    # TODO: how to handle missing values???

    pastframe = 1-length(rec_orig) + clamp(round(Int, metric.horizon/rec_orig.timestep), 0, length(rec_orig)-1)

    # pull true value
    vehicle_index = get_index_of_first_vehicle_with_id(rec_orig, egoid, pastframe)
    veh_ego = rec_orig[vehicle_index, pastframe]
    x_true = veh_ego.state.posG.x
    y_true = veh_ego.state.posG.y

    # pull sim value
    vehicle_index = get_index_of_first_vehicle_with_id(rec_sim, egoid, pastframe)
    veh_ego = rec_sim[vehicle_index, pastframe]
    x_montecarlo = veh_ego.state.posG.x
    y_montecarlo = veh_ego.state.posG.y

    d = (x_true - x_montecarlo)^2 + (y_true - y_montecarlo)^2

    metric.running_sum += d
    metric.n_obs += 1

    metric
end

# Check for outliers, which seem to  occur with controller rollouts
function check_outliers(rec_orig::SceneRecord, rec_sim::SceneRecord, egoid::Int)
    pastframe = 1-length(rec_orig) + clamp(round(Int, 5.0/rec_orig.timestep), 0, length(rec_orig)-1)

    # pull true value
    vehicle_index = get_index_of_first_vehicle_with_id(rec_orig, egoid, pastframe)
    veh_ego = rec_orig[vehicle_index, pastframe]
    x_true = veh_ego.state.posG.x
    y_true = veh_ego.state.posG.y

    # pull sim value
    vehicle_index = get_index_of_first_vehicle_with_id(rec_sim, egoid, pastframe)
    veh_ego = rec_sim[vehicle_index, pastframe]
    x_montecarlo = veh_ego.state.posG.x
    y_montecarlo = veh_ego.state.posG.y

    d = sqrt((x_true - x_montecarlo)^2 + (y_true - y_montecarlo)^2)
    return (d >= 50.0)
end

function rollout!(
    rec::SceneRecord,
    model::DriverModel,
    egoid::Int,
    trajdata::Trajdata,
    time_start::Float64,
    time_end::Float64,
    simparams::Auto2D.SimParams;
    prime_history::Int = 0,
    gru_type::Bool=true,
    )

    # Initialize values
    Δt = rec.timestep
    simstate = simparams.simstates[1]

    # clear rec and make first observations
    #println("fieldnames(model.net)=", fieldnames(model.net))
    #if Symbol("gru") in fieldnames(model.net)
    if gru_type
        #println("rollout::gru!!")
        model.net[:gru].h_prev = zeros(length(model.net[:gru].h_prev))
    end
    empty!(simstate.rec)
    empty!(rec)
    reset_hidden_state!(model)
    for i in 1 : prime_history
        t = time_start - Δt * (prime_history - i + 1)
        update!(simstate.rec, get!(simstate.scene, trajdata, t))
        AutomotiveDrivingModels.observe!(model, simparams, simstate.scene, trajdata.roadway, simstate.egoid)
    end

    # empty playback reactive
    empty!(simstate.playback_reactive_active_vehicle_ids)

    # Run start time step
    t = time_start
    update!(simstate.rec, get!(simstate.scene, trajdata, t))
    update!(model.rec, simstate.scene)
    AutomotiveDrivingModels.observe!(model, simparams, simstate.scene, trajdata.roadway, simstate.egoid)

    # run simulation
    while t < time_end

        # Find action and step forward

        ego_action = rand(model)
        a = clamp(ego_action.a, -5.0, 3.0)
        ω = clamp(ego_action.ω, -0.1, 0.1)
        Auto2D.step(simparams, [a, ω])

        # update record
        AutomotiveDrivingModels.observe!(model, simparams, simstate.scene, trajdata.roadway, simstate.egoid)
        update!(rec, simstate.scene)

        # update time
        t += Δt
    end
    rec
end

function simulate!(
    rec::SceneRecord,
    model::DriverModel,
    egoid::Int,
    trajdata::Trajdata,
    time_start::Float64,
    time_end::Float64,
    simparams::Auto2D.SimParams;
    prime_history::Int = 0,
    )

    # Initialize values
    Δt = rec.timestep
    simstate = simparams.simstates[1]

    # clear rec and make first observations
    empty!(simstate.rec)
    empty!(rec)
    reset_hidden_state!(model)
    for i in 1 : prime_history
        t = time_start - Δt * (prime_history - i + 1)
        update!(simstate.rec, get!(simstate.scene, trajdata, t))
        AutomotiveDrivingModels.observe!(model, simstate.scene, trajdata.roadway, egoid)
    end

    # empty playback reactive
    empty!(simstate.playback_reactive_active_vehicle_ids)

    # Run start time step
    t = time_start
    update!(simstate.rec, get!(simstate.scene, trajdata, t))

    # If controller, set desired speed to be inital speed
    veh_index = get_index_of_first_vehicle_with_id(simstate.scene, egoid)
    if simparams.ego_action_type == LatLonAccel
        set_desired_speed!(model, simstate.scene[veh_index].state.v)
    end
    AutomotiveDrivingModels.observe!(model, simstate.scene, trajdata.roadway, egoid)

    # run simulation
    while t < time_end

        # Find action and step forward
        ego_action = rand(model)
        veh_index = get_index_of_first_vehicle_with_id(simstate.scene, egoid)
        if simstate.scene[veh_index].state.v < 0.0 && (simparams.ego_action_type == LatLonAccel)
            ego_action = LatLonAccel(ego_action.a_lat, max(ego_action.a_lon, 0.0))
        end
        if simparams.ego_action_type == AccelTurnrate
            Auto2D.step(simparams, [ego_action.a, ego_action.ω])
        else
            Auto2D.step(simparams, [ego_action.a_lat, a_lon])
        end
        AutomotiveDrivingModels.observe!(model, simstate.scene, trajdata.roadway, egoid)

        # update record
        update!(rec, simstate.scene)

        # update time
        t += Δt
    end
    rec
end

# Initialize values for simstate
function reset_simstate!(simstate::Auto2D.SimState, seg::TrajdataSegment)
    simstate.egoid = seg.egoid
    simstate.trajdata_index = seg.trajdata_index
    simstate.frame = seg.frame_lo
end

########################################
#        Calculate Metric Scores       #
########################################

function calc_metrics!(
    metrics_df::DataFrame,
    model::DriverModel,
    metrics::Vector{TraceMetricExtractor},
    simparams::Auto2D.SimParams,
    foldset_seg_test::FoldSet; # should match the test segments in evaldata
    n_simulations_per_trace::Int = 10,
    row::Int = foldset_seg_test.fold, # row in metrics_df to write to
    prime_history::Int = 0,
    calc_logl::Bool = true,
    gru_type::Bool = true,
    max_loop::Int = 1000,
    )
    # reset metrics
    for metric in metrics
        try
            reset!(metric)
            #println("reset!(metric)",metric)
        catch
            AutomotiveDrivingModels.reset!(metric)
            #println("AutomotiveDrivingModels.reset!(metric)",metric)
        end
    end

    logl = 0.0
    n_traces = 0

    # Set correct action type
    if Symbol("a") in fieldnames(rand(model))
        simparams.ego_action_type = AccelTurnrate
        #println(" ego_action_type=AccelTurnrate")
    else
        simparams.ego_action_type = LatLonAccel
        println(" ego_action_type=LatLonAccel")
    end

    # simulate traces and perform online metric extraction
    scene = Scene()
    for seg_index in foldset_seg_test
        if (n_traces >= max_loop)
            break
        end
        println(now(), " ", seg_index, "/", max_loop)
        seg = simparams.segments[seg_index]
        trajdata = simparams.trajdatas[seg.trajdata_index]

        rec_orig = pull_record(seg, trajdata, prime_history) # TODO - make efficient
        rec_sim = deepcopy(rec_orig)

        time_start = get_time(trajdata, seg.frame_lo)
        time_end = get_time(trajdata, seg.frame_hi)
        n_traces += 1

        for sim_index in 1 : n_simulations_per_trace
            #debug = open("debug.log", "a"); print(debug, " ", sim_index, "/", n_simulations_per_trace);  close(debug)
            print(" ", sim_index, "/", n_simulations_per_trace)
            reset_simstate!(simparams.simstates[1], seg)

            #println("fieldnames=",fieldnames(model))
            if Symbol("net") in fieldnames(model)
                #debug = open("debug.log", "a"); println(debug, "rollout!()");  close(debug)
                rollout!(rec_sim, model, seg.egoid, trajdata,
                          time_start, time_end, simparams, prime_history=prime_history, gru_type=gru_type)
                #println("rollout!()")
            else
                #debug = open("debug.log", "a"); println(debug, "simulate!()");  close(debug)
                simulate!(rec_sim, model, seg.egoid, trajdata,
                          time_start, time_end, simparams, prime_history=prime_history)
                println("simulate!()")
            end

            for metric in metrics
                if !check_outliers(rec_orig, rec_sim, seg.egoid)
                    try
                        extract!(metric, rec_orig, rec_sim, trajdata.roadway, seg.egoid)
                        #debug = open("debug.log", "a"); println(debug, "extract!()");  close(debug)
                    catch
                        AutomotiveDrivingModels.extract!(metric, rec_orig, rec_sim, trajdata.roadway, seg.egoid)
                        #debug = open("debug.log", "a"); println(debug, "AutomotiveDrivingModels.extract!()");  close(debug)
                    end
                end
            end
        end
        println(" ")
    end

    # compute metric scores
    for metric in metrics
        try
            metrics_df[row, Symbol(metric)] = get_score(metric)
        catch
            metrics_df[row, Symbol(metric)] = AutomotiveDrivingModels.get_score(metric)
        end
    end
    metrics_df[row, :time] = string(now())
    metrics_df[row, :logl] = logl/n_traces

    metrics_df
end

max_stars_repo_name: yoshihikokuwahara/gail-driver

max_issues_repo_name: yoshihikokuwahara/gail-driver

max_forks_repo_name: yoshihikokuwahara/gail-driver


text:
#!/usr/bin/env julia

const detail_script_path = joinpath(dirname(@__FILE__()), "detail",  "parallel_worker.jl")

include(detail_script_path)

function addnodes()
    nodefile = get(ENV,"PBS_NODEFILE","")
    if !ispath(nodefile)
        error("nodefile ($nodefile) does not exist.")
    end
    nodes = Tuple{ASCIIString,Int}[(n,24) for n in split(readall(`sort -u $nodefile`))]
    if length(nodes) == 1
        ps = addprocs(topology=:master_slave)
    else
        ps = addprocs(nodes, topology=:master_slave)
    end
    @sync for p in ps
        @async remotecall_fetch(p, include, detail_script_path)
    end
    return ps
end

function makejobs(jobfile)
    jobs = Job[]
    pbs_id = parse(Int, split(ENV["PBS_JOBID"], ".")[1])
    io = open(jobfile)
    for (jobid, line) in enumerate(eachline(io))
        push!(jobs, Job(line, jobid, pbs_id))
    end
    close(io)
    return jobs
end

if length(ARGS) != 1
    info("usage: parallel.jl <jobfile>")
    exit()
elseif !ispath(ARGS[1])
    error("jobfile $(ARGS[1]) does not exist.")
end

jobs = makejobs(ARGS[1])
addnodes()
pmap(ex, jobs)

max_stars_repo_name: yomichi/SekireiTools

max_issues_repo_name: yomichi/SekireiTools

max_forks_repo_name: yomichi/SekireiTools


text:
"""
    AbstractFitnessStrategy

Base type for fitness strategies.

See [`fitnessfun`](@ref)
"""
abstract type AbstractFitnessStrategy end

"""
    fitnessfun(s::T, x, y) where T <: AbstractFitnessStrategy

Returns an `AbstractFitness` given the features `x` and labels `y`.

Allows for simple control over e.g. train vs validation split.
"""
fitnessfun(s::T, x, y) where T <: AbstractFitnessStrategy = error("Not implemented for $(T)!")

"""
    AbstractEvolutionStrategy

Base type for strategies on how to perform the evolution.

See [`evostrategy`](@ref)
"""
abstract type AbstractEvolutionStrategy end

"""
    evostrategy(s::AbstractEvolutionStrategy, inshape)

Returns an `AbstractEvolution`.

Argument `inshape` is the size of the input feature maps (i.e. how many pixels images are) and may be used to determine which mutation operations are allowed for example to avoid that feature maps accidentally become 0 sized.
"""
evostrategy(s::AbstractEvolutionStrategy, inshape) = evostrategy_internal(s, inshape)

"""
    struct TrainSplitAccuracy{T} <: AbstractFitnessStrategy
        TrainSplitAccuracy(;split, accuracyconfig, accuracyfitness, trainconfig, trainfitness)

Strategy to train model on a subset of the training data and measure fitness as accuracy on the rest.

Size of subset for accuracy fitness is `ceil(Int, split * nobs)` where `nobs` is the size of along the last dimension of the input.

# Arguments
- `accuracyconfig` (default [`BatchedIterConfig()`](@ref)) determine how to iterate over the accuracy subset.
- `accuracyfitness` (default [`AccuracyVsSize`](@ref)) determine how to measure fitness based on `accuracyconfig`.
- `trainconfig` (default [`TrainIterConfig()`](@ref)) determines how to iterate over the training subset.
- `trainfitness` is a function accepting the iterator produced by `trainconfig` and the fitness strategy produced
   by `accuracyfitness` and returns the `AbstractFitness` to be used (default [`TrainThenFitness`](@ref) wrapped in
   [`GpuFitness`](@ref)).
"""
struct TrainSplitAccuracy{S, VC, VF, TC, TF} <: AbstractFitnessStrategy
    split::S
    accconfig::VC
    accfitness::VF
    trainconfig::TC
    trainfitness::TF
end
function TrainSplitAccuracy(;split=0.1,
            accuracyconfig=BatchedIterConfig(),
            accuracyfitness=AccuracyVsSize,
            trainconfig=TrainIterConfig(),
            trainfitness=(iter, accf) -> GpuFitness(TrainThenFitness(StatefulGenerationIter(iter), Flux.Losses.logitcrossentropy, ADAM(), accf, 0.0)))

    return TrainSplitAccuracy(split, accuracyconfig, accuracyfitness, trainconfig, trainfitness)
end
function fitnessfun(s::TrainSplitAccuracy, x, y)
    train_x, acc_x = split_examples(x, ceil(Int, last(size(x)) * s.split))
    train_y, acc_y = split_examples(y, ceil(Int, last(size(y)) * s.split))

    aiter = dataiter(s.accconfig, acc_x, acc_y)
    afitness = s.accfitness(aiter)

    titer = dataiter(s.trainconfig, train_x, train_y)
    return s.trainfitness(titer, afitness)
end

split_examples(a::AbstractArray{T, 1}, splitpoint) where T = a[1:end-splitpoint], a[end-splitpoint+1:end]
split_examples(a::AbstractArray{T, 2}, splitpoint) where T = a[:,1:end-splitpoint], a[:,end-splitpoint+1:end]
split_examples(a::AbstractArray{T, 4}, splitpoint) where T = a[:,:,:,1:end-splitpoint], a[:,:,:,end-splitpoint+1:end]

"""
    AccuracyVsSize(data, accdigits=2, accwrap=identity)

Produces an `AbstractFitness` which measures fitness accuracy on `data` and based on number of parameters.

The two are combined so that a candidate `a` which achieves higher accuracy rounded to the first `accdigits` digits compared to a candidate `b` will always have a better fitness.

If the first `accdigits` of accuracy is the same the candidate with fewer parameters will get higher fitness.

Accuracy part of the fitness is calculated by `accwrap(`[`AccuracyFitness(data)`](@ref)`)`.
"""
AccuracyVsSize(data, accdigits=2; accwrap=identity) = sizevs(accwrap(AccuracyFitness(data)), accdigits)
function sizevs(f::AbstractFitness, accdigits=2)
    truncacc = MapFitness(x -> round(x, digits=accdigits), f)

    size = SizeFitness()
    sizefit = MapFitness(x -> min(10.0^-(accdigits+1), 1 / x), size)

    return AggFitness(+, truncacc, sizefit)
end

"""
    struct TrainAccuracyVsSize <: AbstractFitnessStrategy
    TrainAccuracyVsSize(;trainconfig, trainfitness)

Produces an `AbstractFitness` which measures fitness accuracy on training data and based on number of parameters combined in the same way as is done for [`AccuracyVsSize`](@ref).

# Arguments
- `trainconfig` (default [`TrainIterConfig()`](@ref)) determines how to iterate over the training subset.
- `trainfitness` is a function accepting the iterator produced by `trainconfig` and the fitness strategy produced
   by `accuracyfitness` and returns the `AbstractFitness` to be used (default [`TrainAccuracyFitness`](@ref) wrapped in
   [`GpuFitness`](@ref)).

Beware that fitness as accuracy on training data will make evolution favour overfitted candidates.
"""
struct TrainAccuracyVsSize{TC, TF} <: AbstractFitnessStrategy
    trainconfig::TC
    trainfitness::TF
end
function TrainAccuracyVsSize(;
                        trainconfig=TrainIterConfig(),
                        trainfitness = dataiter -> sizevs(GpuFitness(TrainAccuracyFitness(
                                                                            dataiter=StatefulGenerationIter(dataiter),
                                                                            defaultloss=Flux.Losses.logitcrossentropy, defaultopt = ADAM()))))
        return TrainAccuracyVsSize(trainconfig, trainfitness)
end
function fitnessfun(s::TrainAccuracyVsSize, x, y)
    iter = dataiter(s.trainconfig, x, y)
    return s.trainfitness(iter)
end

"""
    struct BatchedIterConfig{T, V}
    BatchedIterConfig(;batchsize=32, dataaug=identity, iterwrap=GpuIterator)

Configuration for creating batch iterators from array data.

The function `dataiter(s::BatchedIterConfig, x, y)` creates an iterator which returns a tuple of batches from `x` and `y` respectively.

More specifically, the result of `s.iterwrap(zip(s.dataaug(bx), by))` will be returned where `bx` and `by` are `BatchIterator`s.
"""
struct BatchedIterConfig{T, V}
    batchsize::Int
    dataaug::T
    iterwrap::V
end
BatchedIterConfig(;batchsize=32, dataaug=identity, iterwrap=GpuIterator) = BatchedIterConfig(batchsize, dataaug, iterwrap)
dataiter(s::BatchedIterConfig, x, y) = dataiter(x, y, s.batchsize, s.dataaug) |> s.iterwrap

"""
    struct ShuffleIterConfig{T, V}
    ShuffleIterConfig(;batchsize=32, seed=123, dataaug=identity, iterwrap=GpuIterator)

Configuration for creating shuffled batch iterators from array data. Data will be re-shuffled every time the iterator restarts.

The function `dataiter(s::ShuffleIterConfig, x, y)` creates an iterator which returns a tuple of batches from `x` and `y` respectively.

More specifically, the result of `s.iterwrap(Iterators.map(((x,y),) -> (s.dataaug(x), y), iter))` will be returned where `iter` is a `BatchIterator` over `x` and `y` with `shuffle=true`.

Note there there is no upper bound on how many generations are supported as the returned iterator cycles the data indefinitely. Use e.g. `Iterators.take(itr, cld(nepochs * nbatches_per_epoch, nbatches_per_gen))` to limit to `nepochs`
epochs.
"""
struct ShuffleIterConfig{T, V}
    batchsize::Int
    seed::Int
    dataaug::T
    iterwrap::V
end
ShuffleIterConfig(;batchsize=32, seed=123, dataaug=identity, iterwrap=GpuIterator) = ShuffleIterConfig(batchsize, seed, dataaug, iterwrap)
dataiter(s::ShuffleIterConfig, x, y) = dataiter(x, y, s.batchsize, s.seed, s.dataaug) |> s.iterwrap

"""
    struct TrainIterConfig{T}
    TrainIterConfig(nbatches_per_gen, baseconfig)
    TrainIterConfig(;nbatches_per_gen=400, baseconfig=ShuffleIterConfig())

Standard training strategy for creating a batched iterators from data. Data from `baseconfig` is cycled in partitions of `nbatches_per_gen` each generation using a [`RepeatPartitionIterator`](@ref).

"""
struct TrainIterConfig{T}
    nbatches_per_gen::Int
    baseconfig::T
end
TrainIterConfig(;nbatches_per_gen=400, baseconfig=ShuffleIterConfig()) = TrainIterConfig(nbatches_per_gen,baseconfig)
function dataiter(s::TrainIterConfig, x, y)
    baseiter = dataiter(s.baseconfig, x, y)
    return RepeatPartitionIterator(baseiter, s.nbatches_per_gen)
end

function dataiter(x::AbstractArray,y::AbstractArray{T, 1}, args...) where T
    _dataiter(x, Flux.onehotbatch(y, sort(unique(y))), args...)
end
dataiter(x::AbstractArray,y::AbstractArray{T, 2}, args...) where T = _dataiter(x,y, args...)

function _dataiter(x::AbstractArray,y::AbstractArray, bs::Integer, seed::Integer, xwrap, ywrap = identity)
    biter = shuffleiter((x,y), bs, seed)

    # iterates feature/label pairs
    # Cycle is pretty important here. It is what makes it so that each generation sees a new shuffle
    # This is also what makes SeedIterators useful with random data augmentation:
    #   - Even if SeedIterator makes the augmentation itself the same, it will be applied to different images
    # Something ought to be done about the inconsistency vs BatchedIterConfig through...
    citer = Iterators.map(biter) do (x, y)
        xwrap(x), ywrap(y)
    end |> Iterators.cycle
    # Ensures all models see the exact same examples, even when a RepeatStatefulIterator restarts iteration.
    return SeedIterator(citer;rng=biter.rng)
end
shuffleiter(data, batchsize, seed) = BatchIterator(data, batchsize; shuffle=MersenneTwister(seed))

function _dataiter(x::AbstractArray, y::AbstractArray, bs::Integer, xwrap, ywrap = identity)
    biter = BatchIterator((x,y), bs)
    # iterates feature/label pairs
    return Iterators.map(biter) do (x, y)
        xwrap(x), ywrap(y)
    end
end

"""
    struct GlobalOptimizerMutation{S<:AbstractEvolutionStrategy, F} <: AbstractEvolutionStrategy
    GlobalOptimizerMutation(base::AbstractEvolutionStrategy)
    GlobalOptimizerMutation(base::AbstractEvolutionStrategy, optfun)

Maps the optimizer of each candidate in a population through `optfun` (default `randomlrscale()`).

Basically a thin wrapper for `global_optimizer_mutation`.

Useful for applying the same mutation to every candidate, e.g. global learning rate schedules which all models follow.
"""
struct GlobalOptimizerMutation{S<:AbstractEvolutionStrategy, F} <: AbstractEvolutionStrategy
    base::S
    optfun::F
end
GlobalOptimizerMutation(base::AbstractEvolutionStrategy) = GlobalOptimizerMutation(base, NaiveGAflux.randomlrscale())

function evostrategy_internal(s::GlobalOptimizerMutation, inshape)
    base = evostrategy_internal(s.base, inshape)
    return AfterEvolution(base, pop -> NaiveGAflux.global_optimizer_mutation(pop, s.optfun))
end

"""
    struct EliteAndSusSelection <: AbstractEvolutionStrategy
    EliteAndSusSelection(popsize, nelites, evolve)
    EliteAndSusSelection(;popsize=50, nelites=2, evolve = crossovermutate())

Standard evolution strategy.

Selects `nelites` candidates to move on to the next generation without any mutation.

Also selects `popsize - nelites` candidates out of the whole population using [`SusSelection`](@ref) to evolve by applying random mutation.

Mutation operations are both applied to the model itself (change sizes, add/remove vertices/edges) as well as to the optimizer (change learning rate and optimizer algorithm).

Finally, models are renamed so that the name of each vertex of the model of candidate `i` is prefixed with "model`i`".
"""
struct EliteAndSusSelection{F} <: AbstractEvolutionStrategy
    popsize::Int
    nelites::Int
    evolve::F
end
EliteAndSusSelection(;popsize=50, nelites=2, evolve=crossovermutate()) = EliteAndSusSelection(popsize, nelites, evolve)

function evostrategy_internal(s::EliteAndSusSelection, inshape)
    elite = EliteSelection(s.nelites)

    evolve = SusSelection(s.popsize - s.nelites, EvolutionChain(ShuffleCandidates(), s.evolve(inshape)))

    combine = CombinedEvolution(elite, evolve)
    return AfterEvolution(combine, rename_models ∘ clear_redundant_vertices)
end

"""
    struct EliteAndTournamentSelection <: AbstractEvolutionStrategy
    EliteAndTournamentSelection(popsize, nelites, k, p, evolve)
    EliteAndTournamentSelection(;popsize=50, nelites=2; k=2, p=1.0, evolve = crossovermutate())

Standard evolution strategy.

Selects `nelites` candidates to move on to the next generation without any mutation.

Also selects `popsize - nelites` candidates out of the whole population using [`TournamentSelection`](@ref) to evolve by applying random mutation.

Mutation operations are determined by `evolve` both applied to the model itself (change sizes, add/remove vertices/edges) as well as to the optimizer (change learning rate and optimizer algorithm).

Finally, models are renamed so that the name of each vertex of the model of candidate `i` is prefixed with "model`i`".
"""
struct EliteAndTournamentSelection{R,F} <: AbstractEvolutionStrategy
    popsize::Int
    nelites::Int
    k::Int
    p::R
    evolve::F
end
EliteAndTournamentSelection(;popsize=50, nelites=2, k=2, p=1.0, evolve = crossovermutate()) = EliteAndTournamentSelection(popsize, nelites, k, p, evolve)

function evostrategy_internal(s::EliteAndTournamentSelection, inshape)
    elite = EliteSelection(s.nelites)

    evolve = TournamentSelection(s.popsize - s.nelites, s.k, s.p, s.evolve(inshape))

    combine = CombinedEvolution(elite, evolve)
    return AfterEvolution(combine, rename_models ∘ clear_redundant_vertices)
end

"""
    crossovermutate(;pcrossover=0.3, pmutate=0.9)

Return a function which creates an [`EvolutionChain`](@ref) when called with an inputshape.

Crossover will be applied with a probability of `pcrossover` while mutation will be applied with a probability of `pmutate`. Note that these probabilities apply to models only, not optimizers.

Crossover is done using [`CrossoverSwap`](@ref) for models and [`LearningRateCrossover`](@ref) and [`OptimizerCrossover`](@ref) for optimizers.

Mutation is applied both to the model itself (change sizes, add/remove vertices/edges) as well as to the optimizer (change learning rate and optimizer algorithm).
"""
crossovermutate(;pcrossover=0.3, pmutate=0.9) = function(inshape)
    cross = candidatecrossover(pcrossover)
    crossoverevo = AfterEvolution(PairCandidates(EvolveCandidates(cross)), align_vertex_names)

    mutate = candidatemutation(pmutate, inshape)
    mutationevo = EvolveCandidates(mutate)

    return EvolutionChain(crossoverevo, mutationevo)
end

candidatemutation(p, inshape) = evolvemodel(MutationProbability(graphmutation(inshape), p), optmutation())
candidatecrossover(p) = evolvemodel(MutationProbability(graphcrossover(), p), optcrossover())

function clear_redundant_vertices(pop)
    foreach(cand -> NaiveGAflux.model(check_apply, cand), pop)
    return pop
end

function align_vertex_names(pop)
    for i in eachindex(pop)
        m = match(r"model(\d+).\.*", modelname(pop[i]))
        if !isnothing(m)
            pop[i] = rename_model(m.captures[1], pop[i])
        end
    end
    return pop
end

function rename_models(pop)
    for i in eachindex(pop)
        pop[i] = rename_model(i, pop[i])
    end
    return pop
 end

function rename_model(i, cand)
    # No need to copy layers now
    return fmap(cand; walk=(f, x) -> x isa NaiveNASflux.AbstractMutableComp ? x : Functors._default_walk(f, x)) do x
        x isa String ? replace(x, r"^model\d+\.*" => "model$i.") : x
    end
end

function optcrossover(poptswap=0.3, plrswap=0.4)
    lrc = MutationProbability(LearningRateCrossover(), plrswap) |> OptimizerCrossover
    oc = MutationProbability(OptimizerCrossover(), poptswap) |> OptimizerCrossover
    return MutationChain(lrc, oc)
end

function optmutation(p=0.1)
    lrm = LearningRateMutation()
    om = MutationProbability(OptimizerMutation([Descent, Momentum, Nesterov, ADAM, NADAM, ADAGrad]), p)
    return MutationChain(lrm, om)
end

function graphcrossover()
    vertexswap = LogMutation(((v1,v2)::Tuple) -> "Crossover swap between $(name(v1)) and $(name(v2))", CrossoverSwap(0.05; mergefun=NaiveGAflux.default_mergefun(;layerfun=ActivationContributionLow)))
    crossoverop = VertexCrossover(MutationProbability(vertexswap, 0.05), 0.05)
    return LogMutation(((g1,g2)::Tuple) -> "Crossover between $(modelname(g1)) and $(modelname(g2))", crossoverop)
end

function graphmutation(inshape)
    acts = default_actfuns()

    change_nout = NoutMutation(-0.2, 0.2) # Max 20% change in output size
    decrease_nout = NoutMutation(-0.2, 0)
    add_vertex = add_vertex_mutation(acts)
    add_downsampling = AddVertexMutation(ConditionalArchSpace(candownsample(inshape), downsamplingspace(default_layerconf(); outsizes = 2 .^(4:9), activations=acts)))
    rem_vertex = RemoveVertexMutation()
    # [-2, 2] keeps kernel size odd due to CuArrays issue# 356 (odd kernel size => symmetric padding)
    change_kernel = KernelSizeMutation(ParSpace2D([-2, 2]), maxsize=maxkernelsize(inshape))
    decrease_kernel = KernelSizeMutation(ParSpace2D([-2]))
    mutate_act = ActivationFunctionMutation(acts)

    add_edge = AddEdgeMutation(0.1; mergefun=NaiveGAflux.default_mergefun(;layerfun=ActivationContributionLow))
    rem_edge = RemoveEdgeMutation()

    # Create a shorthand alias for MutationProbability
    mpn(m, p) = VertexMutation(MutationProbability(m, p))
    mph(m, p) = VertexMutation(HighUtilityMutationProbability(m, p))
    mpl(m, p) = VertexMutation(LowUtilityMutationProbability(m, p))

    cnout = mpn(LogMutation(v -> "Mutate output size of vertex $(name(v))", change_nout), 0.05)
    dnout = mpl(LogMutation(v -> "Reduce output size of vertex $(name(v))", decrease_nout), 0.05)
    maddv = mph(LogMutation(v -> "Add vertex after $(name(v))", add_vertex), 0.005)
    maddd = mpn(LogMutation(v -> "Add downsampling after $(name(v))", add_downsampling), 0.01)
    mremv = mpl(LogMutation(v -> "Remove vertex $(name(v))", rem_vertex), 0.01)
    ckern = mpn(MutationFilter(allowkernelmutation, LogMutation(v -> "Mutate kernel size of $(name(v))", change_kernel)), 0.02)
    dkern = mpn(MutationFilter(allowkernelmutation, LogMutation(v -> "Decrease kernel size of $(name(v))", decrease_kernel)), 0.02)
    mactf = mpl(LogMutation(v -> "Mutate activation function of $(name(v))", mutate_act), 0.005)
    madde = mph(LogMutation(v -> "Add edge from $(name(v))", add_edge), 0.02)
    mreme = mpn(MutationFilter(v -> length(outputs(v)) > 1, LogMutation(v -> "Remove edge from $(name(v))", rem_edge)), 0.02)

    mremv = MutationFilter(g -> nvertices(g) > 5, mremv)

    # Create two possible mutations: One which is guaranteed to not increase the size:
    dsize = MutationChain(mremv, dnout, dkern, mreme, maddd)
    # ...and another which is not guaranteed to decrease the size
    csize = MutationChain(mremv, cnout, ckern, mreme, madde, maddd, maddv)
    # Add mutation last as new vertices with neuronutility == 0 screws up outputs selection as per https://github.com/DrChainsaw/NaiveNASlib.jl/issues/39

    mgp = VertexMutation(MutationProbability(LogMutation(v -> "Mutate global pool type for $(name(v))", MutateGlobalPool()), 0.1), SelectGlobalPool())

    # If isbig then perform the mutation operation which is guaranteed to not increase the size
    # This is done mostly to avoid OOM and time outs.
    mall = MutationChain(MutationFilter(isbig, dsize), MutationFilter(!isbig, csize), mactf, mgp)

    return LogMutation(g -> "Mutate $(modelname(g))", mall)
end

isbig(g) = nparams(g) > 20e7

candownsample(inshape) = v -> candownsample(v, inshape)
candownsample(v::AbstractVertex, inshape) = is_convtype(v) && !infork(v) && canshrink(v, inshape)

function canshrink(v, inshape)
    # Note assumes stride = 2!
    # Also assumes single output after a global pool and flatten
    allvs = all_in_graph(v)
    gpv = allvs[findfirst(is_globpool, allvs)] |> inputs |> first
    return all(fshape(shapetrace(gpv) |> squashshapes, inshape) .> 1)
end

function infork(v, inputcnt = Dict{AbstractVertex, Int}(inputs(v) .=> 1), seen = Set())
    v in seen && return any(x -> x < 0, values(inputcnt))
    push!(seen, v)

    # How many times do we expect to find v as input to some subsequent vertex if we did not start inside a fork?
    inputcnt[v] = get(inputcnt, v, 0) + length(outputs(v))

    for vi in inputs(v)
        # And now subtract by 1 each time we find it
        inputcnt[vi] = get(inputcnt, vi, 0) - 1
    end

    foreach(vo -> infork(vo, inputcnt, seen), outputs(v))
    # If we have any "unaccounted" for inputs when we hit the output then we started inside a fork
    return any(x -> x < 0, values(inputcnt))
end

maxkernelsize(inshape) = v -> maxkernelsize(v, inshape)
function maxkernelsize(v::AbstractVertex, inshape)
    ks = fshape(shapetrace(v) |> squashshapes, inshape)
    # Kernel sizes must be odd due to CuArrays issue# 356 (odd kernel size => symmetric padding)
    return @. ks - !isodd(ks)
end

allowkernelmutation(v) = allowkernelmutation(NaiveNASflux.layertype(v), v)
allowkernelmutation(l, v) = false
allowkernelmutation(::NaiveNASflux.FluxConvolutional{N}, v) where N = all(isodd, size(NaiveNASflux.weights(layer(v)))[1:N])

function add_vertex_mutation(acts)

    function outselect(vs)
        rss = randsubseq(vs, 0.5)
        return isempty(rss) ? [rand(vs)] : rss
    end

    wrapitup(as) = AddVertexMutation(rep_fork_res(as, 1,loglevel=Logging.Info), outselect)

    add_conv = wrapitup(convspace(default_layerconf(),2 .^(4:9), 1:2:5, acts,loglevel=Logging.Info))
    add_dense = wrapitup(LoggingArchSpace(VertexSpace(default_layerconf(), NamedLayerSpace("dense", DenseSpace(2 .^(4:9), acts)));level = Logging.Info))

    return MutationChain(MutationFilter(is_convtype, add_conv), MutationFilter(!is_convtype, add_dense))
end

is_convtype(v::AbstractVertex) = any(is_globpool.(outputs(v))) || any(is_convtype.(outputs(v)))
is_globpool(v::AbstractVertex) = is_globpool(layer(v))
is_globpool(f) = f isa NaiveGAflux.GlobalPool

struct SelectGlobalPool <:AbstractVertexSelection
    s::AbstractVertexSelection
end
SelectGlobalPool() = SelectGlobalPool(AllVertices())
NaiveGAflux.select(s::SelectGlobalPool, g::CompGraph, ms...) = filter(is_globpool, NaiveGAflux.select(s.s, g, ms...))

struct MutateGlobalPool <: AbstractMutation{AbstractVertex} end
function (::MutateGlobalPool)(v::AbstractVertex)
    # Remove the old vertex and replace it with a new one with a different global pool type
    iv = inputs(v)[]
    RemoveVertexMutation()(v)
    v in outputs(iv) && return # failed to remove for some reason

    newname = first(splitext(name(v)))
    insert!(iv, vx -> GlobalPoolSpace(newpool(layer(v)))(newname, vx))
end

newpool(::GlobalPool{MaxPool}) = MeanPool
newpool(::GlobalPool{MeanPool}) = MaxPool


text:
using DependencyParser.DependencyParsing.ArcEager

# Configuration for "He wrote her a letter." after: SH, LA, RA, RA

function build_configuration()
  sentence = build_sentence()
  stack = Stack{Integer}()
  unshiftable = zeros(Bool, sentence.length)
  push!(stack, 0)
  push!(stack, 2)
  push!(stack, 3)

  buffer = Integer[4, 5, 6]
  tree = build_uncomplete_tree()

  Configuration(buffer, stack, sentence, unshiftable, tree)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
#############################
# Internal Structures
#############################

struct StaticPool <: AbstractThreadPool
    tids :: Vector{Int}
end

#############################
# Constructors
#############################

"""
    StaticPool(init_thrd=1, nthrds=Threads.nthreads())

The main StaticPool object.
"""
function StaticPool(init_thrd::Integer=1, nthrds::Integer=Threads.nthreads())
    thrd0 = min(init_thrd, Threads.nthreads())
    thrd1 = min(init_thrd+nthrds-1, Threads.nthreads())
    return StaticPool(thrd0:thrd1)
end

#############################
# ThreadPool Interface
#############################

function tmap(pool::StaticPool, fn::Function, itr)
    N = length(itr)
    result = Vector{_detect_type(fn, itr)}(undef, N)
    nts = length(pool.tids)
    n = div(N,nts)
    r = N % nts

    _fn = (tind) -> begin
        n0 = (tind-1)*n + 1 + (nts-tind+1 > r ? 0 : tind-nts+1)
        n1 = n0-1 + n + (nts-tind+1 <= r ? 1 : 0)
        for i in n0:n1
            @inbounds result[i] = fn(Base.unsafe_getindex(itr, i))
        end
    end

    Threads.@threads for tid in 1:Threads.nthreads()
        ind = findfirst(t->t==tid, pool.tids)
        isnothing(ind) || _fn(ind)
    end

    return result
end

#


text:
module PropertyUtils

export joinprops, fields, indexes, @with

#-----------------------------------------------------------------------------# joinprops
struct JoinProps{T}
    items::T
end

"""
    joinprops(items...)

Join items into a single struct with shared properties.  The first item that has the requested
property will be used.

# Example

    a = (x = 1, y = 2)
    b = (x = 3, z = 4)
    j = joinprops(a,b)
    j.x == 1
    j.z == 4
"""
joinprops(args...) = JoinProps(args)
Base.propertynames(j::JoinProps) = reduce(union, propertynames.(fields(j).items))
function Base.getproperty(j::JoinProps, x::Symbol)
    for item in getfield(j, :items)
        hasproperty(item, x) && return getproperty(item, x)
    end
    error("$j has no property $x")
end
function Base.setproperty!(j::JoinProps, name::Symbol, x)
    for item in getfield(j, :items)
        hasproperty(item, name) && return setproperty!(item, name, x)
    end
    error("`setproperty!(::JoinProps, args...) cannot create new properties.")
end

#-----------------------------------------------------------------------------# Indexes
struct Indexes{T}
    item::T
end

"""
    indexes(x)

Map `getproperty` to `getindex`.

# Example

    d = Dict(:x => 1, :y => 2)
    id = indexes(d)
    id.x
    id.z = 3
"""
indexes(x) = Indexes(x)
Base.propertynames(i::Indexes) = collect(keys(fields(i).item))
Base.getproperty(i::Indexes, x::Symbol) = getindex(getfield(i, :item), x)
Base.setproperty!(i::Indexes, name::Symbol, x) = setindex!(getfield(i, :item), x, name)

#-----------------------------------------------------------------------------# Fields
struct Fields{T}
    item::T
end

"""
    fields(x)

Map `getproperty` to `getfield`.

# Example
    struct A
        x::Int
    end
    Base.getproperty(a::A, x::Symbol) = "hello"

    a = A(1)
    a.x == "hello"

    fields(a).x == 1
"""
fields(x) = Fields(x)
Base.propertynames(f::Fields{T}) where {T} = fieldnames(T)
Base.getproperty(f::Fields, x::Symbol) = getfield(getfield(f, :item), x)
Base.setproperty!(f::Fields, name::Symbol, x) = setfield!(getfield(f, :item), name, x)

#-----------------------------------------------------------------------------# @with
"""
    @with src expr

For every symbol `x` in `expr`, replace it with `getproperty(src, x)` if `hasproperty(src, x)`.

# Example

    @with (x = 1, y = 2) x + y
"""
macro with(src, ex)
    esc(quote
        eval(PropertyUtils.replace_props!($src, $(Meta.quot(ex))))
    end)
end
replace_props!(df, x) = x
replace_props!(src, x::Symbol) = hasproperty(src, x) ? getproperty(src, x) : x
function replace_props!(src, ex::Expr)
    if ex.head === :call || ex.head === :.
        Expr(ex.head, vcat(ex.args[1], replace_props!.(Ref(src), ex.args[2:end]))...)
    else
        Expr(ex.head, replace_props!.(Ref(src), ex.args)...)
    end
end
end

macro with(src, interp, ex)
    @info interp
end


text:
mutable struct JacobiData{T <: Real, S, L <: LinearOperator} <: PreconditionerDataK2{T, S}
  P::L
  diagQ::S
  invDiagK::S
end

function Jacobi(
  id::QM_IntData,
  fd::QM_FloatData{T},
  regu::Regularization{T},
  D::AbstractVector{T},
  K::LinearOperator{T},
) where {T <: Real}
  invDiagK = (one(T) / regu.δ) .* fill!(similar(fd.c, id.nvar + id.ncon), one(T))
  diagQ = get_diag_Q_dense(fd.Q)
  invDiagK[1:(id.nvar)] .= .-one(T) ./ (D .- diagQ)
  P = opDiagonal(invDiagK)
  return JacobiData{eltype(diagQ), typeof(diagQ), typeof(P)}(P, diagQ, invDiagK)
end

function update_preconditioner!(
  pdat::JacobiData{T},
  pad::PreallocatedData{T},
  itd::IterData{T},
  pt::Point{T},
  id::QM_IntData,
  fd::QM_FloatData{T},
  cnts::Counters,
) where {T <: Real}
  if typeof(pad) <: PreallocatedDataK2_5Krylov
    pad.pdat.invDiagK[1:(id.nvar)] .=
      abs.(one(T) ./ (pad.D .- (pad.pdat.diagQ .* pad.sqrtX1X2 .^ 2)))
  else
    pad.pdat.invDiagK[1:(id.nvar)] .= abs.(one(T) ./ (pad.D .- pad.pdat.diagQ))
  end
  pad.pdat.invDiagK[(id.nvar + 1):end] .= one(T) / pad.regu.δ
end


text:
function varargout = regression(varargin)
#  REGRESSION MATLAB code for regression.fig
#       REGRESSION, by itself, creates a new REGRESSION or raises the existing
#       singleton*.
#
#       H = REGRESSION returns the handle to a new REGRESSION or the handle to
#       the existing singleton*.
#
#       REGRESSION('CALLBACK',hObject,eventData,handles,...) calls the local
#       function named CALLBACK in REGRESSION.M with the given input arguments.
#
#       REGRESSION('Property','Value',...) creates a new REGRESSION or raises the
#       existing singleton*.  Starting from the left, property value pairs are
#       applied to the GUI before regression_OpeningFcn gets called.  An
#       unrecognized property name or invalid value makes property application
#       stop.  All inputs are passed to regression_OpeningFcn via varargin.
#
#       *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
#       instance to run (singleton)".
#
#  See also: GUIDE, GUIDATA, GUIHANDLES

#  Edit the above text to modify the response to help regression

#  Last Modified by GUIDE v2.5 17-Apr-2013 12:46:04

#  Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @regression_OpeningFcn, ...
                   'gui_OutputFcn',  @regression_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
#  End initialization code - DO NOT EDIT

#  --- Executes just before regression is made visible.
function regression_OpeningFcn(hObject, ~, handles, varargin)
#  This function has no output args, see OutputFcn.
#  hObject    handle to figure
#  eventdata  reserved - to be defined in a future version of MATLAB
#  handles    structure with handles and user data (see GUIDATA)
#  varargin   command line arguments to regression (see VARARGIN)

#  Choose default command line output for regression
handles.output = hObject;

#  Update handles structure
guidata(hObject, handles);

#  UIWAIT makes regression wait for user response (see UIRESUME)
#  uiwait(handles.figure1);

# File: regression.jl
# Perform a linear regression
# Author:
# Date: 29-June-2020

# import XL data
[~,~,returns] = getXLData('hfma_matlab_data.xlsx','LS Index & S&P 500');
# S&P 500 (# )
x = returns(:,2);
# Long-Short Index (# )
y = returns(:,1);

# requested statistics
list = {'beta','adjrsquare','rsquare','tstat'};
# perform regression
stats = regstats(y,x,'linear',list);

# output statistics
cfs = stats.beta; # regression coefficients   β  \beta
arsq = stats.adjrsquare; # adjusted R square
rsq = stats.rsquare; # R square
t = stats.tstat.t; # t-statistic
pv = stats.tstat.pval; # p-value

# other statistics
se = cfs./t; # standard error
correl = corr(x,y); # correlation

# OUTPUT
# format data 1
data1 = num2cell([correl;arsq;rsq]);
    for i = 1:numel(data1)
            data1{i} = sprintf('# 12.4f',data1{i});
    end
# set column headings
headings1 = {'Correlation','adjusted R2','R2'};
# set table output
set(handles.table_regression_results_1,'ColumnName','','RowName',headings1,'Data',data1);
# format data 2
data2 = num2cell([cfs,se,t,pv]);
    for i = 1:numel(data2)
        data2{i} = sprintf('# 12.4f',data2{i});
    end
# set size column widths
cwidth = num2cell(repmat(70,1,10));
# set column headings
headings2 = {'Coefficients','SE','t-statistic','p-Value'};
# set table output
set(handles.table_regression_results_2,'ColumnName',headings2,'ColumnWidth',cwidth,'RowName',{'Intercept (a)','Slope (b) [S&P500]'},'Data',data2);

# plot linear regression
p = polyfit(x,y,1);
f = polyval(p,x);

figure()
plot(x,y,'.');
hold on
# linear fit
plot(x,f,'-r');
xlabel('S&P500 (# )');
ylabel('Long-Short Index (# )');
# add linear equation & R-squared
text(-1,-5,['y= ',num2str(cfs(2)),'x + ',num2str(cfs(1)),',  R^2= ',num2str(rsq)]);

#  --- Outputs from this function are returned to the command line.
function varargout = regression_OutputFcn(~, ~, handles)
#  varargout  cell array for returning output args (see VARARGOUT);
#  hObject    handle to figure
#  eventdata  reserved - to be defined in a future version of MATLAB
#  handles    structure with handles and user data (see GUIDATA)

#  Get default command line output from handles structure
varargout{1} = handles.output;


text:
# Silly example with slider
using Makie
sr, or = textslider(10:1000, "resolution", start = 20);

parent = Scene(resolution = (1000, 500))

t = lift(or) do resolution
  LinRange(0.0, 100.0, resolution)
end
y = lift(t) do t; sin.(t); end

pl = lines(t, y)
scene = hbox(sr, pl, parent=parent)


text:
mutable struct ROCK2ConstantCache{T,T2,zType} <: OrdinaryDiffEqConstantCache
  ms::SVector{46, Int}
  fp1::SVector{46, T}
  fp2::SVector{46, T}
  recf::Vector{T2}
  zprev::zType
  mdegprev::Int
  mdeg::Int
  recind::Int
end
@cache struct ROCK2Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  gprev::uType
  gprev2::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  k2::rateType
  constantcache::ROCK2ConstantCache
end

function alg_cache(alg::ROCK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = ROCK2ConstantCache(uEltypeNoUnits, uEltypeNoUnits, u) # WIP: not sure about what type to use in here
  gprev = similar(u)
  gprev2 = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  k2 = zero(rate_prototype)
  ROCK2Cache(u, uprev, gprev, gprev2, tmp, atmp, fsalfirst, k, k2, constantcache)
end

function alg_cache(alg::ROCK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  ROCK2ConstantCache(uEltypeNoUnits, uEltypeNoUnits, u) # WIP: not sure about what type to use in here
end

mutable struct ROCK4ConstantCache{T,T2,T3,T4,zType} <: OrdinaryDiffEqConstantCache
  ms::SVector{50, Int}
  fpa::Vector{T}
  fpb::Vector{T2}
  fpbe::Vector{T3}
  recf::Vector{T4}
  zprev::zType
  mdegprev::Int
  mdeg::Int
  recind::Int
end
@cache struct ROCK4Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  gprev::uType
  gprev2::uType
  gprev3::uType
  gprev4::uType
  gprev5::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  constantcache::ROCK4ConstantCache
end

function alg_cache(alg::ROCK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = ROCK4ConstantCache(real(uEltypeNoUnits), real(uEltypeNoUnits), u)
  gprev = similar(u)
  gprev2 = similar(u)
  gprev3 = similar(u)
  gprev4 = similar(u)
  gprev5 = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  ROCK4Cache(u, uprev, gprev, gprev2, gprev3, gprev4, gprev5, tmp, atmp, fsalfirst, k, k2, k3, k4, k5, constantcache)
end

function alg_cache(alg::ROCK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  ROCK4ConstantCache(real(uEltypeNoUnits), real(uEltypeNoUnits),u)
end


text:
using AbstractGUI
using Test
using Meshes
using WindowAbstractions

include("geometry.jl")
include("widget.jl")

function propagate_event(ws)
    function _propagate_event(ed::EventDetails)
        w = find_target(ws, ed)
        if !isnothing(w)
            execute_callback(w, ed)
        end
    end
end

cb_1 = WidgetCallbacks(on_key_pressed = x -> :rect1)
cb_2 = WidgetCallbacks(on_key_released = x -> :rect2)

rect1 = WRectangle(convert(Quadrangle, Box(P2(0,0),P2(256,256))),1,cb_1)
rect2 = WRectangle(convert(Quadrangle, Box(P2(50,50),P2(70,70))),2,cb_2)

ws = [rect1, rect2]

struct FakeWindowManager <: AbstractWindowManager end
struct FakeWindow <: AbstractWindow end

wm = FakeWindowManager()
win = FakeWindow()

cb = WindowCallbacks(on_key_pressed = propagate_event(ws), on_key_released = propagate_event(ws))

@testset "AbstractGUI.jl" begin
    @test captures_event(rect1, :on_key_pressed)
    @test captures_event(rect2, :on_key_released)

    ed = EventDetails(KeyEvent(:Z02, KeySymbol(:z), 'z', KeyModifierState(), KeyPressed()), (60,60), time(), win)
    @test find_target(ws, ed) == rect1
    @test execute_callback(cb, ed) == :rect1

    ed = EventDetails(KeyEvent(:Z02, KeySymbol(:z), 'z', KeyModifierState(), KeyReleased()), (60,60), time(), win)
    @test find_target(ws, ed) == rect2
    @test execute_callback(cb, ed) == :rect2

    ed = EventDetails(KeyEvent(:Z02, KeySymbol(:z), 'z', KeyModifierState(), KeyPressed()), (120,120), time(), win)
    @test find_target(ws, ed) == rect1
    @test execute_callback(cb, ed) == :rect1

    ed = EventDetails(KeyEvent(:Z02, KeySymbol(:z), 'z', KeyModifierState(), KeyReleased()), (120,120), time(), win)
    @test isnothing(find_target(ws, ed))
    @test isnothing(execute_callback(cb, ed))
end


text:
function j_piece(cur_container::AbstractContainer)
  cur_piece = Piece(
    cur_container,
    Block[],
    'j',
    "blue",
    3,
    0,
    -1,
    -1
  )

  cur_block_coords = [
    (-1, +1),
    (-1, 0),
    (0, 0),
    (+1, 0)
  ]

  make_blocks!(cur_piece, cur_block_coords)

  cur_piece
end


text:
const global options = ["Geology", "Geophysics", "Geobiology", "Geochemistry", "Planetary Science", "ESE"]

# find matching names, return matching score and index
function matchNames(student_name,fac_names)
    mc  = 0
    bestMatch = "" # best match!
    subname = split(strip(student_name)," ")
    for s in fac_names
        cc = 0
        fac_sub = split(strip(s)," ")
        for sub in lowercase.(subname), subf in lowercase.(fac_sub)
            # Could use regex or something, just being lazy...
            cc += occursin(sub, subf)
            cc += occursin(subf, sub)
        end
        if cc > mc
            mc = cc
            bestMatch = s
        end
    end
    return mc, bestMatch
end

"Load student table generated by Julie"
function loadStudents(xlsFile)
    students = Dict()
    for option in options
        students[option] = DataFrame(XLSX.readtable(xlsFile, option)...)
    end
    return students
end

"Check student submissions"
function checkStudents(students, submissions, thresh=3)
    for option in options
        for i = 1:size(students[option],1)
            foundEmail = students[option]."Email"[i] in submissions."Email Address"
            if !foundEmail
                # Try matching names:
                sName = students[option]."Student Name"[i]
                mc, bestMatch = matchNames(sName,submissions."Your Name")
                if mc<thresh
                    println(option, ": ", sName,  " ", students[option]."Year Of Study"[i], " , matchScore: ", mc, ", Best match: ", bestMatch )
                end
            end
        end
    end
end

"find matches for students (for a class and whether faculty and students matched in preference)"
function find_matched_students(df_classes, df_students, row, course)
    ind_matched = findall(df_classes."Class Name".==row[FacTA_class[course]])[1]
    #println(ind_matched)
    df_classes[ind_matched,3] = "X"
    #println("Selectable candidates (chosen by faculty): ", row[FacTA_student[course]])
    students = split(row[FacTA_student[course]],",")
    # @show students
    #println("Student Choices (Name, PriorityRank)")
    for a in keys(StuClass)
        matchedStudent = filter(df -> df[StuClass[a]]  == row[FacTA_class[course]], dropmissing(df_students, StuClass[a]))
        for i=1:size(matchedStudent,1)
            sName = matchedStudent[!,:"Your Name"][i]
            matches, _ = matchNames(sName,students)
            #@show sName

            if matches>1
                if df_classes[ind_matched,4] == ""
                    df_classes[ind_matched,4] = df_classes[ind_matched,4] *  sName
                else
                    df_classes[ind_matched,4] = df_classes[ind_matched,4] *  ", " * sName
                end

                #println(sName, " ", matchedStudent[!,StuInterest[a]][i], ", NameMatches: ", matches)
            else
                if df_classes[ind_matched,5] == ""
                    df_classes[ind_matched,5] = df_classes[ind_matched,5]  * sName
                else
                    df_classes[ind_matched,5] = df_classes[ind_matched,5]  * ", " * sName
                end

            end
        end
    end
end

"find matches for students (for a class, with faculty input missing)"
function find_matched_students(df_classes, df_students, class)
    ind_matched = findall(df_classes."Class Name".==class)[1]
    #println(ind_matched)
    df_classes[ind_matched,3] = "O"

    for a in keys(StuClass)
        matchedStudent = filter(df -> df[StuClass[a]]  == class, dropmissing(df_students, StuClass[a]))
        for i=1:size(matchedStudent,1)
            sName = matchedStudent[!,:"Your Name"][i]
            if df_classes[ind_matched,5] == ""
                df_classes[ind_matched,5] = df_classes[ind_matched,5]  * sName
            else
                df_classes[ind_matched,5] = df_classes[ind_matched,5]  * ", " * sName
            end
        end
    end
end

max_stars_repo_name: cfranken/TA-matching

max_issues_repo_name: cfranken/TA-matching

max_forks_repo_name: cfranken/TA-matching


text:
export activationForward
import Base: unsafe_convert
using Knet.KnetArrays: DevArray
using AutoGrad: AutoGrad, @primitive1

using CUDA.CUDNN: handle,
    cudnnActivationForward,
    cudnnActivationBackward,
    cudnnActivationDescriptor_t,
        cudnnCreateActivationDescriptor,
        cudnnSetActivationDescriptor,
        cudnnGetActivationDescriptor,
        cudnnDestroyActivationDescriptor,
    cudnnActivationMode_t,
        CUDNN_ACTIVATION_SIGMOID,
        CUDNN_ACTIVATION_RELU,
        CUDNN_ACTIVATION_TANH,
        CUDNN_ACTIVATION_CLIPPED_RELU,
        CUDNN_ACTIVATION_ELU,
        CUDNN_ACTIVATION_IDENTITY,
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN,
        CUDNN_PROPAGATE_NAN

@primitive1 activationForward(x; o...),dy,y  activationBackward(x,y,dy; o...)
@primitive1 activationBackward(x,y...;o...)  throw(MethodError(back,activationBackward))

function activationForward(x::R; alpha=1, o...) where {T,R<:DevArray{T}}
    y, td = similar(x), TD(T,(1,1,length(x),1))
    cudnnActivationForward(handle(), AD(; o...), Ref(T(alpha)), td, x, Ref(T(0)), td, y)
    return y
end

function activationBackward(x::R,y::R,dy::R; alpha=1, o...) where {T,R<:DevArray{T}}
    dx, td = similar(x), TD(T,(1,1,length(x),1))
    cudnnActivationBackward(handle(), AD(; o...), Ref(T(alpha)), td, y, td, dy, td, x, Ref(T(0)), td, dx)
    return dx
end

mutable struct ActivationDescriptor; ptr; end
unsafe_convert(::Type{<:Ptr}, ad::ActivationDescriptor)=ad.ptr
function AD(; # Defaults:
            mode::cudnnActivationMode_t=CUDNN_ACTIVATION_RELU,
            reluNanOpt::cudnnNanPropagation_t=CUDNN_NOT_PROPAGATE_NAN,
            coef::Float64=1.0
            )
    ptr = cudnnActivationDescriptor_t[C_NULL]
    cudnnCreateActivationDescriptor(ptr)
    cudnnSetActivationDescriptor(ptr[1], mode, reluNanOpt, coef)
    ad = ActivationDescriptor(ptr[1])
    finalizer(x->cudnnDestroyActivationDescriptor(x.ptr), ad)
    return ad
end


text:
import Base.getindex
import Base.length

abstract type AbstractCorrespondence end

struct Correspondence{N,T} <: AbstractCorrespondence
    x::NTuple{N,T}
end

function getindex(correspondence::Correspondence, index)
    getindex(correspondence.x, index)
end

function length(correspondence::Correspondence)
    length(correspondence.x)
end


text:
module MatrixImpute

using Distributed: pmap, @distributed, @everywhere
using LinearAlgebra: pinv, svd, norm, I, dot, diagm
using StatsBase: sample
using Arpack: svds
using DataFrames: allowmissing
using Statistics: mean

function Impute(A,k;method=:fastImpute,γ = 10e4, λ=0.1, B=nothing, lr= 1/64)

    dims = length(size(A))

    if dims == 2
        if method == :fastImpute
            Aopt = fastImpute(A,k,γ,lr,B, approx = true)
            return Aopt
        else
            error("Method not yet implemented!")
        end
    elseif dims ==3
        if method == :fastImputeT && B == nothing
            Aopt = fastImputeT(A,k,γ,lr)
            return Aopt
        elseif method == :fastImputeT
            error("Tensor with side information is not supported!")
        else
            error("Method not yet implemented!")
        end
    else
        error("Number of Dimensions Incorrect")
    end
end

function fastImpute(A,k,γ,lr,B; approx = true)
# Unknown values should be in missing for A
# A is the matrix that needs to be imputed, size n x m
# B is the feature selector, size m x p
# k is the maximum rank of the resultant matrix
# γ is the regularization parameter
# lr is the learning rate of the gradient descent function
    n = size(A)[1]
    m = size(A)[2]
    chunk = 1000
    if B == nothing
        p = m
        startchunk = collect(1:chunk:m)
        nchunks = length(startchunk)
        Uopt = zeros(Float64,n,k)
        Sopt = zeros(Float64,k,p)
        Aopt = zeros(Float64,n,p)
        if nchunks>1
            for sc in startchunk
                UoptT, SoptT = fastImputeInner(A[:,sc:min(sc+chunk-1,p)],k,γ,lr,B)
                Aopt[:,sc:min(sc+chunk-1,p)] = UoptT * SoptT
            end
            if approx
                return Aopt
            else
                output = svds(Aopt, nsv = k)
                Aopt = output[1].U*diagm(output[1].S)*output[1].Vt
                return Aopt
            end
        else
            Uopt, Sopt = fastImputeInner(A,k,γ,lr,B)
            return Uopt * Sopt
        end
    else
        if m != size(B)[1]
            error("Sizes of A and B must match")
        end
        Uopt, Sopt = fastImputeInner(A,k,γ,lr,B)
        return Uopt * Sopt
    end
end

function fastImputeInner(A,k,γ,lr,B)
    n = size(A)[1]
    m = size(A)[2]
    if B == nothing
        p = m
    elseif B != nothing
        p = size(B)[2]
    end
    M = sum(ismissing.(A))/(n*m)
    Uopt = zeros(n,k)
    W = zeros(n,m)
    for i = 1:n
        W[i,:] = ones(m) - convert(Vector{Float64},ismissing.(A[i,:]))
    end
    println("Preprocessing Complete")
    A[ismissing.(A)] .= 0
    if M<0.5 && B == nothing
        output = svds(A,nsv = k)
        S = diagm(output[1].S)*output[1].Vt
    else
        S = rand(k,p)
        S = S ./ norm(S)
    end
    objbest = Inf
    ∇obj = zeros(k,p)
    i = 0
    j = 0
    if B != nothing
        obj0, ∇obj0 = MatrixDeriv2(A,B,W,S,k,γ,M,j)
    else
        obj0, ∇obj0 = MatrixDeriv(A,W,S,k,γ,M,j)
    end
    noimprov = 0
    raised = false
    while i<50
        if (noimprov>=5) && (!raised)
            j = j + 1
            raised = true
            objbest = obj0
            noimprov = 0
        end
        ∇obj = ∇obj0 + i / (i+3) * ∇obj
        ∇tan = -∇obj + dot(∇obj,S) * S
        updategrad = ∇tan / norm(∇tan)
        S = S * cos(pi * lr) + updategrad*sin(pi * lr)
        S = S / norm(S)
        if B != nothing
            obj0, ∇obj0 = MatrixDeriv2(A,B,W,S,k,γ,M,j)
        else
            obj0, ∇obj0 = MatrixDeriv(A,W,S,k,γ,M,j)
        end
        if obj0 > objbest
            noimprov = noimprov + 1
            if (noimprov % 5 == 0) && (raised)
                raised = false
            end
        else
            noimprov = 0
            objbest = obj0
            raised = false
        end
        i = i + 1
    end
    Sopt = S
    if B == nothing
        X = Sopt'
    else
        X = B * Sopt'
    end
    println("Model Solved")
    function PopulateA(Xtemp,Atemp)
        return pinv(Xtemp'*Xtemp,1e-7)*(Xtemp'*Atemp)
    end
    Xtemp=Array{Float64}[X[W[i,:].==1,:] for i=1:n]
    Atemp=Array{Float64}[A[i,W[i,:].==1] for i=1:n]
    result=pmap(PopulateA,Xtemp,Atemp)
    for i=1:n
        Uopt[i,:]=result[i]
    end
    return Uopt, Sopt
end

function MatrixDeriv(A,W,S,k,γ,M,j)
    # Derivative Calculation for Matrix Completion without Side Information
    n = size(A)[1]
    m = convert(Int,size(A)[2])
    p = m
    nsquare = sqrt(n*m)
    mnew = m
    nnew = min(max(100,Int(round(nsquare * log(nsquare) * k * 2^j/(8 * mnew * (1-M))))),n)
    normfactor = 2 * mnew * nnew
    ∇obj = zeros(k,p)
    samplen=sample(1:n,nnew,replace=false)
    samplem=sample(1:m,mnew,replace=false)
    obj=0
    SmallInv = function (Xrow, Wrow, Arow)
        return Arow - Xrow * (inv(I / γ + Xrow' * Xrow) * (Xrow' * Arow))
    end
    Wpar = Array{Float64}[W[samplen[i],samplem] for i=1:nnew]
    Xpar = Array{Float64}[S[:,samplem[Wpar[i].==1]]' for i=1:nnew]
    Apar = Array{Float64}[A[samplen[i],samplem[Wpar[i].==1]] for i=1:nnew]
    objpar = map(SmallInv,Xpar,Wpar,Apar)
    obj = @distributed (+) for i = 1:nnew
        dot(Apar[i],objpar[i]) / normfactor
    end
    Btemp2 = zeros(k,nnew)
    for i = 1:nnew
        Btemp2[:,i] = Xpar[i]' * objpar[i]
        ∇obj[:,samplem[Wpar[i].==1]] += -2 * γ * Btemp2[:,i] * objpar[i]' / normfactor
    end
    return obj, ∇obj
end

function MatrixDeriv2(A,B,W,S,k,γ,M,j)
    # Derivative Calculation for Matrix Completion with Side Information
    p = size(B)[2]
    n = size(A)[1]
    m = convert(Int,size(A)[2])
    ∇obj = zeros(k,p)
    nsquare = sqrt(n * m)
    mnew = min(4*p,m)
    nnew = min(max(100,Int(round(nsquare * log(nsquare) * k * 2^j / (8 * mnew * (1-M))))),n)
    samplen = sample(1:n,nnew,replace=false)
    samplem = sample(1:m,mnew,replace=false)
    X = B[samplem,:]*S'
    obj = 0
    SmallInv = function (Xrow, Wrow, Arow)
        return Arow - Xrow * (inv(I / γ + Xrow' * Xrow) * (Xrow' * Arow))
    end
    Wpar = Array{Float64}[W[samplen[i],samplem] for i=1:nnew]
    Xpar = Array{Float64}[X[Wpar[i].==1,:] for i=1:nnew]
    Apar = Array{Float64}[A[samplen[i],samplem[Wpar[i].==1]] for i=1:nnew]
    objpar = map(SmallInv,Xpar,Wpar,Apar)
    obj = @distributed (+) for i = 1:nnew
        dot(Apar[i],objpar[i]) / (2*mnew*nnew)
    end
    Btemp1 = zeros(p,nnew)
    Btemp2 = zeros(k,nnew)
    for i = 1:nnew
        Btemp1[:,i] = B[samplem[Wpar[i].==1],:]' * objpar[i]
        Btemp2[:,i] = X[Wpar[i].==1,:]' * objpar[i]
    end
    ∇obj = @distributed (+) for i = 1:nnew
        -2 * γ * Btemp2[:,i] * Btemp1[:,i]' / (2*mnew*nnew)
    end
    return obj, ∇obj
end

function fastImputeT(A,k,γ,lr)
# Unknown values should be in missing for A
# A is the tensor that needs to be imputed, size n x m x p
# k is the maximum CP rank of the resultant matrix
# γ is the regularization parameter
# lr is the learning rate of the gradient descent function
    n,m,p = size(A)
    Amat = reshape(A,(n*m,p))
    W = ones(n,m,p)-ismissing.(A)
    Wmat = ones(n*m,p)-ismissing.(Amat)
    x = rand(n, k)
    y = rand(m, k)
    M = sum(ismissing.(A))/(n*m*p)

    z = [x; y]
    z = z ./ norm(z)
    x = z[1:n,:]
    y = z[n+1:end,:]
    i = 0
    j = 0
    println("Preprocessing Complete")
    obj0, ∇obj0x, ∇obj0y = TensorDeriv(Amat,Wmat,x,y,k,γ,M,j,size(A))

    ∇obj0 = [∇obj0x; ∇obj0y]
    ∇obj = zeros(m+n,k)

    objbest = Inf
    noimprov = 0
    raised = false
    while i<50
        if (noimprov>=5) && (!raised)
            j = j + 1
            raised = true
            objbest = obj0
            noimprov = 0
        end
        ∇obj = ∇obj0 + i / (i+3) * ∇obj
        ∇tan = - ∇obj + dot(∇obj,z) * z
        updategrad = ∇tan / norm(∇tan)
        z = z .* cos(pi * lr) + updategrad * sin(pi * lr)
        z = z ./ norm(z)
        x = z[1:n,:]
        y = z[n+1:end,:]
        obj0, ∇obj0x, ∇obj0y = TensorDeriv(Amat,Wmat,x,y,k,γ,M,j,size(A))
        ∇obj0 = [∇obj0x; ∇obj0y]
        if obj0 > objbest
            noimprov = noimprov + 1
            if (noimprov % 5 == 0) && (raised)
                raised = false
            end
        else
            noimprov = 0
            objbest = obj0
            raised = false
        end
        i = i + 1
    end
    X=zeros(m*n,k)
    function findbeta(Xrow, Arow)
        return inv(I/γ+Xrow'*Xrow)*Xrow'*Arow
    end
    for i=1:k
        X[:,i]=reshape(x[:,i]*y[:,i]',n*m)
    end
    function PopulateA(Xtemp,Atemp)
        return X*(pinv(Xtemp'*Xtemp,1e-7)*(Xtemp'*Atemp))
    end
    Aopt=zeros(n,m,p)
    Wpar=Array{Float64}[Wmat[:,i] for i=1:p]
    Xpar=Array{Float64}[X[Wpar[i].==1,:] for i=1:p]
    Apar=Array{Float64}[Amat[Wpar[i].==1,i] for i=1:p]
    result=map(PopulateA,Xpar,Apar)
    println("Model Solved")
    for i=1:p
        Aopt[:,:,i]=reshape(result[i],(n,m))
    end
    return Aopt
end

function TensorDeriv(Amat,Wmat,x,y,k,γ,M,j,dim)
    n,m,p = dim
    ncube = (n*m*p)^(1/3)
    nsamples = ncube * k * log(ncube)^2 * 4 * 2^j / (1-M)
    α = nsamples / (n*m*p)
    pnew = min(max(100, Int(round(α * p))), p)
    samplep = sample(1:p,pnew,replace=false)
    normfactor = n * m * pnew
    function SmallInv(Xrow, Wrow, Arow)
        return Arow- Xrow*(inv(I/γ+Xrow'*Xrow)*(Xrow'*Arow))
    end
    X = zeros(n*m,k)
    for i=1:k
        X[:,i] = reshape(x[:,i]*y[:,i]',n*m)
    end
    Wpar=Array{Float64}[Wmat[:,samplep[i]] for i=1:pnew]
    Xpar=Array{Float64}[X[Wpar[i].==1,:] for i=1:pnew]
    Apar=Array{Float64}[Amat[Wpar[i].==1,samplep[i]] for i=1:pnew]
    objpar=map(SmallInv,Xpar,Wpar,Apar)
    obj=@distributed (+) for i=1:pnew
        dot(Apar[i],objpar[i]) / normfactor
    end
    ∇objx=zeros(n,k)
    ∇objy=zeros(m,k)
    ∇obj=zeros(n*m,k)
    # Btemp1 = zeros(p, k)
    for i = 1:pnew
        # Btemp1[i, :] = objpar[i]' * Xpar[i]
        ∇obj[Wpar[i].==1,:] += -2 * γ * objpar[i] * (objpar[i]' * Xpar[i]) / normfactor
    end

    ∇obj = reshape(∇obj, (n,m,k))
    for i = 1:k
        ∇objx[:,i]=∇obj[:,:,i] * y[:,i]
    end
    for i = 1:k
        ∇objy[:,i]=∇obj[:,:,i]' * x[:,i]
    end

    return obj, ∇objx, ∇objy
end

export Impute

end # module