text:
module Junk
# TODO clean this up

using ..OpticSim
using BenchmarkTools
import Unitful

# space for junk code which should ultimately be deleted or cleaned up and moved elsewhere

function drawbezier()
    curve = TestData.onespanspline()
    drawcurve(curve, 100, 1000)
end

function drawbsplinecurve()
    curve = TestData.bsplinecurve()
    drawcurve(curve, 100, 1000)
end

function drawbeziersurface()
    vx = 0.0:0.1:1
    vy = 0.0:0.1:1
    surface = TestData.beziersurface()

    f(x, y) = point(surface, x, y)[3]
    scene = Makie.Scene(resolution = (1000, 1000))
    # One way to style the axis is to pass a nested dictionary / named tuple to it.

    # for i in 1:10000
    #     surface.controlpolygon.controlpoints[2,2].point[3] = i/100.0
    #     Makie.wireframe!(vx,vy,zvals)
    # end
    # record(scene, "surface.mp4", 1:255; framerate = 60) do i
    # surface.controlpolygon[2,2][3] = Float64(i)/40.0
    zvals = [f(x, y) for x in vx, y in vy]
    # Makie.wireframe!(vx,vy,zvals,axis = (frame = (linewidth = 2.0,),))
    # Makie.surface!(scene, vx, vy, f, axis = (frame = (linewidth = 2.0,),))
    Makie.wireframe!(vx, vy, zvals)
    # end
end

function plotcurveintersections()
    curve = TestData.curvecoefficients()
    maxtheta = 1.01

    let p
        curvepts = collect((evaluatecurve(curve, theta) for theta in 0.0:0.01:maxtheta))
        let x = Array{Real,1}(undef, 0), y = Array{Real,1}(undef, 0)
            for pt in curvepts
                push!(x, pt[1])
                push!(y, pt[2])
            end
            p = Plots.plot(x, y)
        end

        origin = [0, 0]

        for theta in 0.01:0.2:maxtheta
            endpoint = evaluatecurve(curve, theta)
            line = Ray([0.0, 0.0], endpoint)
            points = intersections(line, curve)
            maxalpha = 0

            for int in points
                pt = int[1]
                theta = int[3]

                if 0 <= theta <= maxtheta # only plot intersection points in the curve parameter range from 0..maxtheta
                    if maxalpha < int[2]
                        maxalpha = int[2]
                    end

                    ptx = [pt[1]]
                    pty = [pt[2]]
                    Plots.scatter!(p, ptx, pty, legend = false)
                end
            end

            line = hcat([0, 0], endpoint .* maxalpha)
            x = line'[:, 1]
            y = line'[:, 2]
            Plots.plot!(p, x, y)
        end
        Plots.display(p)
    end
end

function makiemeshexample5()
    points = [(0, 0, 0), (1, 0, 2), (2, 1, 0), (0, 2, 1)]
    color = [:red, :green, :blue, :yellow]

    # indices interpreted as triangles (every 3 sequential indices)
    indices = [1 2 3; 1 3 4; 1 4 2; 2 3 4]
    for i in 1:1
        Makie.mesh!(points, indices, color = color)
    end
end

function makiemeshexample4()
    points = [(0, 0, 0), (1, 0, 2), (2, 1, 0), (0, 2, 1)]
    color = [:red, :green, :blue, :yellow]

    # indices interpreted as triangles (every 3 sequential indices)
    indices = [1 2 3; 1 3 4; 1 4 2; 2 3 4]
    # replprint(indices)
    # for i in 1:1
    scene = Makie.Scene(resolution = (1000, 1000))
    Makie.mesh!(points, indices, color = color)
    scene
    # end
end

function makiemeshexample3()
    points = [(0, 0, 0), (1, 0, 2), (2, 1, 0), (0, 2, 1)]
    color = [:red, :green, :blue, :yellow]

    # indices interpreted as triangles (every 3 sequential indices)
    indices = [1, 2, 3, 1, 3, 4, 1, 4, 2, 2, 3, 4]
    Makie.mesh(points, indices, color = color)
end

function makiemeshexample2()
    x = [0, 1, 2, 0]
    y = [0, 0, 1, 2]
    z = [0, 2, 0, 1]
    color = [:red, :green, :blue, :yellow]
    i = [0, 0, 0, 1]
    j = [1, 2, 3, 2]
    k = [2, 3, 1, 3]
    # indices interpreted as triangles (every 3 sequential indices)
    indices = [1, 2, 3, 1, 3, 4, 1, 4, 2, 2, 3, 4]
    Makie.mesh(x, y, z, indices, color = color)
end

function makiemeshexample()
    coordinates = [
        0.0 0.0
        0.5 0.0
        1.0 0.0
        0.0 0.5
        0.5 0.5
        1.0 0.5
        0.0 1.0
        0.5 1.0
        1.0 1.0
    ]
    connectivity = [
        1 2 5
        1 4 5
        2 3 6
        2 5 6
        4 5 8
        4 7 8
        5 6 9
        5 8 9
    ]
    color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
    scene = Makie.mesh(coordinates, connectivity, color = color, shading = false)

    Makie.wireframe!(scene[end][1], color = (:black, 0.6), linewidth = 3)
end

# function simplerender(samples::Int = 100)
#     cyl = Cylinder(0.3, 2.0)
#     cyl2 = Cylinder(0.2, 2.0)
#     rbt = translation(0.5, 0.5, 0.0)
#     # generator = csgintersection(leaf(cyl),leaf(cyl2,rbt))
#     generator = csgunion(leaf(cyl), leaf(cyl2, rbt))
#     # generator = leaf(cyl2,rbt)

#     # surf = rbt*AcceleratedParametricSurface(TestData.beziersurface(),10)
#     # generator = csgintersection(leaf(cyl),leaf(surf))

#     csg = generator(identitytransform())
#     view = Vis.CSGViewer([3.0, 3.0, 20.0], 2.0, samples, (csg,))
#     image = Vis.render(view)
#     Vis.show(image)
# end

# function csgrender(samples::Int = 100)
#     # generator = TestData.planecylindercsg()
#     # csg = generator(identitytransform())

#     # generator = TestData.trianglemeshcsg()
#     # csg = generator(RigidBodyTransform(-.5, -.5, 0.0))
#     # generator = TestData.csgobject()
#     # csg = generator(RigidBodyTransform(-.5, -.5, 0.0))

#     generator = TestData.lensarray(2)
#     csg = generator(identitytransform())
#     # view = CSGViewer(SVector{3}(0.0, 60.0, 80.0), .5, samples, (csg,))
#     # generator = leaf(AcceleratedParametricSurface(beziersurface(),1

#     # generator = leaf(Cylinder(.3,10.0))
#     # csg = generator(RigidBodyTransform(0.0,0.0,0.0))

#     # csg = AcceleratedParametricSurface(beziersurface(), 10)
#     view = Vis.CSGViewer([0.0, 4.0, 8.0], 1.4, samples, (csg,))
#     image = Vis.render(view)

#     # Vis.show(image)
# end

function timebezierpoint()
    surf = TestData.beziersurface()
    res = @benchmark point($surf, 0.5, 0.5)
    show(res)
end

function timereflectedray()
    res = @benchmark TestData.reflectedray(rand(3), rand(3))
    println(" Array{Float64,1} time")
    show(res)
    res = @benchmark TestData.reflectedray(SVector{3,Float64}(rand(3)), SVector{3,Float64}(rand(3)))
    println("SVector time")
    show(res)
end

function timebezierrayintersection()
    u = 0.5
    v = 0.5
    surf, r = TestData.beziersurfaceandray(u, v)
    surf = AcceleratedParametricSurface(surf, 5)
    bnch = @benchmark surfaceintersection($surf, $r)
    println(surfaceintersection(surf, r), " =? ", point(surf, u, v))
    println()
    show(bnch)
end

function timebboxrayintersection()
    bbox, r = TestData.bboxray()
    println(doesintersect(bbox, r))
    res = @benchmark doesintersect($bbox, $r)
    show(res)
end

function timecylinderrayintersection()
    cyl = Cylinder(0.5)
    r = Ray([1.0, 1.0, 1.0], [-1.0, -1.0, -1.0])
    bnch = @benchmark surfaceintersection($cyl, $r)
    println(surfaceintersection(cyl, r))
    println()
    show(bnch)
end

function timesphererayintersection()
    sph = Sphere(0.5)
    r = Ray([1.0, 1.0, 1.0], [-1.0, -1.0, -1.0])
    bnch = @benchmark surfaceintersection($sph, $r)
    println(surfaceintersection(sph, r))
    println()
    show(bnch)
end

function timetrianglerayintersection()
    tri, r = TestData.triangleray()
    res = @benchmark surfaceintersection($tri, $r)
    println(surfaceintersection(tri, r))
    println()
    show(res)
end

function bboxrayintersection(numits)
    bbox, r = TestData.bbox()

    for i in 1:numits
        res = surfaceintersection(bbox, r)
    end
end

function timeplanerayintersection()
    pln, r = TestData.planeray()
    res = @benchmark surfaceintersection($pln, $r)
    println(surfaceintersection(pln, r))
    println()
    show(res)
end

function onaxissphereintersection()
    println(surfaceintersection(Sphere(1.0), Ray([0.0, 0.0, -4.0], [0.0, 0.0, 1.0])))
    println(surfaceintersection(Sphere(1.0), Ray([0.0, 0.0, 4.0], [0.0, 0.0, -1.0])))
    println(surfaceintersection(Sphere(1.0), Ray([0.0, 4.0, 0.0], [0.0, -1.0, 0.0])))
    println(surfaceintersection(Sphere(1.0), Ray([0.0, -4.0, 0.0], [0.0, 1.0, 0.0])))
    println(surfaceintersection(Sphere(1.0), Ray([4.0, 0.0, 0.0], [-1.0, 0.0, 0.0])))
    println(surfaceintersection(Sphere(1.0), Ray([-4.0, 0.0, 0.0], [1.0, 0.0, 0.0])))
end

function makeimagedata(samples)
    zeropix = (0.0f0, 0.0f0, 0.0f0)
    return (fill(zeropix, samples, samples), fill(zeropix, samples, samples))
end

function threadingtest(samples)
    image, image2 = makeimagedata(samples)

    function innerloop(x, image, image2, samples, convolutionsize)
        for y in (convolutionsize + 1):samples
            for l in 1:convolutionsize
                for m in 1:convolutionsize
                    pix = image[y - l, x - m]
                    oldpix = image2[y, x]
                    image2[y, x] = (oldpix[1] + pix[1] * 0.8, oldpix[2] + pix[2] * 0.8, oldpix[3] + pix[3] * 0.8)
                end
            end
        end
    end

    convolutionsize = 30
    # println("multithreaded")
    # @sync
    for x in (convolutionsize + 1):samples
        #  for x in convolutionsize + 1:samples
        #    Threads.@spawn
        innerloop(x, image, image2, samples, convolutionsize)
    end
    return image
end

function timeconvolution()
    result = @benchmark threadingtest(1000)
    show(result)
end

function testdrawmesh()
    Vis.draw(TestData.beziersurface())
end

function testintervalsforcsg()
    surf1 = AcceleratedParametricSurface(TestData.beziersurface(), 5)
    surf2 = AcceleratedParametricSurface(TestData.upsidedownbeziersurface(), 5)
    r = Ray([0.5, 0.5, 2.0], [0.0, 0.0, -1.0])
    int1 = DisjointUnion(surfaceintersection(surf1, r)...)
    int2 = DisjointUnion(surfaceintersection(surf2, r)...)
    intsct = intervalintersection(int1, int2)
    if !(intsct isa EmptyInterval)
        s = Vis.scene()
        Vis.draw!(s, surf1)
        Vis.draw!(s, surf2)
        Vis.draw!(s, intsct)
        Vis.display(s)
    else
        throw(ErrorException("This should intersect (testintervalsforcsg)"))
    end
end

function testclippedbeziersurface()
    r = Ray([0.5, 0.5, 3.0], [0.0, 0.0, -1.0])
    gen = TestData.clippedbeziersurface()
    csg = gen(identitytransform())
    intsct = OpticSim.evalcsg(csg, r)
    # println(intsct)
    view = Vis.CSGViewer(SVector{3}(2.0, 2.0, 5.0), 2.0, 1000, (csg,))
    Vis.imshow(Vis.render(view))
end

function testdrawcylinder()
    surf = Cylinder(1.0)
    s = Vis.scene()
    Vis.draw!(s, surf, numdivisions = 100)
    r = Ray([0.0, 2.0, 0.0], [0.0, -1.0, 0.5])
    int = surfaceintersection(surf, r)
    Vis.draw!(s, int)
    Vis.display(s)
end

function testdrawsphere()
    surf = Sphere(1.0)
    s = Vis.scene()
    Vis.draw!(s, surf, numdivisions = 100)
    r = Ray([0.0, 2.0, 0.0], [0.0, -1.0, 0.5])
    int = surfaceintersection(surf, r)
    Vis.draw!(s, int)
    Vis.display(s)
end

function testcsgeval()
    surf1 = AcceleratedParametricSurface(TestData.beziersurface(), 5)
    surf2 = AcceleratedParametricSurface(TestData.upsidedownbeziersurface(), 5)
    r = Ray([0.5, 0.5, 2.0], [0.0, 0.0, -1.0])
    gen = csgintersection(leaf(surf1), leaf(surf2))
    csgtree = gen(identitytransform())
    intsct = surfaceintersection(csgtree, r)
    if !(intsct isa EmptyInterval)
        s = Vis.scene()
        Vis.draw!(s, surf1)
        Vis.draw!(s, surf2)
        Vis.draw!(s, intsct)
        Vis.display(s)
    else
        throw(ErrorException("This should intersect (testcsgeval)"))
    end
end

function testcsgwithmultipleobjects()
    generator, _, surf1, surf2 = TestData.csgobject()
    csg = generator(identitytransform())
    r = Ray{Float64,3}(SVector(0.5, 0.5, 4.0), SVector(0.0, 0.0, -1.0))

    intscts = OpticSim.evalcsg(csg, r)[1]
    # println(intscts)
    s = Vis.scene()
    Vis.draw!(s, surf1)
    Vis.draw!(s, surf2)
    Vis.draw!(s, intscts)
    Vis.display(s)
end

function testmeshcsgsurf1(subdiv = 20)
    _, _, surf, _ = TestData.csgobject()
    m = csgintersection(leaf(surf, translation(-0.5, -0.5, 0.0)), csgintersection(leaf(Cylinder(0.3, 5.0)), leaf(surf, RigidBodyTransform(rotmatd(0, 180, 0), SVector(0.5, -0.5, 0.0)))))()
    b = @benchmark makemesh($m, $subdiv)
    Vis.draw(m)
    return b
end

function testmeshcsgsurf2(subdiv = 20)
    _, _, surf1, surf2 = TestData.csgobject()
    m = csgintersection(leaf(surf1, translation(-0.5, -0.5, 0.0)), csgintersection(leaf(Cylinder(0.3, 5.0)), leaf(surf2, translation(-0.5, -0.5, 0.0))))()
    b = @benchmark makemesh($m, $subdiv)
    Vis.draw(m)
    return b
end

function testmeshcsguni(subdiv = 20)
    m = csgunion(leaf(Cylinder(0.5, 3.0)), (leaf(Sphere(1.0))))()
    b = @benchmark makemesh($m, $subdiv)
    Vis.draw(m)
    return b
end

function testmeshcsgint(subdiv = 20)
    m = csgintersection(leaf(Cylinder(0.5, 3.0)), (leaf(Sphere(1.0))))()
    b = @benchmark makemesh($m, $subdiv)
    Vis.draw(m)
    return b
end

function testmeshcsgsub(subdiv = 20)
    m = csgdifference(leaf(Cylinder(0.5, 3.0)), leaf(Sphere(1.0)))()
    b = @benchmark makemesh($m, $subdiv)
    Vis.draw(m)
    return b
end

function testmeshcsgmla(subdiv = 20)
    topsphere = leaf(Sphere(2.0), translation(0.05, 0.0, -1.5))
    botsphere = leaf(Sphere(3.0), translation(0.0, 0.1, 1.5))
    lens = csgintersection(leaf(Cylinder(0.7, 5.0)), csgintersection(botsphere, topsphere))
    d = 0.8
    MLA = lens
    for x in 0:2
        for y in 0:2
            if x == 0 && y == 0
                continue
            end
            MLA = csgunion(MLA, leaf(lens, RigidBodyTransform(rotmatd(0, 0, rand(Int) % 180), SVector{3}(x * d, y * d, 0.0))))
        end
    end
    m = MLA(identitytransform())
    b = @benchmark makemesh($m, $subdiv)
    Vis.draw(m)
    return b
end

function testorthogonalitymatrix()
    x = [0 4 -4.5 4.5; 0 11.25 -31.5 20.25]
    correctCmatrix = [
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
        4.0 0.0 0.0 0.0 11.25 0.0 0.0 0.0 0.0 1.0 0.0 0.0
        -4.5 4.0 0.0 0.0 -31.5 11.25 0.0 0.0 0.0 0.0 1.0 0.0
        4.5 -4.5 4.0 0.0 20.25 -31.5 11.25 0.0 0.0 0.0 0.0 1.0
        0.0 4.5 -4.5 4.0 0.0 20.25 -31.5 11.25 0.0 0.0 0.0 0.0
        0.0 0.0 4.5 -4.5 0.0 0.0 20.25 -31.5 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 4.5 0.0 0.0 0.0 20.25 0.0 0.0 0.0 0.0
    ]
    cmat = OpticSim.orthogonalitymatrix(x, 3)
    show(IOContext(stdout), "text/plain", cmat)
    for i in CartesianIndices(cmat)
        if cmat[i] != correctCmatrix[i]
            println("did not pass test")
            break
        end
    end
end

function testtriangulate()
    surf = TestData.beziersurface()
    numtris = 20 # number of triangles per row or column
    tris = triangulate(surf, numtris)
    tmesh = TriangleMesh(tris)
    Vis.draw(tmesh)
end

function testcurvetobeziersegments()
    orig = TestData.homogeneousbsplinecurve()
    segments = tobeziersegments(orig)
    allsegments = [BezierCurve{OpticSim.Rational,Float64,3,3}(segment) for segment in segments]
    Vis.draw(orig, allsegments...)
end

function testsurfacetobeziersegments()
    orig = TestData.bsplinesurface()
    segments = tobeziersegments(orig)
    Vis.draw(segments...)
end

function testinsertknot()
    orig = TestData.bsplinecurve()
    curve1 = TestData.bsplinecurve()
    let curve2 = curve1
        for i in 1:2
            curve2 = insertknot(curve1, 5 + i - 1)
            curve1 = curve2
        end
        Vis.draw(orig, curve2)
    end
end

# function timelinalg(vec::Bool = true, mat::Bool = true)
#     # pf(q) = print(BenchmarkTools.prettytime(time(median(q))), "\n")

#     if vec
#         print("VECTOR OPS\n")
#         vec3a = @SVector rand(3)
#         vec3b = @SVector rand(3)
#         vec1000a = @SVector rand(1000)
#         vec1000b = @SVector rand(1000)

#         for (vec1, vec2) in [(vec3a, vec3b), (vec1000a, vec1000b)]
#             print("SIZE: $(length(vec1))\n")
#             # NORMS
#             print("NORMALIZE\n")
#             @btime (OpticSim.mnormalize($vec1))
#             @btime (LinearAlgebra.normalize($vec1))

#             # DOT
#             print("DOT\n")
#             @btime (OpticSim.dot($vec1, $vec2))
#             @btime (LinearAlgebra.dot($vec1, $vec2))

#             # CROSS
#             if length(vec1) == 3
#                 print("CROSS\n")
#                 @btime (OpticSim.cross($vec1, $vec2))
#                 @btime (LinearAlgebra.cross($vec1, $vec2))
#             end
#         end
#     end

#     if mat
#         print("MATRIX OPS\n")
#         mat3a = @SMatrix rand(3, 3)
#         mat3b = @SMatrix rand(3, 3)
#         mat100a = @SMatrix rand(12, 12)
#         mat100b = @SMatrix rand(12, 12)
#         vec100 = @SVector rand(12)

#         for (mat1, mat2) in [(mat3a, mat3b), (mat100a, mat100b)]
#             print("SIZE: $(size(mat1))\n")

#             # DETERMINANT
#             print("DET\n")
#             @btime (OpticSim.det($mat1))
#             @btime (LinearAlgebra.det($mat1))

#             # INVERSE
#             print("INV\n")
#             @btime (OpticSim.minv($mat1))
#             @btime (LinearAlgebra.inv($mat1))

#             # DIVISON
#             print("DIV\n")
#             @btime (OpticSim.mdiv($mat1, $mat2))
#             @btime ($mat1 \ $mat2)
#         end
#     end
# end

function testniandnt()
    r = Ray([0.0, 0.0, 7.0], [0.0, 0.0, -1.0])
    normal = SVector{3,Float64}(0.0, 0.0, 1.0)

end

function testrefractedray()
    elt = TestData.planoconvexelement(OpticSim.GlassCat.SCHOTT.N_BK7, 0.0, 5.0, 60.0)
    lens = OpticSim.LensAssembly(elt.lens)
    r = Ray([0.0, 0.0, 7.0], [0.0, 0.0, -1.0])

    green = 500 * Unitful.u"nm"
    glass = OpticSim.GlassCat.SCHOTT.BAK50

    intsct = OpticSim.closestintersection(lens, r)

    pt = point(intsct[1])
    nml = normal(intsct[1])

    rdir = direction(r)
    (nᵢ, nₜ) = OpticSim.nᵢandnₜ(index(OpticSim.GlassCat.Air, green), index(glass, green), nml, r)

    refracted = OpticSim.refractedray(nᵢ, nₜ, nml, rdir)
end

badray() = Ray([5.000000042187876, 5.000000042187876, 0.6143579891640621], [0.042187876699796796, 0.042187876699796796, -0.9982185963600987])

function printtrace(raytrace)
    for trace in raytrace
        println("RAY***************")
        println(trace)
        println("END RAY***********")
    end
end
function testdoubleconcave()
    system = TestData.doubleconcave()
    green = 500 * Unitful.u"nm"
    r1 = Ray([5.0, 5.0, 2.0], [0.0, 0.0, -1.0])
    r2 = badray()
    elt = system.system.assembly.elements[1]
    # intsct = OpticSim.surfaceintersection(elt,r)
    # println(intsct)
    for r in (r1,)
        allrays = Array{OpticSim.LensTrace{Float64,3},1}(undef, 0)

        res = OpticSim.trace(system, r, green, trackrays = allrays)
        printtrace(allrays)
    end
end

function testplane()
    rparallel = Ray([0.0, 0.0, 2.0], [1.0, 1.0, 0.0])
    rinplane = Ray([0.0, 0.0, 1.0], [1.0, 1.0, 0.0])
    routside = Ray([0.0, 0.0, 2.0], [1.0, 1.0, 1.0])
    rinside = Ray([0.0, 0.0, -1.0], [1.0, 1.0, -1.0])
    routintersects = Ray([0.0, 0.0, 2.0], [0.0, 0.0, -1.0])
    rinintersects = Ray([0.0, 0.0, -1.0], [0.0, 0.0, 1.0])
    pln = Plane(0.0, 0.0, 1.0, 0.0, 0.0, 1.0)

    surfaceintersection(pln, routside)
    if !(nothing === surfaceintersection(pln, rparallel))
        @error "should have returned nothing for ray parallel to plane and outside"
    else
        if !(nothing === surfaceintersection(pln, routside))
            @error "shouldhave returned nothing for entirely outside plane but not parallel"
        else
            # NOTE for coplanar faces to work with visualization, rays in the plane must count as being 'inside'
            res = surfaceintersection(pln, rinplane)
            if !(lower(res) isa RayOrigin && upper(res) isa Infinity)
                @error "should have returned (rayorigin,∞)"
                return false
            else
                res = surfaceintersection(pln, rinside)
                if !(lower(res) isa RayOrigin && upper(res) isa Infinity)
                    @error "should have returned (rayorigin,∞)"
                    return false
                else
                    res = surfaceintersection(pln, routintersects)
                    if !(samepoint(point(lower(res)), [0.0, 0.0, 1.0]) && upper(res) isa Infinity)
                        @error "should have returned ([0.0,0.0,0.0],∞)"
                        return false
                    else
                        res = surfaceintersection(pln, rinintersects)
                        if !(lower(res) isa RayOrigin && samepoint(point(upper(res)), [0.0, 0.0, 1.0]))
                            @error "should have returned(rayorigin,[0.0,0.0,0.0])"
                            return false
                        else
                            return true
                        end
                    end
                end
            end
        end
    end
end

function testintervalcomplement()
    intvl = surfaceintersection((leaf(Sphere(1.0)))(identitytransform()), Ray([2.0, 0.0, 0.0], [-1.0, 0.0, 0.0]))
    comp = OpticSim.intervalcomplement(intvl)
    println(comp)
end

function testcsgoperations()
    pln1 = Plane(SVector{3}(0.0, 0.0, -1.0), SVector{3}(0.0, 0.0, -1.0))
    pln2 = Plane(SVector{3}(0.0, 0.0, 1.0), SVector{3}(0.0, 0.0, 1.0))
    r = Ray{Float64,3}(SVector{3}(0.0, 0, 2.0), SVector{3}(0.0, 0.0, -1.0))
    gen1 = csgintersection(pln1, pln2)
    gen2 = csgunion(pln1, pln2)

    csg = gen1(identitytransform())
    intsct = surfaceintersection(csg, r)

    intvlpt1 = lower(intsct)
    intvlpt2 = upper(intsct)
    println(intvlpt1, "  ", intvlpt2)
    if !(isapprox(point(intvlpt1), SVector{3}(0.0, 0.0, 1.0)) && isapprox(point(intvlpt2), SVector{3}(0.0, 0.0, -1.0)))
        @error "Should have two intersections with the two planes"
        return false
    end
    csg = gen2(identitytransform())
    intsct = surfaceintersection(csg, r)
    if !((lower(intsct) isa RayOrigin) && (upper(intsct) isa Infinity))
        @error "Union space should be half infinite"
        return false
    end
    return true
end

function testinsideplane()
    r = Ray([0.0, 0.0, 1.0], [0.0, 0.0, -1.0])
    p = Plane(0.0, 0.0, 1.0, 0.0, 0.0, 1.5)
    res = surfaceintersection(p, r)
    println(res)
end

function testsinglenegsurface()
    frontradius = -60.0
    frontvertex = 0.0
    semidiameter = 9.0
    r = Ray([5.0, 5.0, 1.0], [0.0, 0.0, -1.0])

    temp = leaf(Sphere(abs(frontradius)), OpticSim.translation(0.0, 0.0, frontvertex - frontradius))
    d = -frontradius - sqrt(frontradius^2 - semidiameter^2) #offset from vertex to cutting plane. Plane is necessary to prevent parts of the complement from showing up in the final CSG
    p₀ = frontvertex + d
    plane = Plane(0.0, 0.0, 1.0, 0.0, 0.0, p₀)
    gen = csgdifference(leaf(plane), temp)

    concave = gen()
    intsct = surfaceintersection(concave, r)
    println(intsct)
end

# function testsinglepossurface()
#     backradius = 60.0
#     frontvertex = 0.0
#     thickness = 10.0
#     semidiameter = 9.0
#     r = Ray([5.000000042187876, 5.000000042187876, 0.6143579891640621], [0.042187876699796796, 0.042187876699796796, -0.9982185963600987])

#     temp = leaf(Sphere(abs(backradius)), OpticSim.translation(0.0, 0.0, (frontvertex - thickness) - backradius))
#     d = backradius - sqrt(backradius^2 - semidiameter^2)
#     p₀ = frontvertex - thickness - d
#     plane = Plane(0.0, 0.0, -1.0, 0.0, 0.0, p₀)
#     # gen = csgintersection(leaf(plane),csgcomplement(temp))
#     gen = csgcomplement(temp)
#     concave = gen()
#     intsct = closestintersection(surfaceintersection(concave, r))
#     println(intsct)
# end

# function testnegcurveintersection()
#     gen = OpticSim.csgcomplement(leaf(Sphere(60.0), OpticSim.translation(0.0, 0.0, 60.0)))
#     concave = gen()

#     r = Ray([0.0, 0.0, 1.0], [0.0, 0.0, -1.0])
#     intsct = surfaceintersection(concave, r)

#     println(intsct)
# end

function testdoubleconvex()
    system = Examples.doubleconvex()
    green = 500 * Unitful.u"nm"
    r = Ray([0.0, 0.0, 10.0], [0.0, 0.0, -1.0])

    res = OpticSim.trace(system, r, green)
    println(res)
end

function testcooketriplet()
    system = Examples.cooketriplet()

    green = 500 * Unitful.u"nm"
    temperature = 20 * Unitful.u"°C"

    r = Ray([0.0, 0.0, 1.0], [0.0, 0.0, -1.0])

    res = OpticSim.trace(system, r, green, temperature)
end

function testopticalsystem(lenssystem::T) where {T<:OpticalSystem}
    # lenssystem = Examples.cooketriplet()
    # lenssystem = Examples.planoconvex()
    lenssystem = Examples.doubleconvex()
    green = 500 * Unitful.u"nm"

    rad = OpticSim.semidiameter(lenssystem)
    points = Array{SVector{3,Float64},1}(undef, 0)

    count = 0

    for x in (-rad):0.5:rad, y in (-rad):0.5:rad
        r = Ray([x, y, 1.0], [0.0, 0.0, -1.0])
        if mod(count, 100) == 0
            println(count)
        end

        count += 1
        res = OpticSim.trace(lenssystem, r, green)

        if !(nothing === res)
            push!(points, point(res.intersection))
        end
    end

    twodpoints = map(x -> (x[1], x[2]), points)

    # scatter(twodpoints, markersize = 3)
end

function testelementsurfaceintersection()
    elt = OpticSim.SphericalLens(OpticSim.GlassCat.SCHOTT.BK6, 0.0, 60.0, Inf64, 5.0, 9.0)
    r = Ray([0.0, 0.0, 10.0], [0.0, 0.0, -1.0])
    res = surfaceintersection(elt, r)
    println(res)
end

function testsurfaceintersection()
    # hit from inside
    surf = TestData.beziersurface()
    accelsurf = AcceleratedParametricSurface(surf)
    TOLERANCE = 1e-9

    r = Ray([0.5, 0.5, 0.2], [0.0, 1.0, 0.0])
    println("1")
    res = surfaceintersection(accelsurf, r)
    lower(res) isa RayOrigin && isapprox(point(upper(res)), [0.5, 0.8996900152986361, 0.2], atol = TOLERANCE)

    # hit from outside
    r = Ray([0.0, 0.5, 0.2], [1.0, 0.0, -1.0])
    println("2")
    res = surfaceintersection(accelsurf, r)
    isapprox(point(lower(res)), [0.06398711204047353, 0.5, 0.13601288795952649], atol = TOLERANCE) && upper(res) isa Infinity

    # two hits from outside
    r = Ray([0.0, 0.5, 0.25], [1.0, 0.0, 0.0])
    println("3")
    res = surfaceintersection(accelsurf, r)
    isapprox(point(lower(res)), [0.12607777053030267, 0.5, 0.25], atol = TOLERANCE) && isapprox(point(upper(res)), [0.8705887077060419, 0.5, 0.25], atol = TOLERANCE)

    surf = TestData.wavybeziersurface()
    accelsurf = AcceleratedParametricSurface(surf)

    # 3 hit starting outside
    r = Ray([0.5, 0.0, 0.0], [0.0, 1.0, 0.0])
    println("4")
    res = surfaceintersection(accelsurf, r)
    isa(res, DisjointUnion) && length(res) == 2 && isapprox(point(lower(res[1])), [0.5, 0.10013694786182059, 0.0], atol = TOLERANCE) && isapprox(point(upper(res[1])), [0.5, 0.49625, 0.0], atol = TOLERANCE) && isapprox(point(lower(res[2])), [0.5, 0.8971357794109067, 0.0], atol = TOLERANCE) && (upper(res[2]) isa Infinity)

    # 3 hit starting inside
    #     r = Ray([0.5, 1.0, 0.0], [0.0, -1.0, 0.0])
    #     println("5")
    #     res = surfaceintersection(accelsurf, r)
    #     isa(res, DisjointUnion) && length(res) == 2 && (lower(res[1]) isa RayOrigin) && isapprox(point(upper(res[1])), [0.5, 0.8971357794109067, 0.0], atol = TOLERANCE) && isapprox(point(lower(res[2])), [0.5, 0.49625, 0.0], atol = TOLERANCE) && isapprox(point(upper(res[2])), [0.5, 0.10013694786182059, 0.0], atol = TOLERANCE)
    # println("after 5")

    surf = TestData.verywavybeziersurface()
    accelsurf = AcceleratedParametricSurface(surf)

    # five hits starting outside
    r = Ray([0.9, 0.0, -0.3], [0.0, 1.0, 0.7])
    println("6")
    res = surfaceintersection(accelsurf, r)
    a = isapprox(point(lower(res[1])), [0.9, 0.03172286522032046, -0.2777939943457758], atol = TOLERANCE) && isapprox(point(upper(res[1])), [0.9, 0.1733979947040411, -0.17862140370717122], atol = TOLERANCE)
    b = isapprox(point(lower(res[2])), [0.9, 0.5335760974594397, 0.07350326822160776], atol = TOLERANCE) && isapprox(point(upper(res[2])), [0.9, 0.7767707607392784, 0.24373953251749486], atol = TOLERANCE)
    c = isapprox(point(lower(res[3])), [0.9, 0.9830891958374246, 0.3881624370861975], atol = TOLERANCE) && (upper(res[3]) isa Infinity)
    isa(res, DisjointUnion) && length(res) == 3 && a && b && c

    # five hits starting inside
    r = Ray([0.9, 1.0, 0.4], [0.0, -1.0, -0.7])
    println("7")
    res = surfaceintersection(accelsurf, r)
    a = (lower(res[1]) isa RayOrigin) && isapprox(point(upper(res[1])), [0.9, 0.9830891958374246, 0.3881624370861975], atol = TOLERANCE)
    b = isapprox(point(lower(res[2])), [0.9, 0.7767707607392784, 0.24373953251749486], atol = TOLERANCE) && isapprox(point(upper(res[2])), [0.9, 0.5335760974594397, 0.07350326822160776], atol = TOLERANCE)
    c = isapprox(point(lower(res[3])), [0.9, 0.17339799470404108, -0.1786214037071712], atol = TOLERANCE) && isapprox(point(upper(res[3])), [0.9, 0.03172286522032046, -0.27779399434577573], atol = TOLERANCE)
    isa(res, DisjointUnion) && length(res) == 3 && a && b && c

    # 4 hits starting inside
    r = Ray([0.1, 0.0, -0.3], [0.0, 1.0, 0.7])
    println("8")
    res = surfaceintersection(accelsurf, r)
    a = (lower(res[1]) isa RayOrigin) && isapprox(point(upper(res[1])), [0.1, 0.2851860296285551, -0.10036977926001144], atol = TOLERANCE)
    b = isapprox(point(lower(res[2])), [0.1, 0.5166793625025807, 0.06167555375180668], atol = TOLERANCE) && isapprox(point(upper(res[2])), [0.1, 0.7770862508789854, 0.24396037561528983], atol = TOLERANCE)
    c = isapprox(point(lower(res[3])), [0.1, 0.98308919558696, 0.3881624369108719], atol = TOLERANCE) && (upper(res[3]) isa Infinity)
    isa(res, DisjointUnion) && length(res) == 3 && a && b && c

    # 4 hits starting outside
    r = Ray([0.9, 0.9, 0.4], [0.0, -0.9, -0.7])
    println("9")
    res = surfaceintersection(accelsurf, r)
    a = isapprox(point(lower(res[1])), [0.9, 0.736072142615238, 0.2725005553674076], atol = TOLERANCE) && isapprox(point(upper(res[1])), [0.9, 0.567439326091764, 0.141341698071372], atol = TOLERANCE)
    b = isapprox(point(lower(res[2])), [0.9, 0.16601081959179267, -0.1708804736508277], atol = TOLERANCE) && isapprox(point(upper(res[2])), [0.9, 0.032434058775915924, -0.274773509840954], atol = TOLERANCE)
    isa(res, DisjointUnion) && length(res) == 2 && a && b

end

function errortest()
    surf = TestData.verywavybeziersurface()
    accelsurf = AcceleratedParametricSurface(surf)

    TOLERANCE = 1e-9
    r = Ray([0.5, 1.0, 0.0], [0.0, -1.0, 0.0])

    res = surfaceintersection(accelsurf, r)
    isa(res, DisjointUnion) && length(res) == 2 && (lower(res[1]) isa RayOrigin) && isapprox(point(upper(res[1])), [0.5, 0.8971357794109067, 0.0], atol = TOLERANCE) && isapprox(point(lower(res[2])), [0.5, 0.49625, 0.0], atol = TOLERANCE) && isapprox(point(upper(res[2])), [0.5, 0.10013694786182059, 0.0], atol = TOLERANCE)
end

function testintervalunion()
    alphasequal(int1, int2) = α(lower(int1)) == α(lower(int2)) && α(upper(int1)) == α(upper(int2))
    intersectionat(a) = TestData.intersectionat(a)
    a = Interval(intersectionat(0.0), intersectionat(0.2))
    aa = Interval(intersectionat(0.1), intersectionat(0.3))

    b = DisjointUnion(aa, Interval(intersectionat(0.5), intersectionat(0.6)))
    res = intervalunion(a, b)

    @assert res == intervalunion(b, a)
    println("alphas equal $(alphasequal(res[1],Interval(intersectionat(0.0), intersectionat(0.3)))) ")

    # @assert res[1] == Interval(intersectionat(0.0), intersectionat(0.3)) this equality test does not work
    @assert res[2] == b[2]
end

function fixequalitybugintests()
    #conclusion: putting BezierSurface as a field in Intersection causes simple value semantics == to fail. Maybe BezierSurface field becomes a ref and is heap allocated?
    a = DisjointUnion(Interval(TestData.intersectionat(0.2), TestData.intersectionat(0.5)), Interval(TestData.intersectionat(0.7), TestData.intersectionat(0.8)))
    b = DisjointUnion(Interval(TestData.intersectionat(0.1), TestData.intersectionat(0.3)), Interval(TestData.intersectionat(0.4), TestData.intersectionat(0.6)))
    res = intervalunion(a, b)

    # println("a∪b $(intervalunion(a,b))")
    # println("identity $(intervalunion(a,b) == intervalunion(a,b))")
    # println("b∪a $(intervalunion(b,a))")
    @assert res == intervalunion(b, a)
    temp = Interval(TestData.intersectionat(0.1), TestData.intersectionat(0.6))
    println(temp)
    println(TestData.intersectionat(0.1))
    @assert res[1] == temp

    @assert res[2] == a[2]

end

function testcylinder()
    sph = Sphere(1.0)
    cyl = Cylinder(0.2)

    csg = csgdifference(leaf(sph), leaf(cyl))
    Vis.draw(csg)
end

function testtransmission()
    lens = TestData.planoplano()
    green = 500 * Unitful.u"nm"

    # result = trace(lens,r,green)
    # println(result)

    r = Ray([10.0, 0.0, 0.0], [-1.0, 0.0, -1.0])
    result = trace(lens, r, green)
    println(result)

    lenselement = lens.system.assembly.elements[1].objecttree
    println("Direct result of surface intersection on csg tree $(OpticSim.evalcsg(lenselement,r))")
    # pln1 = Plane(0.0,0.0,1.0,0.0,0.0,0.0)
    # pln2 = Plane(0.0,0.0,-1.0,0.0,0.0,-10.0)

    pln1 = Rectangle([0.0, 0.0, 1.0], [0.0, 0.0, 0.0], 9.0, 9.0)
    pln2 = Rectangle([0.0, 0.0, -1.0], [0.0, 0.0, -10.0], 9.0, 9.0)

    println("intersection top rectangle $(surfaceintersection(pln1,r)) \n\n\n intersection bottom rectangle $(surfaceintersection(pln2,r))")

    cyl = Cylinder(9.0, 10.0)
    csg = csgintersection(leaf(pln1), csgintersection(leaf(pln2), leaf(cyl)))
    intsct = OpticSim.evalcsg(csg(), r)
    println("rectangle top and bottom CSG intsct $intsct")
    # result = surfaceintersection(cyl, r)
    # println(result)
end

function testgenerateray()
    gen = OpticSim.RayGeneratorUniform(10, 2.0, 10.0)

    for i in gen
        println(i)
    end
end

function testrectangle()
    r = Rectangle([0.0, 0.0, 1.0], [0.0, 0.0, 0.0], 1.0, 1.0)
    ry = Ray([2.0, 0.0, 1.0], [0.0, 0.0, -1.0])
    println(surfaceintersection(r, ry))

end

function randunit()
    let v = rand(SVector{3,Float64})
        while (norm(v) > 1.0)
            v = rand(SVector{3,Float64})
        end
        return normalize(v)
    end
end

function makeeye()
    eye = CSGOpticalSystem(OpticSim.ModelEye(), Rectangle(10.0, 10.0, SVector{3,Float64}(0.0, 0.0, 1.0), SVector{3,Float64}(0.0, 0.0, -26.0)))
    # Vis.drawtracerays(eye)
    ray = OpticalRay([0.0, 0.0, 10.0], [0.0, 0.001, -1.0], 1.0, 0.55)

    return ray, eye
end

function dorays(ray, eye)
    # Threads.@threads for i in 1:1000000
    for i in 1:1000
        trace(eye, ray)
    end
end

function testnearestsquareroots()
    for i in 1:10:1000
        a, b = OpticSim.nearestsqrts(i)
        approx = a * b
        println(approx - i)
    end
end

function hierarchicalimage()
    a = OpticSim.HierarchicalImage{Float32}(100, 200)

    for i in CartesianIndices(a)
        a[i] = 2.0
    end
    return a
end

struct SpecialNumber{T<:Real} <: Real
    val::T

    function SpecialNumber(num::T) where {T<:Real}
        if isnan(num)
            println("here")
        else
            return new{T}(num)
        end
    end
end

# const operators = (:+,:-,:*,:/,:sqrt,:cos,:sin,:exp,:ln)
for op in (:sin, :cos, :tan, :log, :exp, :sqrt, :-)
    eval(quote
        Base.$op(a::SpecialNumber) = SpecialNumber($op(a.val))
    end)
end

for op in (:+, :-, :*, :/)
    eval(quote
        Base.$op(a::SpecialNumber, b::SpecialNumber) = SpecialNumber($op(a.val, b.val))
    end)
end

Base.promote_rule(::Type{T}, ::Type{SpecialNumber{T1}}) where {T<:Real,T1<:Real} = SpecialNumber{promote_type(T, T1)}

function singleemitter(numrays)
    T = Float64
    temp = SVector{3,T}(0.0, 0.0, 0.0)
    # λ = T(0)
    a = Source(UniformSpectrum{T}(), PointOrigin{T}(SVector{3,T}(0.0, 0.0, 0.0)), ConeDistribution(π / 5, numrays), Lambertian{T}())

    # function tracerays()
    #   temp = T(0)
    #         for ray in a
    #             temp = ray.power
    #         end
    #         return temp
    # end

    for ray in a
        temp = OpticSim.direction(a, 1)
        #  temp  = OpticSim.spectrumsample(a)
        # temp = OpticSim.spectrumpower(a.spectrum,.6)
    end

    return temp
end
export singleemitter

printrays(a::Source) =
    for ray in a
        println(ray)
    end

function diagdistributions(numsamples)
    temp = 0.0

    for i in 1:numsamples
        temp = rand(Distributions.Uniform{Float64}(400.0, 680.0))
    end
    return temp
end
export diagdistributions

end # module Junk


text:
# ## Adding LaTeX to plots (2)

#md # ![FILE_NAME.png](images/FILE_NAME.png)

# For pyplot backend

## ref https://jp.mathworks.com/matlabcentral/fileexchange/35237-matlab-plot-gallery-adding-latex-to-plots?focused=6792946&tab=example

using LaTeXStrings, Measures
pyplot()
Plots.PyPlot.rc("text", usetex ="true")

fib = zeros(12);
for i = 1:12
    fib[i] = (((1+sqrt(5))/2)^i - ((1-sqrt(5))/2)^i)/sqrt(5);
end

plot(fib,
    marker=:circle,
    xlabel=L"n",
    ylabel=L"F_n",
    annotation=(5, 100, L"F_n = \frac{1}{\sqrt{5}} \left[\left( \frac{1+\sqrt{5}}{2} \right)^n - \left( \frac{1-\sqrt{5}}{2} \right)^n \right]")
)

ext: jl


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module TikzUtils

export add_to_document!, use_geometry_package!, use_aircraftshapes_package!, wrap_tikzpicture!,
    lualatex_compile

using TikzPictures
using RLESUtils, FileUtils

function add_to_document!(d::TikzDocument, tp::TikzPicture, cap::AbstractString)
    push!(d, tp, caption=cap)
end

function add_to_document!(d::TikzDocument,tps::Array{TikzPicture,1},cap::AbstractString)
    [add_to_document!(d,tp,cap) for tp in tps]
end

function add_to_document!{T<:AbstractString}(d::TikzDocument,tpsTups::Array{Tuple{TikzPicture,T},1})
    [add_to_document!(d,tps...) for tps in tpsTups]
end

function add_to_document!{T<:AbstractString}(d::TikzDocument,
    tpsTups::Array{Tuple{Vector{TikzPicture},T},1})
    [add_to_document!(d,tps...) for tps in tpsTups]
end

function use_geometry_package!(p::TikzPicture; margin::Float64 = 0.5, landscape::Bool = false)
    # Adds margin enforcement to pdf file
    # Adds capability for landscape
    # margin in inches
    # landscape false is portrait
    orientation_str = landscape ? ",landscape" : ""
    #prepend geometry package
    p.preamble = "\\usepackage[margin=$(margin)in" * orientation_str * "]{geometry}\n" * p.preamble
    p
end

function use_aircraftshapes_package!(p::TikzPicture)
    #prepend aircraftshape package
    p.preamble = "\\usepackage{aircraftshapes}\n" * p.preamble
    p
end

function wrap_tikzpicture!(texfile::AbstractString, pre::AbstractString, post::AbstractString)
    foldername = dirname(texfile)
    if isempty(foldername)
        foldername = "."
    end
    begintikz = "\\begin{tikzpicture}"
    endtikz = "\\end{tikzpicture}"
    replace_text(texfile, begintikz, pre * begintikz, foldername, fopen=open, inplace=true)
    replace_text(texfile, endtikz, endtikz * post, foldername, fopen=open, inplace=true)
end

function lualatex_compile(texfile::AbstractString)
    foldername = dirname(texfile)
    if isempty(foldername)
        foldername = "."
    end
    cmd = `lualatex --output-directory=$(foldername) $(texfile)`
    success(cmd)
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
module MultiCellArraysTests

using Test
using Gridap
using Gridap.CellValuesGallery
using Gridap.CachedArrays

l = 10
sv = 1.0
sa = [sv, sv, sv]
sca = ConstantCellArray(sa,l)

sb = [2*sv, sv, 3*sv]
scb = ConstantCellArray(sb,l)

mca = MultiCellArray([sca,scb],[(1,),(3,)])

@test length(mca) == l

i_to_field = mca.fieldids

for ma in mca
  for (i,a) in enumerate(ma)
    ifield, = i_to_field[i]
    @assert isa(a,CachedArray)
  end
  @assert length(ma) == 2
  @assert isa(ma[1],CachedArray)
  @assert isa(ma[2],CachedArray)
end

mca2 = MultiCellVector([sca,scb],[(1,),(3,)])

sa = [sv sv; sv sv]
sca = ConstantCellArray(sa,l)

sb = [2*sv sv; 3*sv sv]
scb = ConstantCellArray(sb,l)

mca2 = MultiCellMatrix([sca,scb],[(1,5),(3,2)])

ids = [1,3,4,2]
cellids = CellValueFromArray(ids)
mca3 = reindex(mca,cellids)

for (ca3,ca) in zip(mca3.blocks,mca.blocks)
  @test ca3 == reindex(ca,cellids)
end

end # module MultiCellArraysTests


text:
module test_massComputation

using  Modia3D
import Modia3D.ModiaMath
using Modia3D.StaticArrays
using Modia3D.Unitful

vmat1 = Modia3D.Material(color="LightBlue", transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

rCM1 = zeros(3)
rCM1[2] = 9.0
rCM2 = zeros(3)
rCM2[1] = 5.0
rCM2[3] = 2.0
I1 = fill(2.0,(3,3))
I2 = fill(3.0,(3,3))
massProp1 = Modia3D.MassProperties(15.0, rCM1, I1)
massProp2 = Modia3D.MassProperties(5.0, rCM2, I2)

rCM3 = zeros(3)
rCM3[2] = 4.0
rCM4 = zeros(3)
rCM4[1] = 3.0
rCM4[3] = 9.0
I3 = fill(45.0,(3,3))
I4 = fill(100.0,(3,3))
massProp3 = Modia3D.MassProperties(30.0, rCM3, I3)
massProp4 = Modia3D.MassProperties(60.0, rCM4, I4)

@signal Sine(;A=1.0, freqHz = 1.0) begin
   y = ModiaMath.RealScalar(causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::Sine, sim::ModiaMath.SimulationState)
    signal.y.value  = signal.A*sin(2*pi*signal.freqHz*sim.time)
end

@assembly DoublePendulum(;length=1.0, width=0.2) begin
   world = Modia3D.Object3D(Modia3D.CoordinateSystem(0.0))

   part1 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBeam(length,width,width), nothing, vmat1),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])
   part2 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBeam(length,width,width), nothing, vmat1),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])

   part3 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBeam(length,width,width), nothing, vmat1),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])
   part4 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBeam(length,width,width), nothing, vmat1),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])

                          #(Lx,Ly,Lz)

   cyl1    = Modia3D.Cylinder(width/2,1.2*width; material=vmat2)
   cyl2    = Modia3D.Cylinder(width/2,1.2*width; material=vmat2)

   cyl11   = Modia3D.Object3D(part1.frames[1], cyl1)
   cyl21   = Modia3D.Object3D(part4.frames[2], cyl1)

   rev1 = Modia3D.Revolute(world          , part1.frames[1])
   Modia3D.connect(part1.frames[2], part2.frames[1])

   Modia3D.connect(part2.frames[2], part3.frames[1])
   Modia3D.connect(part3.frames[2], part4.frames[1])

   sig1 = Sine(A=1.5,freqHz=0.5)
   sig2 = Sine(A=2.5,freqHz=1.0)
   sine1     = Modia3D.SignalToFlangeAngle(sig1.y)
   Modia3D.connect(sine1,   rev1)
end

doublePendulum = DoublePendulum(sceneOptions = Modia3D.SceneOptions(visualizeFrames=false,
                                                                    defaultFrameLength=0.1,
                                                                    enableContactDetection=false))
model = Modia3D.SimulationModel( doublePendulum, analysis=ModiaMath.KinematicAnalysis )
result = ModiaMath.simulate!(model, stopTime=3.0)

ModiaMath.plot(result, ("rev1.phi", "rev2.phi") )

println("... success of test_massComputation.jl!")
end


text:
# Autogenerated wrapper script for Qt5Location_jll for i686-w64-mingw32-cxx11
export libqt5location, libqt5positioning, libqt5positioningquick

using Qt5Declarative_jll
JLLWrappers.@generate_wrapper_header("Qt5Location")
JLLWrappers.@declare_library_product(libqt5location, "Qt5Location.dll")
JLLWrappers.@declare_library_product(libqt5positioning, "Qt5Positioning.dll")
JLLWrappers.@declare_library_product(libqt5positioningquick, "Qt5PositioningQuick.dll")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Declarative_jll)
    JLLWrappers.@init_library_product(
        libqt5location,
        "bin\\Qt5Location.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5positioning,
        "bin\\Qt5Positioning.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5positioningquick,
        "bin\\Qt5PositioningQuick.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    batch_quadrature(
        fs::AbstractVector,
        ms::AbstractVector{<:Real},
        σs::AbstractVector{<:Real},
        num_points::Integer,
    )

Approximate the integrals
```julia
∫ fs[n](x) pdf(Normal(ms[n], σs[n]), x) dx
```
for all `n` in `eachindex(fs)` using Gauss-Hermite quadrature with `num_points`.
"""
function batch_quadrature(
    fs::AbstractVector,
    ms::AbstractVector{<:Real},
    σs::AbstractVector{<:Real},
    num_points::Integer,
)
    # Check that as many bounds are provided as we have functions to integrate.
    length(fs) == length(ms) || throw(error("length(fs) != length(ms)"))
    length(fs) == length(σs) || throw(error("length(fs) != length(σs)"))

    # Construct the quadrature points.
    xs, ws = gausshermite(num_points)

    # Compute the integral.
    return map((f, m, σ) -> _gauss_hermite_quadrature(f, m, σ, xs, ws), fs, ms, σs)
end

Zygote.@nograd gausshermite

# Internal method. Assumes that the gradient w.r.t. xs and ws is never needed, so avoids
# computing it and returns nothing. This is potentially not what you want in general.
function _gauss_hermite_quadrature(f, m::Real, σ::Real, xs, ws)
    t(x, m, σ) = m + sqrt(2) * σ * x
    I = ws[1] * f(t(xs[1], m, σ))
    for j in 2:length(xs)
        I += ws[j] * f(t(xs[j], m, σ))
    end
    return I / sqrt(π)
end

function Zygote._pullback(
    ctx::Zygote.AContext, ::typeof(_gauss_hermite_quadrature), f, m::Real, σ::Real, xs, ws,
)
    function _gauss_hermite_quadrature_pullback(Δ::Real)
        g(f, x, w, m, σ) = w * f(m + sqrt(2) * σ * x)

        _, pb = Zygote._pullback(ctx, g, f, xs[1], ws[1], m, σ)
        _, Δf, _, _, Δm, Δσ = pb(Δ / sqrt(π))
        for j in 2:length(xs)
            _, pb = Zygote._pullback(ctx, g, f, xs[j], ws[j], m, σ)
            _, Δf_, _, _, Δm_, Δσ_ = pb(Δ / sqrt(π))
            Δf = Zygote.accum(Δf, Δf_)
            Δm = Zygote.accum(Δm, Δm_)
            Δσ = Zygote.accum(Δσ, Δσ_)
        end
        return nothing, Δf, Δm, Δσ, nothing, nothing
    end
    return _gauss_hermite_quadrature(f, m, σ, xs, ws), _gauss_hermite_quadrature_pullback
end


text:
module day13

include("inputs.jl")

ex = inputs.ex_d13()
input = inputs.input_d13()

function part1(input)
    depart = parse(Int, input[1])
    buses = [parse(Int, bus) for bus in split(input[2], ",") if bus !="x"]
    waits = [bus - mod(depart, bus) for bus in buses]
    idx = argmin(waits)
    return waits[idx]*buses[idx]
end

@assert part1(ex) == 295
@assert part1(input) == 246
@info "Day 13, Part 1 answer: $(part1(input))"

# Part 2
# A hint pointed me to the Chinese Remainder Theorem. Alg from Rosetta Code

function chineseremainder(n::Array, a::Array)
    Π = prod(n)
    mod(sum(ai * invmod(Π ÷ ni, ni) * (Π ÷ ni) for (ni, ai) in zip(n, a)), Π)
end

function part2(input)
    en_buses = [(i, parse(Int, bus)) for (i, bus) in enumerate(split(input[2], ",")) if bus !="x"]
    offsets = first.(en_buses)
    bus_ids = last.(en_buses)
    m = maximum(offsets)
    return chineseremainder(bus_ids, m .- offsets) - m + 1
end

@assert part2(ex) == 1068781
@assert part2(["None", "67,7,59,61"]) == 754018
@assert part2(["None", "67,x,7,59,61"]) == 779210
@assert part2(["None", "67,7,x,59,61"]) == 1261476
@assert part2(["None", "1789,37,47,1889"]) == 1202161486
@assert part2(input) == 939490236001473
@info "Day 13, Part 2 answer: $(part2(input))"

end


text:
export process_results, process_nsga2_results, process_ucea_results

using DataFrames
using Plots
using ImageFiltering
using OffsetArrays
using BenchmarkTools
using Statistics
using CSV
using TimerOutputs

BASELINES = Dict(
    "assault"=>Dict("MTCGP"=>(890.4, 255), "A3C LSTM"=>14497.9),
    "asteroids"=>Dict("MTCGP"=>(9412,1818)),
    "boxing"=>Dict("MTCGP"=>(38.4,4), "Dueling"=>77.3),
    "breakout"=>Dict("MTCGP"=>(13.2,2), "A3C LSTM"=>766.8),
    "defender"=>Dict("MTCGP"=>(993010,2739)),
    "freeway"=>Dict("MTCGP"=>(28.2,0), "HyperNEAT"=>29),
    "frostbite"=>Dict("MTCGP"=>(782,795), "TPG"=>8144.4),
    "gravitar"=>Dict("MTCGP"=>(2350,50)),
    "private_eye"=>Dict("MTCGP"=>(12702.2,4337), "TPG"=>15028.3),
    "pong"=>Dict("MTCGP"=>(20,0)),
    "riverraid"=>Dict("MTCGP"=>(2914,90), "Prioritized"=>18184.4),
    "solaris"=>Dict("MTCGP"=>(8324,2250)),
    "space_invaders"=>Dict("MTCGP"=>(1001,25), "A3C LSTM"=>23846)
)

gamename_from_romname(rn::String) = titlecase(replace(rn, "_"=>" "))

"""
    function time_monocgp_ms(
        reducer::Reducer,
        cont::CGPInd,
        game_name::String
    )

Time forward pass of complete architecture.
Return time in ms.
"""
function time_monocgp_ms(
    reducer::Reducer,
    cont::CGPInd,
    game_name::String,
    grayscale::Bool,
    downscale::Bool
)
    # Generate input
    g = Game(game_name, 0)
    inp = get_state(g, grayscale, downscale)
    close!(g)
    # Pre-compile
    process(reducer, cont, inp)
    n_iter = 1000
    # Time
    redu_time = 0.0
    flat_time = 0.0
    cont_time = 0.0
    for _ in 1:n_iter
        redu_time += @elapsed features = reducer.reduct(inp, reducer.parameters)
        flat_time += @elapsed features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))
        cont_time += @elapsed cont_out = CartesianGeneticProgramming.process(cont, features_flatten)
    end
    to_ms = 1000 / n_iter
    redu_time *= to_ms
    flat_time *= to_ms
    cont_time *= to_ms
    redu_time, flat_time, cont_time
end

"""
    function time_dualcgp_ms(
        enco::CGPInd,
        reducer::Reducer,
        cont::CGPInd,
        game_name::String
    )

Time forward pass of complete architecture.
Return time in ms.
"""
function time_dualcgp_ms(
    enco::CGPInd,
    reducer::Reducer,
    cont::CGPInd,
    game_name::String,
    grayscale::Bool,
    downscale::Bool
)
    # Generate input
    g = Game(game_name, 0)
    inp = get_state(g, grayscale, downscale)
    close!(g)
    # Pre-compile
    process(enco, reducer, cont, inp)
    n_iter = 1000
    # Time
    enco_time = 0.0
    redu_time = 0.0
    flat_time = 0.0
    cont_time = 0.0
    for _ in 1:n_iter
        enco_time += @elapsed enco_out = CartesianGeneticProgramming.process(enco, inp)
        redu_time += @elapsed features = reducer.reduct(enco_out, reducer.parameters)
        flat_time += @elapsed features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))
        cont_time += @elapsed cont_out = CartesianGeneticProgramming.process(cont, features_flatten)
    end
    to_ms = 1000 / n_iter
    enco_time *= to_ms
    redu_time *= to_ms
    flat_time *= to_ms
    cont_time *= to_ms
    enco_time, redu_time, flat_time, cont_time
end

function add_baselines(graphs::Vector{Plots.Plot{Plots.GRBackend}}, game::String)
    if haskey(BASELINES, game)
        baselines = BASELINES[game]
        for g in graphs
            for k in keys(baselines)
                if k == "MTCGP"
                    score = baselines[k][1]
                    var = baselines[k][2]
                    hline!(g, [score], ribbon=([var]), color=:gray, label=false)
                else
                    score = baselines[k]
                end
                hline!(g, [score], linestyle=:dash, color=:black, label=false)
                annotate!(0, score, text(k, 7, :bottom, :left), :black)
            end
        end
    end
end

function getlabel(labels::Vector{String}, i::Int64, reducer_type::String)
    if i > length(labels)
        return string(reducer_type)
    else
        return labels[i]
    end
end

function getcolor(colors::Vector{Symbol}, i::Int64, reducer_type::String)
    if i > length(colors)
        if reducer_type == "pooling"
            return :skyblue3
        elseif reducer_type == "centroid"
            return :chocolate1
        else
            return :black
        end
    else
        return colors[i]
    end
end

"""
    plot_cibest(
        exp_dirs::Vector{String},
        games::Vector{String},
        gamename::String,
        labels::Vector{String},
        colors::Vector{Symbol};
        ma::Int64=1
    )

Plot the fitness of the best individual for each generation with
confidence intervals across all experiments using the exact same
configuration dictionary.
"""
function plot_cibest(
    exp_dirs::Vector{String},
    games::Vector{String},
    gamename::String,
    labels::Vector{String},
    colors::Vector{Symbol};
    ma::Int64=1
)
    # 1. Link cfgs and indexes to categories
    n_cat = 0
    cfg2cat = Dict{Dict{Any,Any},Int64}() # Link cfg to category
    ind2cat = Dict{Int64,Int64}() # Link index to category
    for i in eachindex(exp_dirs)
        exp_dir = exp_dirs[i]
        cfg = cfg_from_exp_dir(exp_dir)
        if cfg in keys(cfg2cat) # Link index to existing category
            cat = cfg2cat[cfg]
            ind2cat[i] = cat
        else # Add a new category
            n_cat += 1
            cfg2cat[cfg] = n_cat
            ind2cat[i] = n_cat
        end
    end
    # 2. Fetch logs
    cat2logs = Dict{Int64,Vector{Vector{Float64}}}()
    for i in eachindex(exp_dirs)
        exp_dir = exp_dirs[i]
        log = log_from_exp_dir(exp_dir)
        kernel = OffsetArray(fill(1/(2*ma+1), 2*ma+1), -ma:ma)
        best_i = (ma == 1) ? log.best : imfilter(log.best, kernel)
        cat = ind2cat[i]
        if cat ∈ keys(cat2logs)
            push!(cat2logs[cat], best_i)
        else
            cat2logs[cat] = [best_i]
        end
    end
    # 3. Compute mean and confidence intervals
    cat2mean = Dict{Int64,Vector{Float64}}()
    cat2ci = Dict{Int64,Vector{Float64}}()
    for cat in 1:n_cat
        cat2mean[cat] = mean(cat2logs[cat])
        cat2ci[cat] = std(cat2logs[cat])
    end
    # 4. Fetch labels and colors
    @assert length(labels)==n_cat "please provide exactly one label per category"
    @assert length(colors)==n_cat "please provide exactly one color per category"
    cat2lab = Dict{Int64,String}()
    cat2col = Dict{Int64,Symbol}()
    for cat in 1:n_cat
        cat2lab[cat] = labels[cat]
        cat2col[cat] = colors[cat]
        cfg = Dict{Any,Any}()
        for k in keys(cfg2cat)
            if cfg2cat[k] == cat
                cfg = k
            end
        end
        println()
        println("-"^100)
        println("Printing map:")
        println("Number of detected experiments: ", length(cat2logs[cat]))
        println("Label: ", labels[cat])
        println("Color: ", colors[cat])
        println("Configuration: ", cfg)
        println("-"^100)
        println()
    end
    # 5. Plot
    xlabel = "Generation"
    ylabel = string("Best score ", gamename)
    plt = plot(xlabel=xlabel, ylabel=ylabel)
    for cat in 1:n_cat
        y = cat2mean[cat]
        ci = cat2ci[cat]
        lab = cat2lab[cat]
        col = cat2col[cat]
        plot!(plt, y, ribbon=ci, label=lab, color=col)
    end
    #=
    # Re-plot particular experiments
    for i in eachindex(exp_dirs)
        exp_dir = exp_dirs[i]
        cfg = cfg_from_exp_dir(exp_dir)
        log = log_from_exp_dir(exp_dir)
        kernel = OffsetArray(fill(1/(2*ma+1), 2*ma+1), -ma:ma)
        best = (ma == 1) ? log.best : imfilter(log.best, kernel)
        mean = (ma == 1) ? log.mean : imfilter(log.mean, kernel)
        std = (ma == 1) ? log.std : imfilter(log.std, kernel)
        save_gen = cfg["save_gen"]
        x = 1:save_gen:save_gen*length(best)
        plot!(plt, x, best)#, label=label_i, color=color_i)
    end
    =#
    plt
end

function set_graph_dir(
    exp_dirs::Vector{String},
    savedir_index::Int64,
    game::String
)
    exp_dir = exp_dirs[savedir_index]
    set_graph_dir(exp_dir, game)
end

function set_graph_dir(exp_dir::String, game::String)
    exp_name = string(basename(exp_dir)[1:10], "_", game)
    graph_dir = joinpath(dirname(dirname(exp_dir)), "graphs", exp_name)
    mkpath(graph_dir)
    graph_dir
end

function strvec2vec(x::Union{CSV.InlineString,CSV.String31,String})
    x = replace(x, "[" => "")
    x = replace(x, "]" => "")
    [parse(Float64, xi) for xi in split(x, ",")]
end

replace_nan(v::Vector{T}, val::T) where T = map(x -> isnan(x) ? val : x, v)

function retrieve_ucea_statistics(log::CSV.File)
    df = DataFrame(log)
    datadict = Dict()
    datadict["gen"] = Vector{Int64}()
    datadict["best_mean"] = Vector{Float64}()
    datadict["best_means_std"] = Vector{Float64}()
    datadict["best_means_best_score"] = Vector{Float64}()
    datadict["best_means_neval"] = Vector{Int64}()
    datadict["best_score"] = Vector{Float64}()
    datadict["best_scores_neval"] = Vector{Int64}()
    for row in eachrow(df)
        if row.gen_number ∉ datadict["gen"]
            push!(datadict["gen"], row.gen_number)
        end
    end
    for gen in datadict["gen"]
        df_gen = filter(row -> row.gen_number == gen, df)
        fitnesses = [strvec2vec(f) for f in df_gen.fitnesses]
        bests = [maximum(f) for f in fitnesses]
        means = [mean(f) for f in fitnesses]
        nevals = [length(f) for f in fitnesses]
        stds = replace_nan([std(f) for f in fitnesses], 0.0)
        best_mean_index = argmax(means)
        best_score_index = argmax(bests)
        push!(datadict["best_mean"], means[best_mean_index])
        push!(datadict["best_means_std"], stds[best_mean_index])
        push!(datadict["best_means_best_score"], bests[best_mean_index])
        push!(datadict["best_means_neval"], nevals[best_mean_index])
        push!(datadict["best_score"], maximum(bests))
        push!(datadict["best_scores_neval"], nevals[best_score_index])
    end
    datadict
end

function process_ucea_results(
    exp_dirs::Vector{String},
    games::Vector{String},
    colors::Vector{Symbol},
    labels::Vector{String};
    do_display::Bool=true,
    do_save::Bool=true,
    savedir_index::Int64=1
)
    @assert all([g == games[1] for g in games])
    game = games[1]
    gamename = gamename_from_romname(game)
    # Create plots
    xl = "Generation"
    yl_fit = string("Score ", gamename)
    yl_neval = string("Number of evaluations ", gamename)
    plt_fit = plot(ylabel=yl_fit, xlabel=xl, legend=:bottomright)
    plt_neval = plot(ylabel=yl_neval, xlabel=xl, legend=:topleft)
    # Fetch logs and plotpareto_gen
    for k in eachindex(exp_dirs)
        exp_dir = exp_dirs[k]
        cfg = cfg_from_exp_dir(exp_dir)
        log = log_from_exp_dir(exp_dir, log_file="logs/logs.csv",
            header=1, sep=";")
        datadict = retrieve_ucea_statistics(log)
        lb_mean = string(labels[k], " best mean")
        lb_best = string(labels[k], " best score of best mean individual")
        lb_best_ever = string(labels[k], " best score across all individuals")
        lb_best_means_neval = string(labels[k], " individual with best mean")
        lb_best_scores_neval = string(labels[k], " individual with best score")
        plot!(plt_fit, datadict["gen"], datadict["best_mean"],
            ribbon=datadict["best_means_std"],
            label=lb_mean, color=colors[k])
        plot!(plt_fit, datadict["gen"], datadict["best_means_best_score"],
            label=lb_best, color=colors[k], linewidth=2)
        plot!(plt_fit, datadict["gen"], datadict["best_score"],
            label=lb_best_ever, color=:red, linewidth=2)
        plot!(plt_neval, datadict["gen"], datadict["best_means_neval"],
            label=lb_best_means_neval, color=colors[k], linewidth=2)
        plot!(plt_neval, datadict["gen"], datadict["best_scores_neval"],
            label=lb_best_scores_neval, color=:red, linewidth=2)
    end
    # Display
    if do_display
        display(plt_fit)
        display(plt_neval)
    end
    # Save
    if do_save
        graph_dir = set_graph_dir(exp_dirs, savedir_index, game)
        savefig(plt_fit, joinpath(graph_dir, "ucea_fitness.png"))
        savefig(plt_neval, joinpath(graph_dir, "ucea_neval.png"))
    end
end

tofilename(s::String) = replace(lowercase(s), " " => "_")
pareto_key(i::Int64) = string("pareto_", i)
rawfit_max_key(i::Int64) = string("obj_", i, "_best_rawfit")
rawfit_mea_key(i::Int64) = string("obj_", i, "_mean_rawfit")
rawfit_std_key(i::Int64) = string("obj_", i, "_std_rawfit")
norfit_max_key(i::Int64) = string("obj_", i, "_best_norfit")
norfit_mea_key(i::Int64) = string("obj_", i, "_mean_norfit")
norfit_std_key(i::Int64) = string("obj_", i, "_std_norfit")

function retrieve_nsga2_statistics(
    log::CSV.File;
    pareto_gen::Vector{Int64}=Vector{Int64}()
)
    n_obj = length(strvec2vec(log[1].fitness))
    datadict = Dict()
    datadict["gen"] = Vector{Int64}()
    for i in 1:n_obj
        datadict[rawfit_max_key(i)] = Vector{Float64}()
        datadict[rawfit_mea_key(i)] = Vector{Float64}()
        datadict[rawfit_std_key(i)] = Vector{Float64}()
        datadict[norfit_max_key(i)] = Vector{Float64}()
        datadict[norfit_mea_key(i)] = Vector{Float64}()
        datadict[norfit_std_key(i)] = Vector{Float64}()
    end
    for row in log
        if row.gen_number ∉ datadict["gen"]
            push!(datadict["gen"], row.gen_number)
        end
    end
    for gen in datadict["gen"]
        rawfit = [strvec2vec(r.fitness) for r in log[log.gen_number .== gen]]
        norfit = [strvec2vec(r.normalized_fitness) for r in log[log.gen_number .== gen]]
        for i in 1:n_obj
            push!(datadict[rawfit_max_key(i)], maximum([f[i] for f in rawfit]))
            push!(datadict[rawfit_mea_key(i)], mean([f[i] for f in rawfit]))
            push!(datadict[rawfit_std_key(i)], std([f[i] for f in rawfit]))
            push!(datadict[norfit_max_key(i)], maximum([f[i] for f in norfit]))
            push!(datadict[norfit_mea_key(i)], mean([f[i] for f in norfit]))
            push!(datadict[norfit_std_key(i)], std([f[i] for f in norfit]))
        end
    end
    for gen in pareto_gen
        datadict[pareto_key(gen)] = [strvec2vec(r.normalized_fitness) for r in log[log.gen_number .== gen]]
    end
    datadict
end

function process_nsga2_results(
    exp_dirs::Vector{String},
    games::Vector{String},
    objectives_names::Vector{String},
    colors::Vector{Symbol},
    labels::Vector{String};
    pareto_gen::Vector{Int64}=Vector{Int64}(),
    pareto_xlim::Tuple=(0, 1),
    pareto_ylim::Tuple=(0, 1),
    do_display::Bool=true,
    do_save::Bool=true,
    savedir_index::Int64=1
)
    @assert all([g == games[1] for g in games])
    game = games[1]
    gamename = gamename_from_romname(game)
    # Create plots
    xl = "Generation"
    plt_rawfit = Vector{Plots.Plot}()
    plt_norfit = Vector{Plots.Plot}()
    plt_pareto = Vector{Plots.Plot}()
    for i in eachindex(objectives_names)
        yl_raw = string(objectives_names[i], " (", gamename, ")")
        yl_nor = string("Normalized ", objectives_names[i], " (", gamename, ")")
        push!(plt_rawfit, plot(ylabel=yl_raw, xlabel=xl))
        push!(plt_norfit, plot(ylabel=yl_nor, xlabel=xl))
    end
    for gen in pareto_gen
        title = string("Pareto front ", gamename, " (generation: ", gen, ")")
        xl_par = objectives_names[1]
        yl_par = objectives_names[2]
        push!(plt_pareto, plot(ylabel=yl_par, xlabel=xl_par, title=title))
    end
    # Fetch logs and plotpareto_gen
    for k in eachindex(exp_dirs)
        exp_dir = exp_dirs[k]
        cfg = cfg_from_exp_dir(exp_dir)
        log = log_from_exp_dir(exp_dir, log_file="logs/logs.csv",
            header=1, sep=";")
        datadict = retrieve_nsga2_statistics(log, pareto_gen=pareto_gen)
        for i in eachindex(objectives_names)
            best_lb = string(labels[k], " best")
            mean_lb = string(labels[k], " mean")
            plot!(plt_rawfit[i], datadict["gen"], datadict[rawfit_max_key(i)],
                label=best_lb, color=colors[k], linewidth=2)
            plot!(plt_rawfit[i], datadict["gen"], datadict[rawfit_mea_key(i)],
                ribbon=datadict[rawfit_std_key(i)], label=mean_lb,
                color=colors[k])
            plot!(plt_norfit[i], datadict["gen"], datadict[norfit_max_key(i)],
                label=best_lb, color=colors[k], linewidth=2)
            plot!(plt_norfit[i], datadict["gen"], datadict[norfit_mea_key(i)],
                ribbon=datadict[norfit_std_key(i)], label=mean_lb,
                color=colors[k])
        end
        for i in eachindex(pareto_gen)
            x = [nf[1] for nf in datadict[pareto_key(pareto_gen[i])]]
            y = [nf[2] for nf in datadict[pareto_key(pareto_gen[i])]]
            plot!(plt_pareto[i], x, y, seriestype = :scatter, xlims=pareto_xlim,
                ylims=pareto_ylim, label=labels[k], markercolor=colors[k])
        end
    end
    # Display
    if do_display
        for i in eachindex(plt_rawfit)
            display(plt_rawfit[i])
            #display(plt_norfit[i])
        end
        for i in eachindex(pareto_gen)
            display(plt_pareto[i])
        end
    end
    # Save
    if do_save
        graph_dir = set_graph_dir(exp_dirs, savedir_index, game)
        for i in eachindex(objectives_names)
            raw_name = tofilename(string(objectives_names[i], "_best.png"))
            nor_name = tofilename(string(objectives_names[i], "_normalized_best.png"))
            savefig(plt_rawfit[i], joinpath(graph_dir, raw_name))
            savefig(plt_norfit[i], joinpath(graph_dir, nor_name))
        end
        for i in eachindex(pareto_gen)
            par_name = tofilename(string("pareto_gen_", pareto_gen[i], ".png"))
            savefig(plt_pareto[i], joinpath(graph_dir, par_name))
        end
    end
end

"""
    function process_results(
        exp_dirs::Vector{String},
        games::Vector{String},
        dotime::Bool,
        dosave::Bool;
        savedir_index::Int64=1,
        ma::Int64=1,
        baselines::Bool=false,
        labels::Vector{String}=Vector{String}(),
        colors::Vector{Symbol}=Vector{Symbol}(),
        plotci::Bool=false,
        cilabels::Vector{String}=Vector{String}(),
        cicolors::Vector{Symbol}=Vector{Symbol}()
    )

Main results plot/print method.
"""
function process_results(
    exp_dirs::Vector{String},
    games::Vector{String},
    dotime::Bool,
    dosave::Bool;
    savedir_index::Int64=1,
    ma::Int64=1,
    baselines::Bool=false,
    labels::Vector{String}=Vector{String}(),
    colors::Vector{Symbol}=Vector{Symbol}(),
    plotci::Bool=false,
    cilabels::Vector{String}=Vector{String}(),
    cicolors::Vector{Symbol}=Vector{Symbol}()
)
    # Init graphs
    game = games[1]
    gamename = gamename_from_romname(game)
    xl = "Generation"
    yl_best = string("Best score ", gamename)
    yl_mean = string("Mean score ", gamename)
    plt_best = plot(ylabel=yl_best, xlabel=xl)
    plt_mean = plot(ylabel=yl_mean, xlabel=xl)
    if plotci
        plt_mbst = plot_cibest(exp_dirs, games, gamename, cilabels, cicolors, ma=ma)
    end

    for i in eachindex(exp_dirs)
        exp_dir = exp_dirs[i]
        cfg = cfg_from_exp_dir(exp_dir)
        log = log_from_exp_dir(exp_dir)

        # Plots
        reducer_type = cfg["reducer"]["type"]
        label_i = getlabel(labels, i, reducer_type)
        color_i = getcolor(colors, i, reducer_type)
        kernel = OffsetArray(fill(1/(2*ma+1), 2*ma+1), -ma:ma)
        best = (ma == 1) ? log.best : imfilter(log.best, kernel)
        mean = (ma == 1) ? log.mean : imfilter(log.mean, kernel)
        std = (ma == 1) ? log.std : imfilter(log.std, kernel)
        save_gen = cfg["save_gen"]
        x = 1:save_gen:save_gen*length(best)
        plot!(plt_best, x, best, label=label_i, color=color_i)
        plot!(plt_mean, x, mean, ribbon=std, label=label_i, color=color_i)

        # Append info to print
        p = []
        push!(p, ["Game", games[i]])
        push!(p, ["Number gen", length(log) * cfg["save_gen"]])
        push!(p, ["1st gen", ""])
        push!(p, ["  - best", log[1].best])
        push!(p, ["  - mean", log[1].mean])
        push!(p, ["  - std", log[1].std])
        push!(p, ["end gen", ""])
        push!(p, ["  - best", log[end].best])
        push!(p, ["  - mean", log[end].mean])
        push!(p, ["  - std", log[end].std])
        push!(p, ["Reducer", reducer_type])
        if cfg["reducer"]["type"] == "centroid"
            push!(p, ["  - n_centroids", cfg["reducer"]["n_centroids"]])
        else
            push!(p, ["  - features_size", cfg["reducer"]["features_size"]])
            push!(p, ["  - pooling_function", cfg["reducer"]["pooling_function"]])
        end

        # Timer
        if dotime
            push!(p, ["Forward pass timing (ms)", ""])
            is_dual = isfile(joinpath(exp_dirs[i], "logs/encoder.csv"))
            total_time = 0.0
            grayscale = cfg["grayscale"]
            downscale = cfg["downscale"]
            if is_dual
                enco, reducer, cont = get_last_dualcgp(exp_dirs[i], games[i], cfg)
                enco_time, redu_time, flat_time, cont_time = time_dualcgp_ms(
                    enco, reducer, cont, games[i], grayscale, downscale
                )
                push!(p, ["  - Encoder", enco_time])
                total_time += enco_time
            else
                reducer, cont = get_last_monocgp(exp_dirs[i], games[i], cfg)
                redu_time, flat_time, cont_time = time_monocgp_ms(
                    reducer, cont, games[i], grayscale, downscale
                )
            end
            total_time += redu_time
            total_time += flat_time
            total_time += cont_time
            push!(p, ["  - Reducer", redu_time])
            push!(p, ["  - Flattening", flat_time])
            push!(p, ["  - Controller", cont_time])
            push!(p, ["  - Total", total_time])
        end

        # Finish push and print
        l = maximum([length(k) for k in [pr[1] for pr in p]])
        println()
        for k in p
            println(string(k[1], " "^(l-length(k[1])), " : ", k[2]))
        end
    end
    if baselines
        add_baselines([plt_best, plt_mean], game)
    end
    display(plt_best)
    display(plt_mean)
    if plotci
        if baselines
            add_baselines([plt_mbst], game)
        end
        display(plt_mbst)
    end

    if dosave
        graph_dir = set_graph_dir(exp_dirs, savedir_index, game)
        savefig(plt_best, joinpath(graph_dir, "best.png"))
        savefig(plt_mean, joinpath(graph_dir, "mean.png"))
        if plotci
            savefig(plt_mbst, joinpath(graph_dir, "bestci.png"))
        end
    end
end


text:
# Common utilities

###########################################################
#
#   helper functions for dimension checking
#
###########################################################

function get_common_ncols(a::AbstractMatrix, b::AbstractMatrix)
    na = size(a, 2)
    size(b, 2) == na || throw(DimensionMismatch("The number of columns in a and b must match."))
    return na
end

function get_colwise_dims(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
    size(a) == size(b) || throw(DimensionMismatch("The sizes of a and b must match."))
    length(r) == size(a, 2) || throw(DimensionMismatch("Incorrect size of r."))
    return size(a)
end

function get_colwise_dims(r::AbstractArray, a::AbstractVector, b::AbstractMatrix)
    length(a) == size(b, 1) ||
        throw(DimensionMismatch("The length of a must match the number of rows in b."))
    length(r) == size(b, 2) || throw(DimensionMismatch("Incorrect size of r."))
    return size(b)
end

function get_colwise_dims(r::AbstractArray, a::AbstractMatrix, b::AbstractVector)
    size(a, 1) == length(b) ||
        throw(DimensionMismatch("The length of b must match the number of rows in a."))
    length(r) == size(a, 2) || throw(DimensionMismatch("Incorrect size of r."))
    return size(a)
end

function get_pairwise_dims(r::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix)
    ma, na = size(a)
    mb, nb = size(b)
    ma == mb || throw(DimensionMismatch("The numbers of rows in a and b must match."))
    size(r) == (na, nb) || throw(DimensionMismatch("Incorrect size of r."))
    return (ma, na, nb)
end

function get_pairwise_dims(r::AbstractMatrix, a::AbstractMatrix)
    m, n = size(a)
    size(r) == (n, n) || throw(DimensionMismatch("Incorrect size of r."))
    return (m, n)
end

# for metrics with fixed dimension (e.g. weighted metrics)

function get_colwise_dims(d::Int, r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
    size(a, 1) == size(b, 1) == d ||
        throw(DimensionMismatch("Incorrect vector dimensions."))
    length(r) == size(a, 2) || throw(DimensionMismatch("Incorrect size of r."))
    return size(a)
end

function get_colwise_dims(d::Int, r::AbstractArray, a::AbstractVector, b::AbstractMatrix)
    length(a) == size(b, 1) == d ||
        throw(DimensionMismatch("Incorrect vector dimensions."))
    length(r) == size(b, 2) || throw(DimensionMismatch("Incorrect size of r."))
    return size(b)
end

function get_colwise_dims(d::Int, r::AbstractArray, a::AbstractMatrix, b::AbstractVector)
    size(a, 1) == length(b) == d
        throw(DimensionMismatch("Incorrect vector dimensions."))
    length(r) == size(a, 2) || throw(DimensionMismatch("Incorrect size of r."))
    return size(a)
end

function get_pairwise_dims(d::Int, r::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix)
    na = size(a, 2)
    nb = size(b, 2)
    size(a, 1) == size(b, 1) == d || throw(DimensionMismatch("Incorrect vector dimensions."))
    size(r) == (na, nb) || throw(DimensionMismatch("Incorrect size of r."))
    return (d, na, nb)
end

function get_pairwise_dims(d::Int, r::AbstractMatrix, a::AbstractMatrix)
    n = size(a, 2)
    size(a, 1) == d || throw(DimensionMismatch("Incorrect vector dimensions."))
    size(r) == (n, n) || throw(DimensionMismatch("Incorrect size of r."))
    return (d, n)
end

###########################################################
#
#   calculation
#
###########################################################

function sqrt!(a::AbstractArray)
    @simd for i in eachindex(a)
        @inbounds a[i] = sqrt(a[i])
    end
    a
end

function sumsq_percol{T}(a::AbstractMatrix{T})
    m = size(a, 1)
    n = size(a, 2)
    r = Array(T, n)
    for j = 1:n
        aj = slice(a, :, j)
        r[j] = dot(aj, aj)
    end
    return r
end

function wsumsq_percol{T1, T2}(w::AbstractArray{T1}, a::AbstractMatrix{T2})
    m = size(a, 1)
    n = size(a, 2)
    T = typeof(one(T1)*one(T2))
    r = Array(T, n)
    for j = 1:n
        aj = slice(a, :, j)
        s = zero(T)
        @simd for i = 1:m
            @inbounds s += w[i] * abs2(aj[i])
        end
        r[j] = s
    end
    return r
end

function dot_percol!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
    m = size(a,1)
    n = size(a,2)
    size(b) == (m,n) && length(r) == n ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    for j = 1:n
        aj = slice(a,:,j)
        bj = slice(b,:,j)
        r[j] = dot(aj, bj)
    end
    return r
end

dot_percol(a::AbstractMatrix, b::AbstractMatrix) = dot_percol!(Array(Float64, size(a,2)), a, b)


text:
using LOLTools.MatchV5

api_key = get(ENV, "RIOT_TOKEN", "")
region = "asia"

# platformId = "KR"
# summonerName = "Chooooovyyyyyy"
puuid = "2SAuZFcOENdoO5UWE_h9DmGSU5MFkPgs0eQY-YqZKOSLLN7iOkYYCRm7Z1bx2lqq15KziYpKJ3w7-w"
matches = MatchV5.matches_by_puuid(api_key, region, puuid)
@info :matches matches

matchId = "KR_5424985325"

match = MatchV5.match_by_match_id(api_key, region, matchId)
# @info :match match

match_timeline = MatchV5.match_timeline_by_match_id(api_key, region, matchId)
# @info :match_timeline match_timeline


text:
using Printf # Just the same as import a module
using Statistics
using CPUTime

x = collect(1:0.1:100)

y = cos.(x)

println(x)

println(y)

open("teste.txt","w") do io

    println(io,x)

    for i in y
        println(io,i)
    end

end

max_stars_repo_name: arturofburgos/Comparing-Languages


text:
##
# Splits the userinfo portion of an URI in the format user:password and
# returns the components as tuple.
#
# Note: This is just a convenience method, and this form of usage is
# deprecated as of rfc3986.
# See: http://tools.ietf.org/html/rfc3986#section-3.2.1
function userinfo(uri::URI)
    Base.warn_once("Use of the format user:password is deprecated (rfc3986)")
    uinfo = uri.userinfo
    sep = findfirst(equalto(':'), uinfo)
    l = length(uinfo)
    username = uinfo[1:(sep-1)]
    password = ((sep == l) || (sep == 0)) ? "" : uinfo[(sep+1):l]
    (username, password)
end

##
# Splits the path into components and parameters
# See: http://tools.ietf.org/html/rfc3986#section-3.3
function path_params(uri::URI, seps=[';',',','='])
    elems = split(uri.path, '/', keep = false)
    p = Array[]
    for elem in elems
        pp = split(elem, seps)
        push!(p, pp)
    end
    p
end

##
# Splits the query into key value pairs
function query_params(query::AbstractString)
    elems = split(query, '&', keep = false)
    d = Dict{AbstractString, AbstractString}()
    for elem in elems
        pp = split(elem, "=", keep = true)
        if length(pp) == 2
            d[unescape(pp[1])] = unescape(pp[2])
        else
            d[elem] = "" #gracefully degrade for ill formed query params
        end
    end
    d
end
query_params(uri::URI) = query_params(uri.query::AbstractString)

##
# Create equivalent URI without the fragment
defrag(uri::URI) = URI(uri.scheme, uri.host, uri.port, uri.path, uri.query, "", uri.userinfo, uri.specifies_authority)

##
# Validate known URI formats
const uses_authority = ["hdfs", "ftp", "http", "gopher", "nntp", "telnet", "imap", "wais", "file", "mms", "https", "shttp", "snews", "prospero", "rtsp", "rtspu", "rsync", "svn", "svn+ssh", "sftp" ,"nfs", "git", "git+ssh", "ldap"]
const uses_params = ["ftp", "hdl", "prospero", "http", "imap", "https", "shttp", "rtsp", "rtspu", "sip", "sips", "mms", "sftp", "tel"]
const non_hierarchical = ["gopher", "hdl", "mailto", "news", "telnet", "wais", "imap", "snews", "sip", "sips"]
const uses_query = ["http", "wais", "imap", "https", "shttp", "mms", "gopher", "rtsp", "rtspu", "sip", "sips", "ldap"]
const uses_fragment = ["hdfs", "ftp", "hdl", "http", "gopher", "news", "nntp", "wais", "https", "shttp", "snews", "file", "prospero"]

function isvalid(uri::URI)
    scheme = uri.scheme
    isempty(scheme) && error("Can not validate relative URI")
    slash = findfirst(equalto('/'), uri.path)
    hasslash = slash !== nothing && slash > 1 # For 0.6 compatibility
    if ((scheme in non_hierarchical) && hasslash) ||                          # path hierarchy not allowed
       (!(scheme in uses_query) && !isempty(uri.query)) ||                    # query component not allowed
       (!(scheme in uses_fragment) && !isempty(uri.fragment)) ||              # fragment identifier component not allowed
       (!(scheme in uses_authority) && (!isempty(uri.host) || (0 != uri.port) || !isempty(uri.userinfo))) # authority component not allowed
        return false
    end
    true
end


text:
function AverageTemperature(r)
    sum = 0
    counter = 0
    side = Ost
    while !isborder(r, Nord)
        sum, counter = movementsandcheckformarker(r, side, sum, counter)
        sum, counter = checktemperature(r, sum, counter)
        move!(r, Nord)
        side = inverse(side)
    end
    sum, counter = movementsandcheckformarker(r, side, sum, counter)
    sum, counter = checktemperature(r, sum, counter)
    return sum/counter
end

function inverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function movementsandcheckformarker(r, side, s, c)
    while !isborder(r, side)
        if ismarker(r)
            s += temperature(r)
            c += 1
        end
        move!(r, side)
    end
    return s, c
end

function checktemperature(r, s, c)
    if ismarker(r)
        s += temperature(r)
        c += 1
    end
    return s, c
end


text:
int_rules_1_1_1_5 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.1.5*P(x)*(a+b*x)^m*(c+d*x)^n =#
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n'), ~x) => Antiderivative(~Px * (~a * ~c + ~b * ~d * (~x) ^ 2) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (PolyQ(~Px, ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ(~m, ~n) && (IntegerQ(~m) || GtQ(~a, 0) && GtQ(~c, 0)))))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n, ~x) => (((~a + ~b * ~x) ^ FracPart(~m) * (~c + ~d * ~x) ^ FracPart(~m)) / (~a * ~c + ~b * ~d * (~x) ^ 2) ^ FracPart(~m)) * Antiderivative(~Px * (~a * ~c + ~b * ~d * (~x) ^ 2) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (PolyQ(~Px, ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ(~m, ~n) && Not(IntegerQ(~m)))))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n'), ~x) => Antiderivative(PolynomialQuotient(~Px, ~a + ~b * ~x, ~x) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (PolyQ(~Px, ~x) && EqQ(PolynomialRemainder(~Px, ~a + ~b * ~x, ~x), 0))
    @apply_utils Antiderivative((~Px * (~(c') + ~(d') * ~x) ^ ~(n')) / (~(a') + ~(b') * ~x), ~x) => Antiderivative(ExpandIntegrand(1 / sqrt(~c + ~d * ~x), (~Px * (~c + ~d * ~x) ^ (~n + 1 / 2)) / (~a + ~b * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (PolyQ(~Px, ~x) && (ILtQ(~n + 1 / 2, 0) && GtQ(Expon(~Px, ~x), 2)))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand(~Px * (~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (PolyQ(~Px, ~x) && ((IntegersQ(~m, ~n) || IGtQ(~m, -2)) && GtQ(Expon(~Px, ~x), 2)))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~(n'), ~x) => With([Qx = PolynomialQuotient(~Px, ~a + ~b * ~x, ~x), R = PolynomialRemainder(~Px, ~a + ~b * ~x, ~x)], (R * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1)) / ((~m + 1) * (~b * ~c - ~a * ~d)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * ExpandToSum((~m + 1) * (~b * ~c - ~a * ~d) * Qx - ~d * R * (~m + ~n + 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (PolyQ(~Px, ~x) && (ILtQ(~m, -1) && GtQ(Expon(~Px, ~x), 2)))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~(n'), ~x) => With([Qx = PolynomialQuotient(~Px, ~a + ~b * ~x, ~x), R = PolynomialRemainder(~Px, ~a + ~b * ~x, ~x)], (R * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1)) / ((~m + 1) * (~b * ~c - ~a * ~d)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * ExpandToSum((~m + 1) * (~b * ~c - ~a * ~d) * Qx - ~d * R * (~m + ~n + 2), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && (PolyQ(~Px, ~x) && (LtQ(~m, -1) && GtQ(Expon(~Px, ~x), 2)))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n'), ~x) => With([q = Expon(~Px, ~x), k = Coeff(~Px, ~x, Expon(~Px, ~x))], (k * (~a + ~b * ~x) ^ (~m + q) * (~c + ~d * ~x) ^ (~n + 1)) / (~d * (~b) ^ q * (~m + ~n + q + 1)) + (1 / (~d * (~b) ^ q * (~m + ~n + q + 1))) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * ExpandToSum((~d * (~b) ^ q * (~m + ~n + q + 1) * ~Px - ~d * k * (~m + ~n + q + 1) * (~a + ~b * ~x) ^ q) - k * (~b * ~c - ~a * ~d) * (~m + q) * (~a + ~b * ~x) ^ (q - 1), ~x), ~x) <-- NeQ(~m + ~n + q + 1, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (PolyQ(~Px, ~x) && GtQ(Expon(~Px, ~x), 2))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.5 P(x) (a+b x)^m (c+d x)^n.jl


text:
using MPSGE

m = Model()

@parameter(m, endow, 1.0)

@sector(m, X)
@sector(m, Y)
@sector(m, U)

@commodity(m, PX)
@commodity(m, PY)
@commodity(m, PU)
@commodity(m, PL)
@commodity(m, PK)

@consumer(m, RA, benchmark=150.)

@production(m, X, 1, [Output(PX, 100)], [Input(PL, 50), Input(PK, 50)])
@production(m, Y, 1, [Output(PY, 50)], [Input(PL, 20), Input(PK, 30)])
@production(m, U, 1, [Output(PU, 150)], [Input(PX, 100), Input(PY, 50)])

@demand(m, RA, [Demand(PU, 150)], [Endowment(PL, :(70 * $endow)), Endowment(PK, 80.)])

solve!(m, cumulative_iteration_limit=0)

set_fixed!(PX, true)
set_value(endow, 1.1)

solve!(m)

algebraic_version(m)


text:
export
    infoString,
    getDeviceCount,
    setDevice,
    getDevice,
    deviceInfo,
    deviceInfos,
    setSeed,
    sync

immutable AFDevice <: AFImpl
    infoString::Ptr{Void}
    getDeviceCount::Ptr{Void}
    setDevice::Ptr{Void}
    getDevice::Ptr{Void}
    deviceInfo::Ptr{Void}
    setSeed::Ptr{Void}
    sync::Ptr{Void}

    function AFDevice(ptr)
        new(
            try
                Libdl.dlsym(ptr, :af_info_string)
            catch
                C_NULL
            end,
            Libdl.dlsym(ptr, :af_get_device_count),
            Libdl.dlsym(ptr, :af_set_device),
            Libdl.dlsym(ptr, :af_get_device),
            Libdl.dlsym(ptr, :af_device_info),
            Libdl.dlsym(ptr, :af_set_seed),
            Libdl.dlsym(ptr, :af_sync)
        )
    end
end

function infoString(af::ArrayFire)
    if (af.device.infoString != C_NULL)
        result = Ref{Ptr{Cchar}}()
        ccall(af.device.infoString, Cint, (Ptr{Ptr{Cchar}}, Cuchar), result, false)
        bytestring(result[])
    else
        info = deviceInfo(af)
        convert(AbstractString, info)
    end
end

function getDeviceCount(af::ArrayFire)
    count = Ref{Cint}(0)
    err = ccall(af.device.getDeviceCount, Cint, (Ref{Cint},), count)
    assertErr(err)
    count[]
end

function setDevice(af::ArrayFire, device)
    err = ccall(af.device.setDevice, Cint, (Cint,), device)
    assertErr(err)
end

function getDevice(af::ArrayFire)
    device = Ref{Cint}(0)
    err = ccall(af.device.getDevice, Cint, (Ref{Cint},), device)
    assertErr(err)
    device[]
end

function deviceInfo(af::ArrayFire)
    name = Array(Cchar, 256)
    platform = Array(Cchar, 256)
    toolkit = Array(Cchar, 256)
    compute = Array(Cchar, 256)
    err = ccall(
        af.device.deviceInfo,
        Cint,
        (Ptr{Cchar}, Ptr{Cchar}, Ptr{Cchar}, Ptr{Cchar}),
        name, platform, toolkit, compute)
    assertErr(err)
    DeviceInfo(
        getDevice(af),
        bytestring(pointer(name)),
        bytestring(pointer(platform)),
        bytestring(pointer(toolkit)),
        bytestring(pointer(compute)))
end

function deviceInfos(af::ArrayFire)
    curr = getDevice(af)
    count = getDeviceCount(af)
    result = Array(DeviceInfo, count)
    try
        for id = 0:(count - 1)
            setDevice(af, id)
            result[id + 1] = deviceInfo(af)
        end
    finally
        setDevice(af, curr)
    end
    result
end

function setSeed(af::ArrayFire, seed)
    err = ccall(af.device.setSeed, Cint, (UInt64,), UInt64(seed))
    assertErr(err)
end

function sync(af::ArrayFire)
    err = ccall(af.device.sync, Cint, (Int32,), getDevice(af))
    assertErr(err)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# test using data
using Polylogarithms
using SpecialFunctions
using Test
using DataFrames, CSV
import Base.MathConstants: π, pi, ℯ, e, γ, eulergamma, catalan, φ, golden
include("test_defs.jl")
using Printf

# desired accuracy is 1.0e-12, but we get a few points above this so
≈(a,b) = near_equal(a,b,1.0e-10)

L = Symbol("Li_s(z)")
@testset "Polylogarithm polylog on data" begin
    for C=1:3
        filename = @sprintf("../data/polylog_test_data_rand_%d.csv", C)
        data1 = CSV.read(filename, DataFrame; delim=",", type=String)

        # has trouble reading in numbers like "2." so read all into strings, and parse
        data1[!,:s] = parse.(Complex{Float64}, data1[!,:s] )
        data1[!,:z] = parse.(Complex{Float64}, data1[!,:z] )
        data1[!, L] = parse.(Complex{Float64}, data1[!, L] )

        m = size(data1,1)
        Li = data1[!, L]
        s  = data1[!,:s]
        z  = data1[!,:z]

        for i=1:m
            @test polylog( s[i], z[i] ) ≈ Li[i]
        end

    end
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 15b6d7c8-f787-11ea-0fd9-4738ba47533d
using Pkg, DrWatson

# ╔═╡ 15b712ce-f787-11ea-245f-55620c5d1b05
begin
	@quickactivate "StatisticalRethinkingTuring"
	using StatisticalRethinking
end

# ╔═╡ bd4d25cc-f786-11ea-18fd-3326d3896012
md"## Fig3.3t.jl"

# ╔═╡ 15b79bfe-f787-11ea-010f-5d6a113c0f19
begin
	p_grid = range(0, step=0.000001, stop=1)
	prior = ones(length(p_grid))
	likelihood = pdf.(Binomial.(3, p_grid), 3)
	posterior = likelihood .* prior
	posterior = posterior / sum(posterior)
	samples = sample(p_grid, Weights(posterior), length(p_grid));
end;

# ╔═╡ 15bf2d08-f787-11ea-1316-e71570fded77
begin
	b1 = quantile(samples, [0.25, 0.75])
	b2 = hpdi(samples, alpha=0.5)

	p1 = plot_density_interval(samples, b1;
	  xlab="Proportion water (p)", title="50% PI");
	p2 = plot_density_interval(samples, b2;
	  xlab="Proportion water (p)", title="50% HPDI");
end;

# ╔═╡ 15bf932c-f787-11ea-295b-6f68b4c100ec
plot(p1, p2, layout=(1, 2))

# ╔═╡ 15cd1a88-f787-11ea-33fa-4f5d7cebab01
md"## End of Fig3.3t.jl"

# ╔═╡ Cell order:
# ╟─bd4d25cc-f786-11ea-18fd-3326d3896012
# ╠═15b6d7c8-f787-11ea-0fd9-4738ba47533d
# ╠═15b712ce-f787-11ea-245f-55620c5d1b05
# ╠═15b79bfe-f787-11ea-010f-5d6a113c0f19
# ╠═15bf2d08-f787-11ea-1316-e71570fded77
# ╠═15bf932c-f787-11ea-295b-6f68b4c100ec
# ╟─15cd1a88-f787-11ea-33fa-4f5d7cebab01


text:
include("findLargestInCol.jl")
include("swapRowsV2.jl")

# we write another function that composes A into A = P^{-1}LU where P is stored as a vector and swapRows is more efficient
function PALUdecompV2(A)

    # we get the size of A
    m, n = size(A)

    # we check that A is indeed a square matrix
    if m != n
        error("Input must be a square matrix")
    end

    # we initialize L and U to zero, P to {1,2, ..., n}
    L = zeros(n, n)
    U = zeros(n, n)
    p = collect(1:n)

    # we make a copy of A to work with
    B = copy(A)

    for k = 1:n

        # we find the pivot
        maxk = findLargestInCol(B, k)
        if maxk != k
            swapRowsV2!(B, maxk, k)
            swapRowsV2!(L, maxk, k)

            # we swap two elements of P instead of two rows
            temp = p[maxk]
            p[maxk] = p[k]
            p[k] = temp
        end

        # we proceed with the normal LU decomposition algorithm
        for j = k:n
            U[k,j] = B[k,j]
        end

        # we check that we don't divide by zero
        if U[k,k] == 0
            error("** A^(k-1)_{k,k} == 0 in PALU decomp")
        end

        for i = k:n
           L[i,k] = B[i,k] / U[k,k]
        end

        for i = k:n
            for j = k:n
                B[i,j] -= L[i,k] * U[k,j];
            end
        end
    end

    P = zeros(n,n)
    for i = 1:n
        P[i,p[i]] = 1
    end

    # we return P, L and U
    return P, L, U
end

max_stars_repo_path: 1a/LU with pivot/PALUdecompV2.jl

max_issues_repo_path: 1a/LU with pivot/PALUdecompV2.jl

max_forks_repo_path: 1a/LU with pivot/PALUdecompV2.jl


text:
# https://github.com/grasingerm/lbxflow/blob/master/inc/lattice.jl
# http://mafija.fmf.uni-lj.si/seminar/files/2012_2013/Lattice_Boltzmann_method.pdf

# const original_order = [1, 2 ,3, 4]

struct D2Q4{
    T <: Real,
    Abscissaes <: AbstractMatrix{<:Integer},
    Weights <: AbstractVector{<:Real},
} <: Quadrature
    abscissae::Abscissaes
    weights::Weights
    speed_of_sound_squared::T
end

D2Q4() = D2Q4(
    [
        1 0 -1 0
        0 1 0 -1
    ],
    [1 / 4, 1 / 4, 1 / 4, 1 / 4],
    2.0,
)

order(q::D2Q4) = 3
function opposite(q::D2Q4, idx::Int64)
    if idx <= 2
        return idx + 2
    end
    return idx - 2
end

Base.show(io::IO, q::D2Q4) = show(io, "D2Q4")
Base.string(q::D2Q4) = "D2Q4"


text:
using Diffractor
using Diffractor: var"'", ∂⃖
using ChainRules
using ChainRulesCore
using ChainRules: ZeroTangent, NoTangent
using Symbolics

using Test

# Unit tests
function tup2(f)
    a, b = ∂⃖{2}()(f, 1)
    c, d = b((2,))
    e, f = d(ZeroTangent(), 3)
    f((4,))
end

@test tup2(tuple) == (NoTangent(), 4)

my_tuple(args...) = args
ChainRules.rrule(::typeof(my_tuple), args...) = args, Δ->Core.tuple(NoTangent(), Δ...)
@test tup2(my_tuple) == (ZeroTangent(), 4)

# Check characteristic of exp rule
@variables ω α β γ δ ϵ ζ η
(x1, c1) = ∂⃖{3}()(exp, ω)
@test simplify(x1 == exp(ω)).val
((_, x2), c2) = c1(α)
@test simplify(x2 == α*exp(ω)).val
(x3, c3) = c2(ZeroTangent(), β)
@test simplify(x3 == β*exp(ω)).val
((_, x4), c4) = c3(γ)
@test simplify(x4 == exp(ω)*(γ + (α*β))).val
(x5, c5) = c4(ZeroTangent(), δ)
@test simplify(x5 == δ*exp(ω)).val
((_, x6), c6) = c5(ϵ)
@test simplify(x6 == ϵ*exp(ω) + α*δ*exp(ω)).val
(x7, c7) = c6(ZeroTangent(), ζ)
@test simplify(x7 == ζ*exp(ω) + β*δ*exp(ω)).val
(_, x8) = c7(η)
@test simplify(x8 == (η + (α*ζ) + (β*ϵ) + (δ*(γ + (α*β))))*exp(ω)).val

# Minimal 2-nd order forward smoke test
@test Diffractor.∂☆{2}()(Diffractor.ZeroBundle{2}(sin),
    Diffractor.TangentBundle{2}(1.0, (1.0, 1.0, 0.0)))[Diffractor.CanonicalTangentIndex(1)] == sin'(1.0)

function simple_control_flow(b, x)
    if b
        return sin(x)
    else
        return cos(x)
    end
end

function myprod(xs)
    s = 1
    for x in xs
      s *= x
    end
    return s
end

function mypow(x, n)
    r = one(x)
    while n > 0
        n -= 1
        r *= x
    end
    return r
end

function times_three_while(x)
    z = x
    i = 3
    while i > 1
        z += x
        i -= 1
    end
    z
end

isa_control_flow(::Type{T}, x) where {T} = isa(x, T) ? x : T(x)

# Simple Reverse Mode tests
let var"'" = Diffractor.PrimeDerivativeBack
    # Integration tests
    @test @inferred(sin'(1.0)) == cos(1.0)
    @test @inferred(sin''(1.0)) == -sin(1.0)
    @test sin'''(1.0) == -cos(1.0)
    @test sin''''(1.0) == sin(1.0)
    @test sin'''''(1.0) == cos(1.0)
    @test sin''''''(1.0) == -sin(1.0)

    f_getfield(x) = getfield((x,), 1)
    @test f_getfield'(1) == 1
    @test f_getfield''(1) == 0
    @test f_getfield'''(1) == 0

    # Higher order mixed mode tests

    complicated_2sin(x) = (x = map(sin, Diffractor.xfill(x, 2)); x[1] + x[2])
    @test @inferred(complicated_2sin'(1.0)) == 2sin'(1.0)
    @test @inferred(complicated_2sin''(1.0)) == 2sin''(1.0)
    @test @inferred(complicated_2sin'''(1.0)) == 2sin'''(1.0)
    @test @inferred(complicated_2sin''''(1.0)) == 2sin''''(1.0)

    # Control flow cases
    @test @inferred((x->simple_control_flow(true, x))'(1.0)) == sin'(1.0)
    @test @inferred((x->simple_control_flow(false, x))'(1.0)) == cos'(1.0)
    @test (x->sum(isa_control_flow(Matrix{Float64}, x)))'(Float32[1 2;]) == [1.0 1.0;]
    @test times_three_while'(1.0) == 3.0

    pow5p(x) = (x->mypow(x, 5))'(x)
    @test pow5p(1.0) == 5.0
end

# Simple Forward Mode tests
let var"'" = Diffractor.PrimeDerivativeFwd
    recursive_sin(x) = sin(x)
    ChainRulesCore.frule(∂, ::typeof(recursive_sin), x) = frule(∂, sin, x)

    # Integration tests
    @test recursive_sin'(1.0) == cos(1.0)
    @test recursive_sin''(1.0) == -sin(1.0)
    @test recursive_sin'''(1.0) == -cos(1.0)
    @test recursive_sin''''(1.0) == sin(1.0)
    @test recursive_sin'''''(1.0) == cos(1.0)
    @test recursive_sin''''''(1.0) == -sin(1.0)

    # Test the special rules for sin/cos/exp
    @test sin''''''(1.0) == -sin(1.0)
    @test cos''''''(1.0) == -cos(1.0)
    @test exp''''''(1.0) == exp(1.0)
end

# Some Basic Mixed Mode tests
function sin_twice_fwd(x)
    let var"'" = Diffractor.PrimeDerivativeFwd
            sin''(x)
    end
end
let var"'" = Diffractor.PrimeDerivativeFwd
    @test sin_twice_fwd'(1.0) == sin'''(1.0)
end

# Regression tests
@test gradient(x -> sum(abs2, x .+ 1.0), zeros(3))[1] == [2.0, 2.0, 2.0]

const fwd = Diffractor.PrimeDerivativeFwd
const bwd = Diffractor.PrimeDerivativeFwd

function f_broadcast(a)
    l = a / 2.0 * [[0. 1. 1.]; [1. 0. 1.]; [1. 1. 0.]]
    return sum(l)
end
@test fwd(f_broadcast)(1.0) == bwd(f_broadcast)(1.0)

# Make sure that there's no infinite recursion in kwarg calls
g_kw(;x=1.0) = sin(x)
f_kw(x) = g_kw(;x)
@test bwd(f_kw)(1.0) == bwd(sin)(1.0)

function f_crit_edge(a, b, c, x)
    # A function with two critical edges. This used to trigger an issue where
    # Diffractor would fail to insert edges for the second split critical edge.
    y = 1x
    if a && b
        y = 2x
    end
    if b && c
        y = 3x
    end

    if c
        y = 4y
    end

    return y
end
@test bwd(x->f_crit_edge(false, false, false, x))(1.0) == 1.0
@test bwd(x->f_crit_edge(true, true, false, x))(1.0) == 2.0
@test bwd(x->f_crit_edge(false, true, true, x))(1.0) == 12.0
@test bwd(x->f_crit_edge(false, false, true, x))(1.0) == 4.0

# Issue #27 - Mixup in lifting of getfield
let var"'" = bwd
    @test (x->x^5)''(1.0) == 20.
    @test (x->x^5)'''(1.0) == 60.
end

# Issue #38 - Splatting arrays
@test gradient(x -> max(x...), (1,2,3))[1] == (0.0, 0.0, 1.0)
@test gradient(x -> max(x...), [1,2,3])[1] == [0.0, 0.0, 1.0]

# Issue #40 - Symbol type parameters not properly quoted
@test Diffractor.∂⃖recurse{1}()(Val{:transformations})[1] === Val{:transformations}()

include("pinn.jl")