text:
using CCBlade
using FLOWMath
using SpecialFunctions: gamma
import LinearAlgebra: dot
using Snopt
import ForwardDiff
import ReverseDiff
import SparseArrays, SparseDiffTools
import FiniteDiff
import NLsolve

function residual2d(as, aps, rotor, section, op)

    phi = atan((1 + as)*op.Vx/((1 - aps)*op.Vy))
    _, outputs = CCBlade.residual(phi, rotor, section, op)

    R = [-outputs.a - as; 10*(-outputs.ap - aps)]

    return R, outputs

end

function solve2d(rotor, section, op)

    # ----- solve residual function ------

    # wrapper to residual function
    R(x) = residual2d(x[1], x[2], rotor, section, op)[1]

    a0 = -1.0/3
    ap0 = -0.1
    x0 = [a0; ap0]

    results = NLsolve.nlsolve(R, x0, method=:trust_region, ftol=1e-14, autodiff=:forward)
    if !NLsolve.converged(results)
        println("not converged")
        println(rotor)
        println(section)
        println(op)
    end
    x = results.zero

    _, outputs = residual2d(x[1], x[2], rotor, section, op)

    return outputs
end

function residualadjoint(R, ymin, ymax, x)
    y, _ = FLOWMath.brent(y -> R(x, y), ymin, ymax)
    return y
end

function residualadjoint(R, ymin, ymax, d::Vector{D}) where D <: ForwardDiff.Dual

    # solve root
    x = ForwardDiff.value.(d)
    y = residualadjoint(R, ForwardDiff.value(ymin), ForwardDiff.value(ymax), x)

    # compute derivatives
    wrap(xy) = R(xy[1:end-1], xy[end])
    g = Vector{Float64}(undef, length(x) + 1)
    ForwardDiff.gradient!(g, wrap, [x; y])
    drdx = g[1:end-1]
    drdy = g[end]
    dydx = -drdx/drdy

    # assemble dual number
    b_in = zip(collect.(ForwardDiff.partials.(d))...)
    b_arr = map(x->dot(dydx, x), b_in)
    p = ForwardDiff.Partials((b_arr...,))
    return D(y, p)
end

function solveadjoint(rotor, section, op)

    # error handling
    if typeof(section) <: Vector
        error("You passed in an vector for section, but this funciton does not accept an vector.\nProbably you intended to use broadcasting (notice the dot): solve.(Ref(rotor), sections, ops)")
    end

    # check if we are at hub/tip
    if isapprox(section.r, rotor.Rhub, atol=1e-6) || isapprox(section.r, rotor.Rtip, atol=1e-6)
        return Outputs()  # no loads at hub/tip
    end

    # parameters
    npts = 10  # number of discretization points to find bracket in residual solve

    # unpack
    Vx = op.Vx
    Vy = op.Vy
    theta = section.theta + op.pitch

    # ---- determine quadrants based on case -----
    Vx_is_zero = isapprox(Vx, 0.0, atol=1e-6)
    Vy_is_zero = isapprox(Vy, 0.0, atol=1e-6)

    # quadrants
    epsilon = 1e-6
    q1 = [epsilon, pi/2]
    q2 = [-pi/2, -epsilon]
    q3 = [pi/2, pi-epsilon]
    q4 = [-pi+epsilon, -pi/2]

    if Vx_is_zero && Vy_is_zero
        return Outputs()

    elseif Vx_is_zero

        startfrom90 = false  # start bracket at 0 deg.

        if Vy > 0 && theta > 0
            order = (q1, q2)
        elseif Vy > 0 && theta < 0
            order = (q2, q1)
        elseif Vy < 0 && theta > 0
            order = (q3, q4)
        else  # Vy < 0 && theta < 0
            order = (q4, q3)
        end

    elseif Vy_is_zero

        startfrom90 = true  # start bracket search from 90 deg

        if Vx > 0 && abs(theta) < pi/2
            order = (q1, q3)
        elseif Vx < 0 && abs(theta) < pi/2
            order = (q2, q4)
        elseif Vx > 0 && abs(theta) > pi/2
            order = (q3, q1)
        else  # Vx < 0 && abs(theta) > pi/2
            order = (q4, q2)
        end

    else  # normal case

        startfrom90 = false

        if Vx > 0 && Vy > 0
            order = (q1, q2, q3, q4)
        elseif Vx < 0 && Vy > 0
            order = (q2, q1, q4, q3)
        elseif Vx > 0 && Vy < 0
            order = (q3, q4, q1, q2)
        else  # Vx[i] < 0 && Vy[i] < 0
            order = (q4, q3, q2, q1)
        end

    end

    # ----- solve residual function ------

    # # wrapper to residual function to accomodate format required by fzero
    R(phi) = CCBlade.residual(phi, rotor, section, op)[1]

    function R2(x, y)
        section2 = Section(section.r, x[2], x[3], section.af)
        op2 = OperatingPoint(x[4], x[5], op.rho, x[1], op.mu, op.asound)

        return CCBlade.residual(y, rotor, section2, op2)[1]
    end
    x2 = [op.pitch, section.chord, section.theta, op.Vx, op.Vy]

    success = false
    for j = 1:length(order)  # quadrant orders.  In most cases it should find root in first quadrant searched.
        phimin, phimax = order[j]

        # check to see if it would be faster to reverse the bracket search direction
        backwardsearch = false
        if !startfrom90
            if phimin == -pi/2 || phimax == -pi/2  # q2 or q4
                backwardsearch = true
            end
        else
            if phimax == pi/2  # q1
                backwardsearch = true
            end
        end

        # force to dual numbers if necessary
        phimin = phimin*one(section.chord)
        phimax = phimax*one(section.chord)

        # find bracket
        success, phiL, phiU = CCBlade.firstbracket(R, phimin, phimax, npts, backwardsearch)

        # once bracket is found, solve root finding problem and compute loads
        if success
            # phistar, _ = FLOWMath.brent(R, phiL, phiU)
            phistar = residualadjoint(R2, phiL, phiU, x2)
            _, outputs = CCBlade.residual(phistar, rotor, section, op)
            return outputs
        end
    end

    # it shouldn't get to this point.  if it does it means no solution was found
    # it will return empty outputs
    # alternatively, one could increase npts and try again

    @warn "Invalid data (likely) for this section.  Zero loading assumed."
    return Outputs()
end

function shape(cspline, tspline)

    Rhub = 1.5
    Rtip = 63.0
    r = [2.8667, 5.6000, 8.3333, 11.7500, 15.8500, 19.9500, 24.0500,
        28.1500, 32.2500, 36.3500, 40.4500, 44.5500, 48.6500, 52.7500,
        56.1667, 58.9000, 61.6333]

    rspline = [0.00, 0.25, 0.5, 0.75, 1.0]*(Rtip-Rhub) .+ Rhub
    chord = akima(rspline, cspline, r)

    rspline2 = [0.11111, 0.4, 0.7, 1.0]*(Rtip-Rhub) .+ Rhub
    theta = [[tspline[1], tspline[1]]; akima(rspline2, tspline, r[3:end])]  # fixed root twist b/c it is a cylinder

    return r, chord, theta, Rhub, Rtip
end

function turbineopt(cspline, tspline, tsr, pitch, method)

    # define pitch to start at zero, add root pitch to twist instead.
    p0 = pitch[1]
    pitch .-= p0
    tspline .+= p0

    # ------- NREL 5MW --------

    B = 3
    turbine = true
    precone = 2.5*pi/180

    r, chord, theta, Rhub, Rtip = shape(cspline, tspline)
    rotor = Rotor(Rhub, Rtip, B; turbine, precone)

    # chord = [3.542, 3.854, 4.167, 4.557, 4.652, 4.458, 4.249, 4.007, 3.748,
    #     3.502, 3.256, 3.010, 2.764, 2.518, 2.313, 2.086, 1.419]
    # theta = pi/180*[13.308, 13.308, 13.308, 13.308, 11.480, 10.162, 9.011, 7.795,
    #     6.544, 5.361, 4.188, 3.125, 2.319, 1.526, 0.863, 0.370, 0.106]

    # Define airfoils.  In this case we have 8 different airfoils that we load into an array.
    # These airfoils are defined in files.
    aftypes = Array{AlphaAF}(undef, 8)
    aftypes[1] = AlphaAF("airfoils/Cylinder1.dat", radians=false)
    aftypes[2] = AlphaAF("airfoils/Cylinder2.dat", radians=false)
    aftypes[3] = AlphaAF("airfoils/DU40_A17.dat", radians=false)
    aftypes[4] = AlphaAF("airfoils/DU35_A17.dat", radians=false)
    aftypes[5] = AlphaAF("airfoils/DU30_A17.dat", radians=false)
    aftypes[6] = AlphaAF("airfoils/DU25_A17.dat", radians=false)
    aftypes[7] = AlphaAF("airfoils/DU21_A17.dat", radians=false)
    aftypes[8] = AlphaAF("airfoils/NACA64_A17.dat", radians=false)

    # indices correspond to which airfoil is used at which station
    af_idx = [1, 1, 2, 3, 4, 4, 5, 6, 6, 7, 7, 8, 8, 8, 8, 8, 8]

    # create airfoil array
    airfoils = aftypes[af_idx]

    sections = Section.(r, chord, theta, airfoils)

    # operating point for the turbine
    yaw = 0.0*pi/180
    tilt = 5.0*pi/180
    hubHt = 90.0
    shearExp = 0.2

    rotorR = Rtip*cos(precone)
    azimuth = 0.0*pi/180
    rho = 1.225  #*one(pitch[1])

    nV = length(pitch)
    Vin = 3.0
    Vout = 25.0
    V = range(Vin, Vout, length=nV)
    # Omega_min = 0.0
    Omega_max = 12.0*pi/30.0
    Omega = min.(collect(V).*tsr/rotorR, Omega_max)  # collect needed for ForwardDiff
    P = similar(V, eltype(chord[1]))
    Q = similar(V, eltype(chord[1]))
    T = similar(V, eltype(chord[1]))

    if method == "traditional"
        solver = solve2d
    elseif method == "adjoint"
        solver = solveadjoint
    else
        solver = solve
    end

    for i = 1:nV
        ops = windturbine_op.(V[i], Omega[i], pitch[i], r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho)
        outputs = solver.(Ref(rotor), sections, ops)

        T[i], Q[i] = thrusttorque(rotor, sections, outputs)
        P[i] = Q[i]*Omega[i]
    end

    i = 30  # near rated
    ops = windturbine_op.(V[i], Omega[i], zero(rho), r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho)
    outputs = solver.(Ref(rotor), sections, ops)
    Np = outputs.Np

    k = 2.0  # weibull shape
    Vbar = 6.0  # mean speed of weibull
    A = Vbar / gamma(1.0 + 1.0/k)
    cdf = 1.0 .- exp.(-(V/A).^k)
    AEP = trapz(cdf, P) * 365*24

    return AEP, V, T, P, Q, Np
end
# --------------------------

function run(method, nP)

    # starting point
    cspline = [3.542, 4.6035, 3.748, 2.8255, 1.419]
    tspline = [13.308, 8.0, 4.1, 0.05]*pi/180
    tsr = 7.55
    pitch = 20*pi/180*ones(nP)

    x0 = [cspline; tspline; tsr; pitch]
    # bounds
    lb = [0.5*ones(5); -5*pi/180*ones(4); 1.0; zeros(nP)]
    ub = [10*ones(5); 20*pi/180*ones(4); 15.0; 30.0*pi/180*ones(nP)]
    lb[1] = 3.542  # lower bound on root chord, for pitch bearing

    Prated = 5e6

    function wrappersparse!(y, x)

        cspline2 = x[1:5]
        tspline2 = x[6:9]
        tsr2 = x[10]
        idx = 10
        pitch2 = x[idx+1:idx+nP]

        AEP, V, T, P, Q, Np = turbineopt(cspline2, tspline2, tsr2, pitch2, method)

        # constraints
        # Qmax = Prated/Omega_max
        Tmax = 600e3
        y .= [P/Prated .- 1; T/Tmax .- 1; pitch2[1:end-1] .- pitch2[2:end]; Np/6500 .- 1]  #; Omega[1:end-1] .- Omega[2:end]]

        return -AEP/1e10  # maximize AEP
    end

    function wrapper!(y, x)
        yv = @view y[2:end]
        y[1] = wrappersparse!(yv, x)
    end

    y = zeros(1 + 2*nP + nP-1 + 17)
    config = ForwardDiff.JacobianConfig(wrapper!, y, x0)
    J = zeros(length(y), length(x0))
    dfdx = zeros(length(x0))

    if method == "sparse"
        ysp = zeros(2*nP + nP-1 + 17)
        Jspstart = zeros(length(ysp), length(x0))
        configsp = ForwardDiff.JacobianConfig(wrappersparse!, ysp, x0)
        ForwardDiff.jacobian!(Jspstart, wrappersparse!, ysp, x0, configsp)
        Jsp = SparseArrays.sparse(Jspstart)
        colors = SparseDiffTools.matrix_colors(Jsp)
        cache = SparseDiffTools.ForwardColorJacCache(wrappersparse!, x0, dx=ysp, colorvec=colors, sparsity=Jsp)

        # setup for AEP derivative
        # Prated = 5e6
        Vin = 3.0
        Vout = 25.0
        Vvec = range(Vin, Vout, length=nP)
        k = 2.0  # weibull shape
        Vbar = 6.0  # mean speed of weibull
        A = Vbar / gamma(1.0 + 1.0/k)
        cdf = 1.0 .- exp.(-(Vvec/A).^k)
    end

    function optimize(x)

        wrapper!(y, x)
        f = y[1]
        g = y[2:end]

        if method == "fd" || method == "traditional"
            FiniteDiff.finite_difference_jacobian!(J, wrapper!, x)
            dfdx .= J[1, :]
            dgdx = J[2:end, :]
        elseif method == "dense" || method == "adjoint"
            ForwardDiff.jacobian!(J, wrapper!, y, x, config)
            dfdx .= J[1, :]
            dgdx = J[2:end, :]
        elseif method == "sparse"
            SparseDiffTools.forwarddiff_color_jacobian!(Jsp, wrappersparse!, x, cache)
            dgdx = Jsp
            dPdx = Prated*Jsp[1:nP, :]
            dfdx .= 0.0
            for i = 1:nP-1
                dfdx .+= (cdf[i+1]-cdf[i])*0.5*(dPdx[i, :] + dPdx[i+1, :]) *365*24/-1e10
            end
        end

        fail = false

        return f, g, dfdx, dgdx, fail
    end

    options = Dict{String, Any}()
    options["Derivative option"] = 1
    # options["Verify level"] = 1

    xopt, fopt, info = snopt(optimize, x0, lb, ub, options)
    println(info)

    return xopt, fopt
end

nP = 80

# cspline = [3.542, 4.6035, 3.748, 2.8255, 1.419]
# tspline = [13.308, 8.0, 4.1, 0.05]*pi/180
# tsr = 7.55
# pitch = 20*pi/180*ones(nP)

# times from snopt output file
xopt, fopt = run("sparse", nP)  # 4.97 seconds  |  5.0 seconds
# xopt, fopt = run("adjoint", nP)  12.35 seconds
# xopt, fopt = run("dense", nP)  #  24.33 seconds
# xopt, fopt = run("fd", nP)  # 66.65 seconds
# xopt, fopt = run("traditional", nP)  # 111.86 seconds

# [3.542, 4.760045358177065, 2.9965715561082655, 2.3398138560921455, 1.2066787649473403, 0.23217027516809968, 0.04756986696388409, -0.007796061979556598, -0.06376494241659889, 9.103517237624589, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.04947636361503202, 0.05797161687843067, 0.06729546915857863, 0.07611472657716972, 0.084545916198257, 0.10238931381185744, 0.12193484245765156, 0.13840114476780663, 0.1528780920046593, 0.1661148262225162, 0.17844270542770568, 0.19005323615677938, 0.2010597330827229, 0.21154712594386327, 0.22161200026472808, 0.23132359234067848, 0.24071807104749224, 0.24983317401699293, 0.2587050277588686, 0.26736142750083625, 0.275818390328649, 0.2840903067458611, 0.29219122920245194, 0.30013125450993156, 0.3079195503900656, 0.31556651476196457, 0.3230785743491701, 0.33045441669536774, 0.33769172890657534, 0.34479738484667255, 0.3517824782971913, 0.35866274762843314, 0.3654573834688694, 0.3721824701081526, 0.3788485527218031, 0.3854667071235377, 0.3920411537844105, 0.3985651750312578, 0.4050348112927709, 0.4114413018179495, 0.4177705818285269, 0.424010858860524, 0.43016182781798307, 0.43623015160814443, 0.44222484500848497, 0.4481533871498382, 0.45402229272571143, 0.4598371291637423, 0.4656041453005771, 0.47133119764929, 0.4770195808895922, 0.4826679312122287, 0.48827517063016895, 0.4938451545400422, 0.4993741840234126, 0.5048596100513677, 0.5103039515039579]
# [3.542, 4.760692986206912, 2.9966783734883844, 2.3398300971260957, 1.2066823140213758, 0.22486275809414039, 0.04025369182394284, -0.015121128599521107, -0.07109139144185356, 9.10344339126887, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.056802086563028875, 0.06529939395248592, 0.0746233170108889, 0.08344266230613664, 0.09187398455874778, 0.10971722887313595, 0.12926238507060298, 0.14572835507190782, 0.16020510529330254, 0.17344180882707613, 0.185769643429235, 0.19738011403457678, 0.20838655097548162, 0.21887391496737257, 0.22893877752369987, 0.23865036183185256, 0.248044834033399, 0.2571599064738573, 0.2660317059572908, 0.2746880750712266, 0.2831450233816709, 0.2914169310236453, 0.2995178525344254, 0.3074578820449117, 0.3152461848941429, 0.32289315744208563, 0.33040522505092784, 0.3377810791699985, 0.34501840647341464, 0.3521240805874463, 0.359109193185985, 0.3659894849247997, 0.37278414016763617, 0.3795092409983166, 0.38617533446656865, 0.392793502288733, 0.3993679619147863, 0.4058919815432788, 0.4123616227379301, 0.4187681259334287, 0.42509741709573906, 0.4313377061687048, 0.4374886904949092, 0.4435570338497668, 0.4495517482205652, 0.45548031214038714, 0.4613492397910199, 0.4671641062716087, 0.4729311514263918, 0.4786582305690257, 0.4843466467170885, 0.4899950254275976, 0.4956022905346512, 0.5011722971998086, 0.5067013431822566, 0.5121867834941963, 0.5176311406324784]
# [3.542, 4.760747241419858, 2.9966625418796187, 2.3398271539260125, 1.2066810840710032, 0.22503666496505845, 0.04042321256815731, -0.014951527586733918, -0.07092131728433311, 9.103418446605376, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.056631997258220754, 0.06512948262810273, 0.0744534159072778, 0.0832727709815761, 0.09170410415865786, 0.10954735239460169, 0.1290925193198949, 0.14555851409497375, 0.16003528352034246, 0.17327199763476067, 0.18559984471682742, 0.1972103276774164, 0.20821677632993488, 0.21870415013865807, 0.22876901949953213, 0.238480610864331, 0.24787509086368356, 0.2569901670989968, 0.26586196647295973, 0.2745183372071666, 0.28297529017334205, 0.2912472029169533, 0.2993481297377437, 0.30728816488127186, 0.31507647343257833, 0.3227234508079986, 0.3302355230935136, 0.33761138092488335, 0.344848711345666, 0.35195438915003024, 0.35893950592318025, 0.36581980273914405, 0.37261446369619056, 0.3793395712226177, 0.38600567202645814, 0.3926238471871718, 0.3991983132541477, 0.4057223359746752, 0.4121919796750226, 0.4185984843296204, 0.42492777583049546, 0.43116806480057174, 0.43731904936252003, 0.44338739364408275, 0.44938210986803434, 0.4553106765514783, 0.4611796077758426, 0.4669944790197102, 0.47276152983256087, 0.47848861472247267, 0.4841770372462505, 0.4898254207202264, 0.4954326906840798, 0.5010027021538521, 0.50653175123733, 0.5120171944853323, 0.5174615544318769]
# [3.542, 4.760640547347785, 2.9966473316771842, 2.3398212666796043, 1.2066945723590519, 0.24268093081920974, 0.05807326350597842, 0.002696531844767569, -0.05326937474340002, 9.103431500321744, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.038982613992732104, 0.04747960751292016, 0.056803519249397874, 0.06562284868821903, 0.07405415281081855, 0.0918974830513985, 0.1114426835897507, 0.1279087208338638, 0.14238550926832594, 0.15562222031103515, 0.16795006081478, 0.1795605349425147, 0.19056698155453236, 0.20105435160205584, 0.21111921382505322, 0.22083080113894826, 0.23022527856463476, 0.23934035714942842, 0.24821216001440732, 0.2568685309560612, 0.265325481370162, 0.273597390326586, 0.28169831135061674, 0.2896383398675135, 0.2974266415611344, 0.3050736121831604, 0.3125856774534014, 0.31996152384512805, 0.32719884212739436, 0.3343045104327401, 0.3412896219185127, 0.3481699148002591, 0.35496457450723745, 0.36168968874323887, 0.36835579897111215, 0.37497398012342953, 0.3815484476991695, 0.3880724695655256, 0.3945421049490535, 0.4009485948898726, 0.4072778717975746, 0.41351814729819303, 0.41966911951255403, 0.4257374513283227, 0.43173215759325456, 0.43766071784905874, 0.4435296460362384, 0.44934451665195557, 0.4551115697167973, 0.4608386561342819, 0.46652707741059735, 0.47217545640139064, 0.47778272243116404, 0.48335273038634896, 0.4888817757682528, 0.49436721548220375, 0.49981156900129603]

cspline = xopt[1:5]
tspline = xopt[6:9]
tsr = xopt[10]
idx = 10
pitch = xopt[idx+1:idx+nP]
# idx += nP
# Omega = xopt[idx+1:idx+nP]

AEP, V, T, P, Q, Np = turbineopt(cspline, tspline, tsr, pitch, "dense")

using PyPlot
close("all")

figure()
plot(V, P/1e6)
xlabel("V (m/s)")
ylabel("P (MW)")

figure()
plot(V, T)
xlabel("V (m/s)")
ylabel("T (N)")

# figure()
# plot(V, Q)

# p0 = pitch[1]
# println("p0 = ", p0*180/pi)
# pitch .-= p0

figure()
plot(V, pitch*180/pi)
xlabel("V (m/s)")
ylabel("pitch (deg)")

# figure()
# plot(V, Omega)

r, chord, theta, Rhub, Rtip = shape(cspline, tspline)

# theta .+= p0

figure()
plot(r, Np)
xlabel("r (m)")
ylabel(L"N^\prime (N/m)")

figure()
plot(r, chord)
xlabel("r (m)")
ylabel("chord (m)")

figure()
plot(r, theta*180/pi)
xlabel("r (m)")
ylabel("theta (deg)")


text:
"""
Data structure used by open-loop layer to update current high-level plan.

Attributes:
    - `curr_graph::SimpleVListGraph{OpenLoopVertex{D,C,AD}}` A list of vertices generated so far by the bridge sampler
    - `curr_start_idx::Int64` The vertex index of the current start vertex in the graph, from which planning will happen. Updated each time open_loop_plan is called
    - `curr_goal_idx::Int64` The vertex index of current goal vertex. Updated each time a goal state vertex is popped.
    - `curr_soln_path_idxs::Vector{Int64}` The in-order list of indices from current start to current goal
    - `num_samples::Int64` The N parameter for the number of bridge samples to generate per mode switch
    - `bookkeeping::B` A data structure for tracking information pertaining to the graph search
"""
mutable struct GraphTracker{D,C,AD,M,B,RNG <: AbstractRNG}
    curr_graph::SimpleVListGraph{OpenLoopVertex{D,C,AD,M}}
    curr_start_idx::Int64
    curr_goal_idx::Int64
    curr_soln_path_idxs::Vector{Int64}
    num_samples::Int64
    bookkeeping::B
    rng::RNG
end

metadatatype(::GraphTracker{D,C,AD,M,B,RNG}) where {D,C,AD,M,B,RNG <: AbstractRNG} = M
bookkeepingtype(::GraphTracker{D,C,AD,M,B,RNG}) where {D,C,AD,M,B,RNG <: AbstractRNG} = B

function GraphTracker{D,C,AD,M,B,RNG}(N::Int64, bookkeeping::B, rng::RNG=Random.GLOBAL_RNG) where {D,C,AD,M,B,RNG <: AbstractRNG}
    return GraphTracker{D,C,AD,M,B,RNG}(SimpleVListGraph{OpenLoopVertex{D,C,AD,M}}(),
                        0, 0, Vector{Int64}(undef, 0), N, bookkeeping, rng)
end

"""
Open loop plan from current state from scratch. Return a graph tracker with solution

Arguments:
    - `cmssp::CMSSP{D,C,AD,AC}` The CMSSP instance
    - `s_t::CMSSPState{D,C}` The current state
    - `curr_time::Int64` The current time step value
    - `edge_weight::Function` A generic edge weight function (handed down by hhpc solver)
    - `heuristic::Function` A generic heuristic function (handed down by hhpc solver)
    - `goal_modes::Vector{D}` A set of modes that are goal conditions
    - `graph_tracker::GraphTracker{D,C}` The instance of the object that maintains the current open-loop graph
    - `start_metadata::M` Default metadata for a new start vertex

Returns:
    Updates the `graph_tracker` object in place.
"""
function open_loop_plan!(cmssp::CMSSP, s_t::CMSSPState,
                         curr_time::Int64,
                         edge_weight::Function,
                         heuristic::Function,
                         goal_modes::Vector{D},
                         graph_tracker::GraphTracker{D,C,AD,M,B,RNG},
                         start_metadata::M) where {D,C,AD,M,B,RNG <: AbstractRNG}

    # First update the context
    update_vertices_with_context!(cmssp, graph_tracker, curr_time)

    # Create start vertex and insert in graph
    # Don't need explicit goal - visitor will handle
    add_vertex!(graph_tracker.curr_graph, OpenLoopVertex{D,C,AD,M}(s_t, cmssp.mode_actions[1], TPDistribution([curr_time], [1.0]), start_metadata))
    graph_tracker.curr_start_idx = num_vertices(graph_tracker.curr_graph)
    graph_tracker.curr_goal_idx = 0

    # Clear out current solution
    empty!(graph_tracker.curr_soln_path_idxs)

    # Obtain path and current cost with A*
    a_star_path_soln = a_star_implicit_shortest_path!(graph_tracker.curr_graph, edge_weight,
                                                      graph_tracker.curr_start_idx,
                                                      GoalVisitorImplicit(graph_tracker, cmssp, goal_modes))

    # If unreachable, report warning
    if graph_tracker.curr_goal_idx == 0
        @warn "A goal state is currently unreachable"
        return
    end

    # Check that goal state is indeed terminal
    # And that A* has non-Inf cost path to it
    @assert a_star_path_soln.dists[graph_tracker.curr_goal_idx] < Inf "Path to goal is Inf cost!"

    ## Walk path back to goal
    # Insert goal in soln vector
    pushfirst!(graph_tracker.curr_soln_path_idxs, graph_tracker.curr_goal_idx)
    curr_vertex_idx = graph_tracker.curr_goal_idx

    while curr_vertex_idx != graph_tracker.curr_start_idx
        prev_vertex_idx = a_star_path_soln.parent_indices[curr_vertex_idx]
        pushfirst!(graph_tracker.curr_soln_path_idxs, prev_vertex_idx)
        curr_vertex_idx = prev_vertex_idx
    end
    # @show graph_tracker.curr_soln_path_idxs
    # for idx in graph_tracker.curr_soln_path_idxs
    #     println(graph_tracker.curr_graph.vertices[idx].state," at ",graph_tracker.curr_graph.vertices[idx].tp)
    #     hor_val = graph_tracker.curr_graph.vertices[idx].tp.vals[1]
    #     if hor_val < typemax(Int64)
    #         display_context_future(context_set,hor_val)
    #     end
    # end
end

"""
The visitor object for the implicit A* search. Attributes are obvious.
"""
struct GoalVisitorImplicit <: AbstractDijkstraVisitor
    graph_tracker::GraphTracker
    cmssp::CMSSP
    goal_modes::Vector
end

function Graphs.include_vertex!(vis::GoalVisitorImplicit,
                                u::OpenLoopVertex, v::OpenLoopVertex, d::Float64, nbrs::Vector{Int64})

    # @show v
    # If popped vertex is terminal, then stop
    if isterminal(vis.cmssp, v.state) == true
        vis.graph_tracker.curr_goal_idx = vertex_index(vis.graph_tracker.curr_graph, v)
        return false
    end

    popped_mode = v.state.mode
    popped_cont = v.state.continuous

    # If goal mode but NOT goal state, add samples from goal
    if popped_mode in vis.goal_modes

        # Generate goal sample set
        (vertices_to_add, nbrs_to_add) = generate_goal_vertex_set!(vis.cmssp, v, vis.graph_tracker, vis.graph_tracker.rng)

        for nbr_idx in nbrs_to_add
            push!(nbrs, nbr_idx)
        end

        num_new_vertices = length(vertices_to_add)

        if num_new_vertices > 0

            # Add vertices to graph and to nbrs
            for vtx in vertices_to_add
                add_vertex!(vis.graph_tracker.curr_graph, vtx)
                push!(nbrs, num_vertices(vis.graph_tracker.curr_graph))
            end
        end
    end

    # Now add for next modes
    next_valid_modes = generate_next_valid_modes(vis.cmssp, popped_mode)

    # @show next_valid_modes

    for (action, nvm) in next_valid_modes

        # Generate bridge samples
        (vertices_to_add, nbrs_to_add) = generate_bridge_vertex_set!(vis.cmssp, v,
                                                                    (popped_mode, nvm),
                                                                    vis.graph_tracker,
                                                                    action,
                                                                    vis.graph_tracker.rng)
        for nbr_idx in nbrs_to_add
            push!(nbrs, nbr_idx)
        end

        num_new_vertices = length(vertices_to_add)
        # @show num_new_vertices

        if num_new_vertices > 0

            # Add vertices to graph and to nbrs
            for vtx in vertices_to_add
                add_vertex!(vis.graph_tracker.curr_graph, vtx)
                push!(nbrs, num_vertices(vis.graph_tracker.curr_graph))
            end
        end
    end

    # @show v
    # @show nbrs
    # readline()
    return true
end

ext: jl


text:
"""
Liu bursting neuron script
"""

using NeuronBuilder, ModelingToolkit, OrdinaryDiffEq, Plots

const Cm = 10.0 # specific capacitance cₘ is a biological constant (around) 10 nF/mm^2
Liu_conv = Cm

channels = [Liu.Na(700.0 * Liu_conv), Liu.CaS(4.0 * Liu_conv), Liu.CaT(2.0 * Liu_conv), Liu.Ka(50.0 * Liu_conv), Liu.KCa(40.0 * Liu_conv),
    Liu.Kdr(70.0 * Liu_conv), Liu.H(0.03 * Liu_conv), Liu.leak(0.01 * Liu_conv)]

τCa = 20.0
Ca∞ = 0.05
fxarea = 14.96 * 0.0628

dynamics = Dict(Voltage => BasicVoltageDynamics(),
    Calcium => Liu.CalciumDynamics(τCa, Ca∞, fxarea),
    Reversal{Calcium} => Liu.CaReversalDynamics())

somatic_parameters = Dict(
    Reversal{Sodium} => 50.0,
    Reversal{Potassium} => -80.0,
    Reversal{Leak} => -50.0,
    Reversal{Proton} => -20.0,
    Voltage => -60.0,
    Calcium => 0.05,
    Reversal{Calcium} => 0.0)

b = BasicNeuron(NoGeometry(Cm), dynamics, somatic_parameters, channels, :test_Liu) #point neuron

neur = b()

prob = ODEProblem(neur, [], (0.0, 5000.0), []; jac=true)

sol = @time solve(prob, Tsit5())
plot(sol)


text:
function TwoSum(a::Float64,b::Float64)
	x = a+b
	z = x-a
	y = (a-(x-z))+(b-z)
	return (x,y)
end

function Split(a::Float64)
	p = 52
	s = 26
	factor = 67108865
	c = factor*a
	x = c - (c-a)
	y = a - x
    return(x,y)
end

function Split(a::Float32)
	p = 23
	s = 12
	factor = 4097
	c = factor*a
	x = c - (c-a)
	y = a - x
    return (x,y)
end

function TwoProduct(a::Float64,b::Float64)
	x = a*b
	[a1,a2] = Split(a)
	[b1,b2] = Split(b)
	y = (a2*b2−(((x−a1*b1)−a2*b1)−a1*b2))
end

function TwoSumComplex(x::Complex128,y::Complex128)
	a = real(x)
	b = imag(x)
	c = real(y)
	d = imag(y)
	(s1,e1) = TwoSum(a,c)
	(s2,e2) = TwoSum(b,d)
	s = s1 + im*s2
	e = e1 + im*e2
	return (s,e)
end

function TwoProductComplex(x::Complex128,y::Complex128)
	a = real(x)
	b = imag(x)
	c = real(y)
	d = imag(y)
	(z1,h1) = TwoProduct(a,c)
	(z2,h2) = TwoProduct(b,d)
	(z3,h3) = TwoProduct(a,d)
	(z4,h4) = TwoProduct(b,c)
	(z5,h5) = TwoSum(z1,-z2)
	(z6,h6) = TwoSum(z3,z4)
	p = z5+im*z6
	e = h1+im*h3
	f = -h2+im*h4
	g = h5+im*h6
	return (p,e,f,g)
end

function TwoProductComplexFast(x::Complex128,y::Complex128)
	a = real(x)
	b = imag(x)
	c = real(y)
	d = imag(y)
	(a1,a2) = Split(a)
	(b1,b2) = Split(b)
	(c1,c2) = Split(c)
	(d1,d2) = Split(d)
	z1 = a*c
	z2 = b*d
	z3 = a*d
	z4 = b*c
	h1 = (a2*c2-(((z1-a1*c1)-a2*c1)-a1*c2))
	h2 = (b2*d2-(((z2-b1*d1)-b2*d1)-b1*d2))
	h3 = (a2*d2-(((z3-a1*d1)-a2*d1)-a1*d2))
	h4 = (b2*c2-(((z4-b1*c1)-b2*c1)-b1*c2))
	(z5,h5) = TwoSum(z1,-z2)
	(z6,h6) = TwoSum(z3,z4)
	p =  z5 + im*z6
	e =  h1 + im*h3
	f = -h2 + im*h4
	g =  h5 + im*h6
	return (p,e,f,g)
end


text:
## header to provide surface and context
using Cairo
c = CairoARGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.5,0.5,0.5);
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);
save(cr);

## paint the following to a Recording Surface
s1 = CairoRecordingSurface()
c1 = CairoContext(s1)

move_to(c1, 128.0, 25.6);
line_to(c1, 230.4, 230.4);
rel_line_to(c1, -102.4, 0.0);
curve_to(c1, 51.2, 230.4, 51.2, 128.0, 128.0, 128.0);
close_path(c1);

move_to(c1, 64.0, 25.6);
rel_line_to(c1, 51.2, 51.2);
rel_line_to(c1, -51.2, 51.2);
rel_line_to(c1, -51.2, -51.2);
close_path(c1);

set_line_width(c1, 10.0);
set_source_rgb(c1, 0, 0, 1);
fill_preserve(c1);
set_source_rgb(c1, 0, 0, 0);
stroke(c1);

## write to script

io = IOBuffer()
s2 = Cairo.CairoScript(io)
Cairo.script_from_recording_surface(s2,s1);
destroy(s2)

seek(io,0)
r = String(read(io))

r1 = split(r,"\n")

t1 = 24.0
for x in r1
    move_to(cr,0.0,t1)
    show_text(cr,x)
    t1 += 12.0
end

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_script0.png");


text:
conic_form!(s::Positive, x::AbstractVariable, unique_conic_forms) = conic_form!(x>=0, unique_conic_forms)
conic_form!(s::Negative, x::AbstractVariable, unique_conic_forms) = conic_form!(x<=0, unique_conic_forms)


text:
@testset "discreteshapebound" begin
    img = generate_surface(SynthSphere(50), 0.5, [0.5,0.0,0.9])
    Z, p, q = DiscreteShapeBound(smoothness=1000, albedo = 1.0, illumination_direction = [0.5,0.0,0.9])(img)
    @test Z[38,38] ≈ 0.00028515016913445543
    @test Z[1,1] ≈ 9.035191974767473e-5
    @test isapprox(Z[76,76], 0.00034007074121476805; atol = 2e-15)
    @test Z[45,45] ≈ 0.00010731186855283559
    @test Z[50,18] ≈ 0.0003244761469170222
    @test Z[11, 68] ≈ 7.131234669714447e-5
    @test Z[42,42] ≈ 0.0002467784663847426
end


text:
using DelimitedFiles
using MLDataUtils

function traintest()
  a = readdlm("data/travel.txt", ',')
  a = shuffleobs(a, obsdim=1)
  train, test = splitobs(a, at = 0.7, obsdim=1)
  train, test
  open("data/train.txt", "w") do f
    writedlm(f, train, ",")
  end
  open("data/test.txt", "w") do f
    writedlm(f, test, ",")
  end
end
traintest()


text:
##############################################################################
##
## Iteration: eachrow, eachcol
##
##############################################################################

# Iteration by rows
"""
    DataFrameRows{D<:AbstractDataFrame} <: AbstractVector{DataFrameRow{D,S}}

Iterator over rows of an `AbstractDataFrame`,
with each row represented as a `DataFrameRow`.

A value of this type is returned by the [`eachrow`](@ref) function.
"""
struct DataFrameRows{D<:AbstractDataFrame,S} <: AbstractVector{DataFrameRow{D,S}}
    df::D
end

Base.summary(dfrs::DataFrameRows) = "$(length(dfrs))-element DataFrameRows"
Base.summary(io::IO, dfrs::DataFrameRows) = print(io, summary(dfrs))

"""
    eachrow(df::AbstractDataFrame)

Return a `DataFrameRows` that iterates a data frame row by row,
with each row represented as a `DataFrameRow`.

# Examples
```jldoctest
julia> df = DataFrame(x=1:4, y=11:14)
4×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 11    │
│ 2   │ 2     │ 12    │
│ 3   │ 3     │ 13    │
│ 4   │ 4     │ 14    │

julia> eachrow(df)
4-element DataFrameRows:
 DataFrameRow (row 1)
x  1
y  11
 DataFrameRow (row 2)
x  2
y  12
 DataFrameRow (row 3)
x  3
y  13
 DataFrameRow (row 4)
x  4
y  14

julia> copy.(eachrow(df))
4-element Array{NamedTuple{(:x, :y),Tuple{Int64,Int64}},1}:
 (x = 1, y = 11)
 (x = 2, y = 12)
 (x = 3, y = 13)
 (x = 4, y = 14)

julia> eachrow(view(df, [4,3], [2,1]))
2-element DataFrameRows:
 DataFrameRow (row 4)
y  14
x  4
 DataFrameRow (row 3)
y  13
x  3
```
"""
eachrow(df::AbstractDataFrame) = DataFrameRows{typeof(df), typeof(index(df))}(df)

Base.IndexStyle(::Type{<:DataFrameRows}) = Base.IndexLinear()
Base.size(itr::DataFrameRows) = (size(parent(itr), 1), )

Base.@propagate_inbounds function Base.getindex(itr::DataFrameRows, i::Int)
    df = parent(itr)
    DataFrameRow(df, index(df), i)
end

Base.@propagate_inbounds function Base.getindex(itr::DataFrameRows{<:SubDataFrame}, i::Int)
    sdf = parent(itr)
    DataFrameRow(parent(sdf), index(sdf), rows(sdf)[i])
end

Base.getproperty(itr::DataFrameRows, col_ind::Symbol) = getproperty(parent(itr), col_ind)
# Private fields are never exposed since they can conflict with column names
Base.propertynames(itr::DataFrameRows, private::Bool=false) = propertynames(parent(itr))

# Iteration by columns

"""
    DataFrameColumns{<:AbstractDataFrame} <: AbstractVector{AbstractVector}

An `AbstractVector` that allows iteration over columns of an `AbstractDataFrame`.
Indexing into `DataFrameColumns` objects using integer or symbol indices
returns the corresponding column (without copying).
"""
struct DataFrameColumns{T<:AbstractDataFrame} <: AbstractVector{AbstractVector}
    df::T
end

Base.summary(dfcs::DataFrameColumns)= "$(length(dfcs))-element DataFrameColumns"
Base.summary(io::IO, dfcs::DataFrameColumns) = print(io, summary(dfcs))

"""
    eachcol(df::AbstractDataFrame)

Return a `DataFrameColumns` that is an `AbstractVector`
that allows iterating an `AbstractDataFrame` column by column.
Additionally it is allowed to index `DataFrameColumns` using column names.

# Examples
```jldoctest
julia> df = DataFrame(x=1:4, y=11:14)
4×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 11    │
│ 2   │ 2     │ 12    │
│ 3   │ 3     │ 13    │
│ 4   │ 4     │ 14    │

julia> collect(eachcol(df))
2-element Array{AbstractArray{T,1} where T,1}:
 [1, 2, 3, 4]
 [11, 12, 13, 14]

julia> map(eachcol(df)) do col
           maximum(col) - minimum(col)
       end
2-element Array{Int64,1}:
 3
 3

julia> sum.(eachcol(df))
2-element Array{Int64,1}:
 10
 50
```
"""
eachcol(df::AbstractDataFrame) = DataFrameColumns(df)

Base.size(itr::DataFrameColumns) = (size(parent(itr), 2),)
Base.IndexStyle(::Type{<:DataFrameColumns}) = Base.IndexLinear()

@inline function Base.getindex(itr::DataFrameColumns, j::Int)
    @boundscheck checkbounds(itr, j)
    @inbounds parent(itr)[!, j]
end

Base.getindex(itr::DataFrameColumns, j::Symbol) = parent(itr)[!, j]

Base.getproperty(itr::DataFrameColumns, col_ind::Symbol) = getproperty(parent(itr), col_ind)
# Private fields are never exposed since they can conflict with column names
Base.propertynames(itr::DataFrameColumns, private::Bool=false) = propertynames(parent(itr))

"""
    keys(dfc::DataFrameColumns)

Get a vector of column names of `dfc`.
"""
Base.keys(itr::DataFrameColumns) = names(parent(itr))

"""
    pairs(dfc::DataFrameColumns)

Return an iterator of pairs associating the name of each column of `dfc`
with the corresponding column vector, i.e. `name => col`
where `name` is the column name of the column `col`.
"""
Base.pairs(itr::DataFrameColumns) = Base.Iterators.Pairs(itr, keys(itr))

"""
    mapcols(f::Union{Function,Type}, df::AbstractDataFrame)

Return a `DataFrame` where each column of `df` is transformed using function `f`.
`f` must return `AbstractVector` objects all with the same length or scalars.

Note that `mapcols` guarantees not to reuse the columns from `df` in the returned
`DataFrame`. If `f` returns its argument then it gets copied before being stored.

# Examples
```jldoctest
julia> df = DataFrame(x=1:4, y=11:14)
4×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 11    │
│ 2   │ 2     │ 12    │
│ 3   │ 3     │ 13    │
│ 4   │ 4     │ 14    │

julia> mapcols(x -> x.^2, df)
4×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 121   │
│ 2   │ 4     │ 144   │
│ 3   │ 9     │ 169   │
│ 4   │ 16    │ 196   │
```
"""
function mapcols(f::Union{Function,Type}, df::AbstractDataFrame)
    # note: `f` must return a consistent length
    vs = AbstractVector[]
    seenscalar = false
    seenvector = false
    for v in eachcol(df)
        fv = f(v)
        if fv isa AbstractVector
            if seenscalar
                throw(ArgumentError("mixing scalars and vectors in mapcols not allowed"))
            end
            seenvector = true
            push!(vs, fv === v ? copy(fv) : fv)
        else
            if seenvector
                throw(ArgumentError("mixing scalars and vectors in mapcols not allowed"))
            end
            seenscalar = true
            push!(vs, [fv])
        end
    end
    DataFrame(vs, _names(df), copycols=false)
end

Base.parent(itr::Union{DataFrameRows, DataFrameColumns}) = getfield(itr, :df)
Base.names(itr::Union{DataFrameRows, DataFrameColumns}) = names(parent(itr))

function Base.show(io::IO, dfrs::DataFrameRows;
                   allrows::Bool = !get(io, :limit, false),
                   allcols::Bool = !get(io, :limit, false),
                   splitcols = get(io, :limit, false),
                   rowlabel::Symbol = :Row,
                   summary::Bool = true)
    df = parent(dfrs)
    summary && print(io, "$(nrow(df))×$(ncol(df)) DataFrameRows")
    _show(io, df, allrows=allrows, allcols=allcols, splitcols=splitcols,
          rowlabel=rowlabel, summary=false)
end

Base.show(io::IO, mime::MIME"text/plain", dfrs::DataFrameRows;
          allrows::Bool = !get(io, :limit, false),
          allcols::Bool = !get(io, :limit, false),
          splitcols = get(io, :limit, false),
          rowlabel::Symbol = :Row,
          summary::Bool = true) =
    show(io, dfrs, allrows=allrows, allcols=allcols, splitcols=splitcols,
         rowlabel=rowlabel, summary=summary)

Base.show(dfrs::DataFrameRows;
          allrows::Bool = !get(stdout, :limit, true),
          allcols::Bool = !get(stdout, :limit, true),
          splitcols = get(stdout, :limit, true),
          rowlabel::Symbol = :Row,
          summary::Bool = true) =
    show(stdout, dfrs, allrows=allrows, allcols=allcols, splitcols=splitcols,
         rowlabel=rowlabel, summary=summary)

function Base.show(io::IO, dfcs::DataFrameColumns;
                   allrows::Bool = !get(io, :limit, false),
                   allcols::Bool = !get(io, :limit, false),
                   splitcols = get(io, :limit, false),
                   rowlabel::Symbol = :Row,
                   summary::Bool = true)
    df = parent(dfcs)
    summary && print(io, "$(nrow(df))×$(ncol(df)) DataFrameColumns")
    _show(io, parent(dfcs), allrows=allrows, allcols=allcols, splitcols=splitcols,
          rowlabel=rowlabel, summary=false)
end

Base.show(io::IO, mime::MIME"text/plain", dfcs::DataFrameColumns;
          allrows::Bool = !get(io, :limit, false),
          allcols::Bool = !get(io, :limit, false),
          splitcols = get(io, :limit, false),
          rowlabel::Symbol = :Row,
          summary::Bool = true) =
    show(io, dfcs, allrows=allrows, allcols=allcols, splitcols=splitcols,
         rowlabel=rowlabel, summary=summary)

Base.show(dfcs::DataFrameColumns;
          allrows::Bool = !get(stdout, :limit, true),
          allcols::Bool = !get(stdout, :limit, true),
          splitcols = get(stdout, :limit, true),
          rowlabel::Symbol = :Row,
          summary::Bool = true) =
    show(stdout, dfcs, allrows=allrows, allcols=allcols, splitcols=splitcols,
         rowlabel=rowlabel, summary=summary)


text:
module CommonFunctions

    export nbasis, nnodes, nodes, order, degree

    nbasis() = nothing
    nnodes() = nothing
    nodes()  = nothing
    order()  = nothing
    degree() = nothing

    export evaluate, evaluate!

    evaluate() = nothing
    evaluate!() = nothing

    export derivative, integral

    derivative() = nothing
    integral() = nothing

    export write_to_hdf5

    write_to_hdf5() = nothing

    export get_solution, get_solution!, set_solution!, reset!, cut_periodic_solution!

    get_solution() = nothing
    get_solution!() = nothing
    set_solution!() = nothing
    reset!() = nothing
    cut_periodic_solution!() = nothing

end


text:
using Test
using LogicCircuits

function little_2var()
    lin = LogicalΔNode[]

    v = Var(2)
    pos = LiteralNode( var2lit(v))
    push!(lin, pos)
    neg = LiteralNode(-var2lit(v))
    push!(lin, neg)
    or1 = ⋁Node([pos,neg])
    push!(lin, or1)
    or2 = ⋁Node([pos,neg])
    push!(lin, or2)

    v = Var(1)
    pos = LiteralNode( var2lit(v))
    push!(lin, pos)
    neg = LiteralNode(-var2lit(v))
    push!(lin, neg)

    and1 = ⋀Node([pos, or1])
    and2 = ⋀Node([neg, or2])
    root = ⋁Node([and1, and2])

    push!(lin, and1)
    push!(lin, and2)
    push!(lin, root)
end

function little_3var()
    lin = little_2var()
    or1 = lin[end]
    v = Var(3)

    pos = LiteralNode( var2lit(v))
    push!(lin, pos)
    neg = LiteralNode(-var2lit(v))
    push!(lin, neg)

    or2 = ⋁Node(children(or1))
    push!(lin, or2)

    and1 = ⋀Node([pos, or1])
    and2 = ⋀Node([neg, or2])
    root = ⋁Node([and1,and2])
    push!(lin, and1)
    push!(lin, and2)
    push!(lin, root)
    lin
end

function little_4var()
    lin = LogicalΔNode[]
    n = 4
    ors = map(1:n) do v
        v = Var(v)
        pos = LiteralNode( var2lit(v))
        push!(lin, pos)
        neg = LiteralNode(-var2lit(v))
        push!(lin, neg)
        or = ⋁Node([pos,neg])
        push!(lin, or)
        or
    end
    and1 = ⋀Node(ors[1:2])
    and2 = ⋀Node(ors[3:4])
    or = ⋁Node([and1, and2])
    push!(lin, and1)
    push!(lin, and2)
    push!(lin, or)
    lin
end

function little_5var()
    lin = little_4var()
    v = Var(5)

    pos = LiteralNode( var2lit(v))
    push!(lin, pos)
    neg = LiteralNode(-var2lit(v))
    push!(lin, neg)
    or = ⋁Node([pos,neg])
    push!(lin, or)

    and = ⋀Node([lin[15], lin[18]])
    push!(lin, and)
    or = ⋁Node([lin[19]])
    push!(lin, or)
    lin
end

@testset "Condition test" begin
    c1 = load_logical_circuit(zoo_sdd_file("random.sdd"))[end]

    lit = Lit(num_variables(c1) + 1)
    @test c1 == condition(c1, lit)

    lit1 = Lit(1)
    c2 = condition(c1, lit1)
    dict = Logical.literal_nodes(c2)
    @test haskey(dict, lit1)
    @test !haskey(dict, -lit1)
    c3 = condition(c2, -lit1)
    @test is_false(c3)
    c4 = condition(c2, lit1)
    @test c4 == c2

    #
    c1 = little_2var()
    c2 = condition(c1, lit1)
    @test num_nodes(c2) == 6
    @test num_edges(c2) == 5

    lit2 = Lit(2)
    c3 = condition(c1, lit2)
    @test num_nodes(c3) == 6
    @test num_edges(c3) == 6

    c4 = condition(c2, lit2)
    @test num_nodes(c4) == 4
    @test num_edges(c4) == 3
end

@testset "Split test" begin
    c0 = little_5var()
    or = c0[end]
    and = children(or)[1]
    v = Var(1)
    c1, _ = split(c0, (or, and), v)
    @test num_nodes(c1) == 22
    @test num_edges(c1) == 24
    @test c1[end].children[1].children[2] == c1[end].children[2].children[2]
    @test c1[end].children[1].children[1].children[2] == c1[end].children[2].children[1].children[2]

    or = c1[end].children[1].children[1]
    and = children(or)[2]
    v = Var(4)
    c2, _ = split(c1, (or, and), v)
    @test num_nodes(c2) == 24
    @test num_edges(c2) == 29
    n1 = c2[end].children[2].children[1].children[2]
    n2 = c2[end].children[1].children[1].children[2]
    @test n1 in c1
    @test !(n2 in c1)
    @test n1.children[1] == n2.children[1]
end

@testset "Copy test" begin
    n0 = little_2var()[end]
    n1 = copy(n0, 0)
    @test n0 === n1
    n2 = copy(n0, 1)
    @test n2 !== n0
    @test all(children(n2) .=== children(n0))

    for depth in [5, 10, 100, 1000]
        n3 = load_logical_circuit(zoo_sdd_file("random.sdd"))[end]
        n4 = copy(n3, depth)
        @test num_nodes(n3) == num_nodes(n4)
        @test num_edges(n3) == num_edges(n4)
    end

end

@testset "Clone test" begin
    n0 = little_3var()[end]
    and = n0.children[1].children[2].children[1]
    or1 = n0.children[1].children[2]
    or2 = n0.children[2].children[2]
    n1 = clone(n0, or1, or2, and; depth=1)
    @test num_nodes(n1) == num_nodes(n0) + 1
    @test num_edges(n1) == num_edges(n0) + num_children(and)
    @test length(Set([node2dag(n1); node2dag(n0)])) == num_nodes(n0) + 4
    n2 = clone(n0, or1, or2, and; depth=2)
    @test num_nodes(n2) == num_nodes(n0) + 2
    @test num_edges(n2) == num_edges(n0) + 4
    @test length(Set([node2dag(n2); node2dag(n0)])) == num_nodes(n0) + 5
end


text:
#########
#Created Date: Wednesday April 28th 2021
#Author: Domenico Di Gangi,  <digangidomenico@gmail.com>
#-----
#Last Modified: Tuesday June 1st 2021 7:55:51 pm
#Modified By:  Domenico Di Gangi
#-----
#Description:
#-----
########

#region import and models
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using DrWatson
using JLD2
using Distributed
using SharedArrays
using ScoreDrivenERGM
using Logging

begin
nWorkers = 10
addprocs(nWorkers - nprocs())
@sync @everywhere begin
    using Pkg
    Pkg.activate(".")

    using ScoreDrivenERGM
    import ScoreDrivenERGM:StaticNets, DynNets
    import ScoreDrivenERGM.DynNets:SdErgm,SdErgmDirBin0Rec0, simulate_and_estimate_parallel
    using ScoreDrivenERGM.Utilities

    end
end

#endregion

model = DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges +  gwesp(decay = 0.25,fixed = TRUE)", true), indTvPar = [true, true], scoreScalingType="FISH_D")

dgpSettings = (type = "AR", opt = (B =[1], sigma = [0, 0.18],  minVals = [-3.0, 0], maxVals = [-0.1, 2] ))

sigmaVals = [0.005, 0.1, 0.2]
dgpList = [deepcopy(dgpSettings) for i in 1:length(sigmaVals)]
[dgpList[ind].opt.sigma[2] = sigma for (ind, sigma) in enumerate(sigmaVals)]

#region simulations
c= Dict{String, Any}()
c["model"] =[ DynNets.SdErgmPml(staticModel = StaticNets.NetModeErgmPml("edges +  gwesp(decay = 0.25,fixed = TRUE)", true), indTvPar = [false, true], scoreScalingType="FISH_D")]
c["T"] = [73]
c["N"] = [100]
c["dgpSettings"] = dgpList
c["nSample"] = 50

list = sort(dict_list(c), by=x->(string(x["model"])))

for d in list
    nSample = d["nSample"]
    T = d["N"]
    timeSim = @elapsed allObsT, allvEstSdResPar, allfVecT_filt, allParDgpT, allConvFlag, allftot_0,  allfVecT_filt_SS =  ScoreDrivenERGM.DynNets.simulate_and_estimate_parallel(d["model"], d["dgpSettings"], d["T"], d["N"],  d["nSample"];)

    modelTag = string(d["model"])

    res1 =  (;modelTag, allObsT, allvEstSdResPar, allfVecT_filt, allParDgpT, allConvFlag, allftot_0, allfVecT_filt_SS) |>DrWatson.ntuple2dict |> DrWatson.tostringdict

    estDict = merge(res1, d)

    saveName = replace.( savename(d;allowedtypes = (Real, String, Symbol, NamedTuple, Tuple, ScoreDrivenERGM.DynNets.SdErgm) ), r"[\"]" => "")

    estDict["saveName"] = saveName

    timeSave = @elapsed save( datadir("sims", "dgpAR_Fil_SS", string(hash(saveName))) *".jld2", estDict)

    Logging.@info("Time sim = $timeSim ,  time save = $timeSave ")

end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
using ParallelAccelerator

@acc function kernel_doitgen(nr::Int, nq::Int, np::Int,
    A::Array{Float64, 3}, C4::Matrix{Float64})

    sum = Array{Float64}(np)
    for r = 1:nr
        for q = 1:nq
            for p = 1:np
                sum[p] = 0.0
                for s = 1:np
                    sum[p] += A[r, q, s] * C4[s, p]
                end
            end
            for p = 1:np
                A[r, q, p] = sum[p]
            end
        end
    end
end

@inbounds function kernel_doitgen_o(nr::Int, nq::Int, np::Int,
    A::Array{Float64, 3}, C4::Matrix{Float64})

    tic()
    sum = Array{Float64}(np)
    for r = 1:nr
        for q = 1:nq
            for p = 1:np
                sum[p] = 0.0
                for s = 1:np
                    sum[p] += A[r, q, s] * C4[s, p]
                end
            end
            for p = 1:np
                A[r, q, p] = sum[p]
            end
        end
    end
    println("Julia $(toq())")
end

nr = 300
nq = nr
np = nr

At = rand(nr, nq, np)
C4 = rand(np, np)

Ao = copy(At)

kernel_doitgen(nr, nq, np, At, C4)
kernel_doitgen_o(nr, nq, np, Ao, C4)

function compare(A, B)
    for i in 1:nr
        for j in 1:nq
            for k in 1:np
                if A[i, j, k] != B[i, j, k]
                    println("different $((i, j, k)) : $(A[i, j, k]) != $(B[i, j, k])")
                    return false
                end
            end
        end
    end
    return true
end

println(compare(At, Ao))

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
###################################################
## geo/data/nyctaxi.jl
## Everything specific to the NYC taxi trip data
###################################################

"""
    `fromNYCTaxiData` : extract one month of trip information from original trip CSV
    Return GeoData object
"""
function fromNYCTaxiCSV(fileName::AbstractString)
    println("Counting lines...")
    f = open(fileName)
    NLINES = countlines(f) - 1
    close(f)
    println("$NLINES trips to parse")

    trips = Array{GeoTrip}(undef, NLINES)
    dateFormat = DateFormat("y-m-d H:M:S")
    f = open(fileName)
    names = split(strip(readline(f)),",")
    PLON  = findfirst(isequal("pickup_longitude"), names)
    PLAT  = findfirst(isequal("pickup_latitude"), names)
    DLON  = findfirst(isequal("dropoff_longitude"), names)
    DLAT  = findfirst(isequal("dropoff_latitude"), names)
    PTIME = findfirst(isequal("tpep_pickup_datetime"), names)
    DTIME = findfirst(isequal("tpep_dropoff_datetime"), names)

    println("Beginning trip parsing...")
    for (i,ln) in enumerate(eachline(f))
        if i%10_000 == 0
            @printf("\r%.2f%% trips parsed     ",100*i/NLINES)
        end
        s = split(strip(ln),",")
        pTime = DateTime(s[PTIME], dateFormat)
        dTime = DateTime(s[DTIME], dateFormat)

        trips[i] = GeoTrip(
            parse(Float32,s[PLON]),
            parse(Float32,s[PLAT]),
            parse(Float32,s[DLON]),
            parse(Float32,s[DLAT]),
            pTime,
            ((dTime-pTime).value)/1000.f0
        )
    end
    print("\r100.00% trips parsed     ")
    return trips
end


text:
""" Assertions on the variables used in the program. Helps locating errors right away """

# All functions are public

# All debug functions run if the global parameter CONFIG.debug is set to true.
# These functions are expensive. When launching benchmarks, make sure to set
# CONFIG.debug to false.

" `Assertion` : `correct`, `feasible`, `integer`, `sorted z1 increasing`, `unicity`, `non dominated solutions`"
function DEBUG_LB(prob::BiOKP, LB::LowerBound)

    # Change type of LB.sols
    if typeof(LB.sols) != Vector{Sol}
        # convert into array (easier to work with indices)
        sols = Vector{Sol}()
        for sol in LB.sols
            push!(sols, sol)
        end
    else
        sols = LB.sols
    end

    for i in 1:length(sols)

        # feasible
        DEBUG_feasibleBinarySolution(prob, sols[i])

        for j in (i+1):length(sols)

            # unicity
            if prob.isInteger
                @assert sols[i].y != sols[j].y "The LB contains duplicates\n$sols"
            else
                @assert differentFloatSols(sols[i],sols[j]) "The LB contains duplicates ($(arrayLB[i].y) and $(arrayLB[j])) in : \n$sols"
            end

            # non dominated solutions
            @assert !dominate(sols[i],sols[j]) "Solutions are not all nondominated, $(sols[i]) dominates $(sols[j])\n$sols"
            @assert !dominate(sols[j],sols[i]) "Solutions are not all nondominated, $(sols[j]) dominates $(sols[i])\n$sols"

            # sorted z1 increasing
            @assert sols[i].y[1] < sols[j].y[1] "The LB isn't sorted !\n$(sols)"

        end
    end
end

" `Assertion` : `correct`, `feasible`, `sorted z1 increasing`, `unicity`, `non dominated solutions` "
function DEBUG_UB(prob::BiOKP, UB::UpperBound)

    # Change type of UB.sols
    if typeof(UB.sols) != Vector{Sol}
        # convert into array (easier to work with indices)
        sols = Vector{Sol}()
        for sol in UB.sols
            push!(sols, sol)
        end
    else
        sols = UB.sols
    end

    # check solutions
    for i in 1:length(sols)

        # correct
        DEBUG_feasibleSolution(prob, sols[i])

        for j in (i+1):length(sols)

            # unicity
            if prob.isInteger
                @assert sols[i].y != sols[j].y "The UB contains duplicates\n$sols"
            else
                @assert differentFloatSols(sols[i],sols[j]) "The UB contains duplicates ($(arrayLB[i].y) and $(arrayLB[j])) in : \n$sols"
            end

            # non dominated solutions
            @assert !dominate(sols[i],sols[j]) "Solutions are not all nondominated, $(sols[i]) dominates $(sols[j])\n$sols"
            @assert !dominate(sols[j],sols[i]) "Solutions are not all nondominated, $(sols[j]) dominates $(sols[i])\n$sols"

            # sorted z1 increasing
            @assert sols[i].y[1] < sols[j].y[1] "The UB isn't sorted !\n$(sols)"

        end
    end

    # check segments
    if UB.explicit
        @assert (length(UB.segments.b) == 0 && length(sols) < 1) || length(UB.segments.b) == length(sols)+1 "Incorrect number of segments : $(length(UB.segments.b)) for a $(length(sols)) solutions"
        @assert (size(UB.segments.A)[1] == 0 && length(sols) < 1) || size(UB.segments.A)[1] == length(sols)+1 "Incorrect number of segments : $(size(UB.segments.A)[1]) for a $(size(UB.segments.A)[1]) solutions"

        nbConstraints = length(sols)+1
        A = zeros(Float64, nbConstraints, 2)
        b = zeros(Float64, nbConstraints)

        max1 = sols[1].y[1]
        max2 = sols[1].y[2]

        for iter in 1:length(sols)-1
            solL = sols[iter]
            solR = sols[iter+1]

            A[iter, 1] = solL.y[2] - solR.y[2]
            A[iter, 2] = solR.y[1] - solL.y[1]

            b[iter] = (solL.y[2] - solR.y[2]) * solR.y[1] + (solR.y[1] - solL.y[1]) * solR.y[2]

            if solL.y[1] > max1
                max1 = solL.y[1]
            end
            if solR.y[1] > max1
                max1 = solR.y[1]
            end
            if solL.y[2] > max2
                max2 = solL.y[2]
            end
            if solR.y[2] > max2
                max2 = solR.y[2]
            end
        end

        # add the constraint extreme horizontal and vertical to finish the polytope

        # vertical (on the right)
        A[nbConstraints-1, 1] = 1
        A[nbConstraints-1, 2] = 0

        b[nbConstraints-1] = max1

        # horizontal (on the left)
        A[nbConstraints, 1] = 0
        A[nbConstraints, 2] = 1

        b[nbConstraints] = max2

        @assert A == UB.segments.A "Matrix A in segments is wrong : $(UB.segments.A) vs $(A)"
        @assert b == UB.segments.b "Vector b in segments is wrong : $(UB.segments.b) vs $(b)"
    end

end

" `Assertion` : `correct` "
function DEBUG_correctSolution(prob::BiOKP, sol::Sol)

    trueProfit = [0.0,0.0]
    trueWeight = 0.0
    floatNumber = false

    for i in 1:prob.nbVar
        elt = sol.x[i]

        if sol.isBinary
            @assert (elt != 0.0 || elt != 1.0) "Solution is not truly binary : $(sol.x), $(sol.isBinary)"
        else
            if elt != 0.0 && elt != 1.0
                floatNumber = true
            end
        end

        if sol.x[i] != 0.0
            trueWeight += sol.x[i] * round(Int8,prob.weights[i])
            trueProfit[1] += sol.x[i] * prob.profits[1,i]
            trueProfit[2] += sol.x[i] * prob.profits[2,i]
        end
    end

    sol.isBinary && @assert trueWeight == sol.w "Solution's weight is wrong : true = $trueWeight, actual = $(sol.w)"
    sol.isBinary && @assert trueProfit == sol.y "Solution's profit is wrong : true = $trueProfit, actual = $(sol.y)"
    !sol.isBinary && @assert floatNumber "Solution is TRULY binary : $floatNumber $(sol.x), $(sol.isBinary)"
end

" `Assertion` : `correct`, `feasible`, not always `binary` !"
function DEBUG_feasibleSolution(prob::BiOKP, sol::Sol)

    DEBUG_correctSolution(prob, sol)

    @assert sol.w <= prob.maxWeight "Solution is too heavy (trueWeight = $trueWeight), $(sol.w) for the problem $(prob.maxWeight)"

end

" `Assertion` : `correct`, `feasible`, `binary` "
function DEBUG_feasibleBinarySolution(prob::BiOKP, sol::Sol)
    DEBUG_feasibleSolution(prob, sol)

    @assert sol.isBinary "Solution is not binary ! $sol"
end

" `Assertion` : make sure the given `sol` isn't dominating any of the solutions in `sols` "
function DEBUG_isntDominatingAny(sol::Sol, sols::LinkedList{Sol})
    # convert into array (easier to work with indices)
    vec = Vector{Sol}()
    for sol in sols
        push!(vec, sol)
    end

    for i in 1:length(vec)
        @assert !dominate(sol, vec[i]) "$sol is dominating $(vec[i])"
    end
end

" `Assertion` : `correct`, `feasible`, `sorted increasing`, `unicity`, `non dominated solutions` "
function DEBUG_front(prob::BiOKP, front::LinkedList{Sol})
    # change linkedlist into array
    front_vec = Vector{Sol}()
    for sol in front
        push!(front_vec, sol)
    end

    DEBUG_front(prob, front_vec)
end

" `Assertion` : `correct`, `feasible`, `sorted increasing`, `unicity`, `non dominated solutions` "
function DEBUG_front(prob::BiOKP, front::Vector{Sol})

    for i in 1:length(front)

        DEBUG_feasibleSolution(prob, front[i])

        for j in (i+1):length(front)

            # unicity
            if prob.isInteger
                @assert front[i].y != front[j].y "The front contains duplicates\n$front"
            else
                @assert differentFloatSols(front[i],front[j]) "The front contains duplicates ($(arrayLB[i].y) and $(arrayLB[j])) in : \n$front"
            end

            # non dominated solutions
            @assert !dominate(front[i],front[j]) "Solutions are not all nondominated, $(front[i]) dominates $(front[j])\n$front"
            @assert !dominate(front[j],front[i]) "Solutions are not all nondominated, $(front[j]) dominates $(front[i])\n$front"

            # sorted
            @assert front[i].y[1] < front[j].y[1] "The lower bound isn't sorted !\n$(front)"

        end
    end
end

" `Assertion` : `correct`, `feasible`, `binary`, `left.z1 < right.z1`, `sorted increasing`, `unicity` "
function DEBUG_nadirs(prob::BiOKP,NP::Vector{PairOfSolution})

    for i in 1:(length(NP)-1)

        # feasible
        DEBUG_feasibleBinarySolution(prob, NP[i].solL)

        # pair correct : left.z1 < right.z1
        @assert NP[i].solL.y[1] < NP[i].solR.y[1] && NP[i].solL.y[2] > NP[i].solR.y[2] "Pair of nadir points is wrong\nhere : $(NP[i])\n$NP"

        # sorted and unicity
        @assert NP[i].solL.y[1] < NP[i+1].solL.y[1] "List of nadirs not sorted ! \n$NP"
    end

end

" `Assertion` : `correct`, `feasible`, `left.z1 < right.z1`, `sorted increasing`, `unicity` "
function DEBUG_pairsOfSolutions(prob::BiOKP, sols::Vector{PairOfSolution})

    for i in 1:length(sols)
        # feasible
        DEBUG_feasibleSolution(prob, sols[i].solL)
        DEBUG_feasibleSolution(prob, sols[i].solR)

        # pair correct : left.z1 < right.z1
        @assert sols[i].solL.y[1] < sols[i].solR.y[1] && sols[i].solL.y[2] > sols[i].solR.y[2] "Pair of points is wrong\nhere : $(sols[i])\n$sols"
    end

    for i in 1:length(sols)-1
        # sorted and unicity
        @assert sols[i].solL.y[1] < sols[i+1].solL.y[1] "List of points not sorted ! \n$sols"
    end
end

"""
Assertion on knownSols : `front`, `sol.x[index] == value for every sol in list`.

Assertion on previousIndicies : `sorted increasing`, `unicity`.
"""
function DEBUG_parentToChild(prob::BiOKP, PTC::ParentToChild; index = -1, value = -1)

    # secure knownSols
    DEBUG_front(prob, PTC.knownSols)

    if index != -1 && value != -1
        for elt in PTC.knownSols
            @assert elt.x[index] == value "The list doesn't contain the right elements"
        end
    end

    lenPrev = length(PTC.previousIndicies)

    @assert lenPrev == length(PTC.knownSols) "Known sols and previous indicies aren't of the same length"

    # secure previousIndicies
    prevArray = Vector{Int}(undef,lenPrev)
    iter = 1
    for elt in PTC.previousIndicies
        prevArray[iter] = elt
        iter += 1
    end

    for i in 1:lenPrev
        for j in (i+1):lenPrev
            @assert prevArray[i] != prevArray[j] "No unicity"
            @assert prevArray[i] < prevArray[j] "Indicies not sorted"
        end
    end
end


text:
using REALMSupport.RegisterCore: CenterIndexedArray

@testset "Translation function test" begin
    outer = reshape(1:480, 24, 20)
    A = outer[1:16, 3:18]
    B = outer[3:18, 2:17]
    C = outer[5:20, 1:16]

    mm = CenterIndexedArray(zeros(11, 11, 3))
    mm[:, :, -1] = translate(A, A, (5, 5))
    mm[:, :, 0] = translate(A, B, (5, 5))
    mm[:, :, 1] = translate(A, C, (5, 5))

    move = [0 1 2; 0 0 0]
    matrix = Float64[1 0 0 1]
    result = translate_optim(mm, matrix, move; g_abstol = 1e-14)

    @test minimum(result) < 1e-3
    #display(reshape(Optim.minimizer(result), 2, 2))

    #3d translation test
    outer = reshape(1:9600, 24, 20, 20)
    A = outer[3:16, 7:12, 10:14]
    B = outer[5:18, 6:11, 10:14]
    C = outer[7:20, 5:10, 10:14]
    D = outer[3:16, 6:11, 10:14]
    E = outer[3:16, 5:10, 10:14]

    mm = CenterIndexedArray(zeros(7, 7, 3, 5))
    mm[:, :, :, -2] = translate(A, A, (3, 3, 1))
    mm[:, :, :, -1] = translate(A, B, (3, 3, 1))
    mm[:, :, :, 0] = translate(A, C, (3, 3, 1))
    mm[:, :, :, 1] = translate(A, D, (3, 3, 1))
    mm[:, :, :, 2] = translate(A, E, (3, 3, 1))

    move = [0 1 2 0 0; 0 0 0 1 2; 0 0 0 0 0]
    matrix = Float64[1 1 0 1 1 0 1 1 0]
    result = translate_optim(mm, matrix, move; g_abstol = 1e-14)

    @test minimum(result) < 1e-3
    #display(reshape(Optim.minimizer(result), 3, 3))
end


text:
module ImageTileArrays

using Requests
using ..DVIDArrays

export ImageTileArray

type ImageTileArray <: AbstractDVIDArray
    base        ::DVIDBase
end

function ImageTileArray( ip::String, port::Integer, node::AbstractString )
    ImageTileArray( DVIDBase(IPv4(ip), port, node) )
end

function ImageTileArray( ip::IPv4, port::Integer, node::AbstractString )
    ImageTileArray( DVIDBase(ip, port, node) )
end

function Base.ndims(A::ImageTileArray)
    3
end

function Base.eltype( A::ImageTileArray )
    UInt8
end

function Base.getindex(A::ImageTileArray, idxes::Union{Integer, UnitRange} ...)
    sz = map(x->length(x), idxes)
    # origin = map(first, idxes)
    # notethat the final coordinates are offset rather than julia 1-based coordinate
    link = "http://$(A.base.address):$(A.base.port)/api/node/$(A.base.node)"*
            "/grayscale/raw/0_1_2/$(sz[1])_$(sz[2])_$(sz[3])/"*
            "$(first(idxes[1]-1))_$(first(idxes[2]-1))_$(first(idxes[3]-1))"
    # @show link
    resp = Requests.get(link)
    reshape(resp.data, sz)
end

function Base.setindex!{T,N}(A::ImageTileArray, x::Array{T,N},
                        idxes::Union{Integer, UnitRange} ...)
    error("not working now!")
    sz = map(x->length(x), idxes)
    # origin = map(first, idxes)
    link = "http://$(A.base.address):$(A.base.port)/api/node/$(A.base.node)"*
            "/grayscale/raw/0_1_2/$(sz[1])_$(sz[2])_$(sz[3])/"*
            "$(first(idxes[1]))_$(first(idxes[2]))_$(first(idxes[3]))"
    resp = Requests.post()
end

end # end of module


text:
function neq(T::Real, mass::Real, g::Real, spin2::Int)
    x = mass / T
    η = spin2 % 2 == 0 ? 1 : -1
    # g * x^2 * T^3 / (2pi^2) * besselk(2, x)

    f(y) = y * sqrt(y^2 - 1) / (exp(y * x) - η)
    g * mass^3 / (2π^2) * quadgk(f, 1, Inf)[1]
end

function yeq(T::Real, mass::Real, g::Real, spin2::Int)
    x = mass / T
    η = spin2 % 2 == 0 ? 1 : -1
    # 45g * x^2 / (4π^4 * sm_heff(T)) * besselk(2, x)

    f(y) = y * sqrt(y^2 - 1) / (exp(y * x) - η)
    45g * x^3 / (4π^4 * sm_heff(T)) * quadgk(f, 1, Inf)[1]
end

function weq(T::Real, mass::Real, g::Real, spin2::Int)
    x = mass / T
    η = spin2 % 2 == 0 ? 1 : -1

    #log(45g * x^2 / (4π^4 * sm_heff(T)) * besselk(2, x))

    f(y) = y * sqrt(y^2 - 1) * exp(-x * (y - 1)) / (1 - η * exp(-y * x))
    -x + log(45g * x^3 / (4π^4 * sm_heff(T)) * quadgk(f, 1, Inf)[1])
end

using QuadGK

max_stars_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_issues_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_forks_repo_name: LoganAMorrison/AsymptoticBoltzmann


text:
function filter_gauss(image, sigma::Number=1.0, mask_order::Number=3)
  # Input description:
  #   image:      The image, wich is going to be blurred. The pixel-values must be an float matrix - range 0~1.
  #   sigma:      the sigma value for the mask creation, just as in the normal distribuition.
  #   mask_order: As the name says, it is basically the order of the mask, wich will be used to
  #               do the convolution trough the given image.
  #               mask_order must be a positive even number, bigger or equal to 3.

  tipo = typeof(image)

  try
    image = convert(Array{Float32}, image)
  catch er
    error(er," --- Problem converting the single band image into an array.")
  end

  d = size(image)

  img1 = zeros(d[1]-mask_order+1,d[2]-mask_order+1)

  # Creates the mask wich will be used
  factor = (1/(2*pi*(sigma^2)))
  a = floor(Int,mask_order/2)

  mask1 = ones(mask_order, mask_order)

  for x=-a:a
    for y=-a:a
      mask1[x+a+1,y+a+1] = factor*exp(-(x^2 + y^2)/(2*(sigma^2)))
    end
  end
  mask = mask1/sum(mask1)

  for i = 1:(d[1]-mask_order+1)
    for j = 1:(d[2]-mask_order+1)
      ai=i
      ia=i+mask_order-1
      aj=j
      ja=j+mask_order-1

      img1[i,j] = sum(mask.*image[ai:ia,aj:ja])
      if img1[i,j] > 1.0
          img1[i,j] = 1.0
      end
    end
  end

  img1 = convert(tipo, img1)

  return img1
end


text:
module BenchmarkTools

using JSON
using Base.Iterators

using Logging: @logmsg, LogLevel
using Statistics
using UUIDs: uuid4
using Printf

const BENCHMARKTOOLS_VERSION = v"0.7.0"

##############
# Parameters #
##############

include("parameters.jl")

export loadparams!

##############
# Trial Data #
##############

include("trials.jl")

export gctime,
       memory,
       allocs,
       params,
       ratio,
       judge,
       isinvariant,
       isregression,
       isimprovement,
       median,
       mean,
       rmskew!,
       rmskew,
       trim

##################
# Benchmark Data #
##################

include("groups.jl")

export BenchmarkGroup,
       invariants,
       regressions,
       improvements,
       @tagged,
       addgroup!,
       leaves

######################
# Execution Strategy #
######################

include("execution.jl")

export tune!,
       warmup,
       @ballocated,
       @benchmark,
       @benchmarkable,
       @belapsed,
       @btime

#################
# Serialization #
#################

include("serialization.jl")

end # module BenchmarkTools


text:
# This file was generated, do not modify it. # hide
using CSV, DataFrames, Latexify, Formatting # hide
D = CSV.read("voc-stats.tsv.gz", DataFrame)
println(latexify(D, latex=false, env=:mdtable, fmt = p -> format(p, commas=true ) * " "))# #hide


text:
#############################################################################
#Copyright (c) 2016 Charles Le Losq
#
#The MIT License (MIT)
#
#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.
#
#############################################################################

function gcvspl_julia(x::Array{Float64,1},y::Array{Float64,1},ese::Array{Float64,1},SmoothSpline::Float64;SplineOrder::Int32 = Int32(2),SplineMode::Int32 = Int32(3))
	"""

    c, wk, ier = gcvspline(x,y,ese,SplineSmooth,SplineOrder,SplineMode)

    INPUTS:

    x: Float64 Array, the independent variables
    y: Float64 Array, the observations (we assume here that you want to use this spline only on 1 dataset... see gcvspl.f if not)
    ese: Float64 Array, the errors on y
    SplineSmooth: Float64, the smoothing factor
    SplineOrder (M parameter in gcvspl.f): Int32, the half order of the required B-splines. default: splorder = 2 (cubic)
    SplineOrder = 1,2,3,4 correspond to linear, cubic, quintic, and heptic splines, respectively.
    SplineMode (MD parameter in gcvspl.f) is the Optimization mode switch:
	    default:   SplineMode = 3 (True predicted mean-squared error)
	               SplineMode = 1: Prior given value for p in VAL
	                         (VAL.ge.ZERO). This is the fastest
	                         use of GCVSPL, since no iteration
	                         is performed in p.
	               SplineMode = 2: Generalized cross validation.
	               SplineMode = 3: True predicted mean-squared error,
	                         with prior given variance in VAL.
	               SplineMode = 4: Prior given number of degrees of
	                         freedom in VAL (ZERO.le.VAL.le.N-M).
	               SplineMode  < 0: It is assumed that the contents of
	                         X, W, M, N, and WK have not been
	                         modified since the previous invoca-
	                         tion of GCVSPL. If MD < -1, WK(4)
	                         is used as an initial estimate for
	                         the smoothing parameter p.  At the
	                         first call to GCVSPL, MD must be > 0.
	               Other values for |MD|, and inappropriate values
	               for VAL will result in an error condition, or
	               cause a default value for VAL to be selected.
	               After return from MD.ne.1, the same number of
	               degrees of freedom can be obtained, for identical
	               weight factors and knot positions, by selecting
	               |MD|=1, and by copying the value of p from WK(4)
	               into VAL. In this way, no iterative optimization
	               is required when processing other data in Y.

    OUPUTS:

		c: the spline coefficients
		WK: work vector, see gcvspl.f
		IER: error parameter.
		IER = 0: Normal exit
		IER = 1: M.le.0 .or. N.lt.2*M
		IER = 2: Knot sequence is not strictly
		increasing, or some weight
		factor is not positive.
		IER = 3: Wrong mode parameter or value.

    SEE gcvspl.f FOR MORE INFORMATION

	"""
	if size(y,2)>1
		error("This function does not accept multiple dataset. Please provide only one dataset at a time.")
		#elseif size(x,1) != size(y,1) || size(ese,1) != size(y,1) || size(x,1) != size(ese,1)
		#error("This function only accepts x, y and ese arrays of the same lengths.")
	end

	WX = 1. ./((SmoothSpline.*ese).^2) # relative variance of observations
	WY = zeros([1])+1. # systematic errors... not used so put them to 1
	VAL = (SmoothSpline.*ese).^2

	# M = SplineOrder # No need to redefine but just this comment to keep record
	N = size(x,1) #same as length(y)
	K = 1 # number of y columns
	# MD = SplineMode # No need to redefine but just this comment to keep record
	NC = length(y)

	c = ones(N,NC)
	WK = ones(6*(N*SplineOrder+1)+N,1)
	IER=Int32[1]
	#pathname = abspath("../Dependencies/gcvspline/libgcvspl.so")
	#libgcvspl = Libdl.dlopen(abspath("../Dependencies/gcvspline/libgcvspl.so"))
	#const pathname = "/Users/charles/.julia/v0.4/Spectra/Dependencies/gcvspline/libgcvspl.so"
	#mylibvar = dlopen("mylib")
	ccall( (:gcvspl_,gcvspl), Void, (Ptr{Float64},Ptr{Float64},Ptr{Cint},Ptr{Float64},Ptr{Float64},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Float64},Ref{Float64},Ptr{Cint},Ref{Float64},Ref{Cint}),x,y,&N,WX,WY,&SplineOrder,&N,&K,&SplineMode,VAL,c,&NC,WK,IER)
	return c, WK, IER
end

function splderivative_julia(xfull::Array{Float64},xparse::Array{Float64},cparse::Array{Float64};SplineOrder::Int32 = Int32(2), L::Int32 = Int32(1), IDER::Int32 = Int32(0))
    """
    Wrapper to the SPLDER function of gcvspl.f, for interpolation purpose

	INPUTS:

	    xfull: Float64 Array, contains the entire x range where the spline has to be evaluated
	    xparse: Float64 Array, contains the x values of interpolation regions
	    WARNING!!! => xparse[0] <= xfull[0] <= xparse[n]
	    cparse: Float64 Array, is the evaluated spline coefficients returned by gcvspl for xparse

    OPTIONS:

        splineorder (integer): is the spline order,
        default: splineorder = 2 (cubic)
        L (integer): see gcvspl.f for details, default: L = 1
        IDER: the Derivative order required, with 0.le.IDER
        and IDER.le.2*M. If IDER.eq.0, the function
        value is returned; otherwise, the IDER-th
        derivative of the spline is returned.

    SEE gcvspl.f FOR MORE INFORMATION
    """
	# Note: In the following t, x and c in the fortran code are renamed  xfull, xparse and cparse for clarity
    N = Int32(size(xparse,1))
    q = zeros(2*SplineOrder,1)  # working array
    y_calc::Array{Float64} = zeros(size(xfull,1),1) # Output array
    #pathname = abspath("../Dependencies/gcvspline/libgcvspl.so")
	#pathname = "/Users/charles/.julia/v0.4/Spectra/Dependencies/gcvspline/libgcvspl.so"
	#libgcvspl = Libdl.dlopen(abspath("../Dependencies/gcvspline/libgcvspl.so"))
    # we loop other xfull to create the output values
	for i =1:size(y_calc,1)
	    y_calc[i] = ccall( (:splder_,gcvspl), Float64, (Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{Cint},Ptr{Float64}),&IDER, &SplineOrder, &N, &xfull[i], xparse, cparse, &L, q)
	end
	return y_calc
end


text:
using .DBInterface
struct Relational <: DatabaseKind end
function clean_sql(s::AbstractString) :: String
    return String(replace(strip(s), r" +"=>" "))
end

function unmarshal(mapper::DBMapper, ::Type{Relational}, dest::Type{ForeignKey{T}}, orig) where T <:Model
    id_field_name = idfield(mapper, T)
    params = Dict{Symbol, Any}(id_field_name=>orig)
    return ForeignKey{T}(data=T(;params...), loaded=false)
end

database_column_type(dbtype, d::Union{Symbol, Type}) = raise("Not implemented")

primary_key_type(dbtype::DataType, x) = x

function create_table_field(mapper::DBMapper, field::Field, table::Table, dbtype::DataType)
    field_name = "$(field.name)"
    if field.type <: ForeignKey
        referenced_table = mapper.tables[element_type(field.type)]
        referenced_field = referenced_table.primary_key.field[1]
        db_field_type = string(database_column_type(dbtype, element_type(referenced_field.type)))
    elseif field.type <: DBId
        pk_type = primary_key_type(dbtype, element_type(field.type))
        db_type = database_column_type(dbtype, pk_type)
        db_field_type = string(db_type)
    else
        db_field_type  = string(database_column_type(dbtype, element_type(field.type)))
    end
    primary_key = isprimarykey(field, table) ? "PRIMARY KEY" : ""
    nullable = field.nullable ? "" : "NOT NULL"
    return strip("$field_name $db_field_type $primary_key $nullable")
end

on_event_action(t::OnEventAction) = throw("Not implemented")
on_event_action(t::Cascade) = "CASCADE"
on_event_action(t::DoNothing) = "NO ACTION"
on_event_action(t::SetNull) = "SET NULL"
on_event_action(t::Restrict) = "RESTRICT"

function create_reference_field(r::Relation)
    return strip("FOREIGN KEY($(r.local_field)) " *
            "REFERENCES $(r.referenced_table)($(r.referenced_field)) " *
            "ON DELETE $(on_event_action(r.on_delete)) " *
            "ON UPDATE $(on_event_action(r.on_update))")

end

function create_table_query(mapper::DBMapper, T::Type{<:Model}; if_not_exists::Bool=true) :: String
    if mapper.dirty == true
        analyze_relations(mapper)
    end
    table = mapper.tables[T]
    create_table_fields = []
    for field in fieldlist(table)
        push!(create_table_fields,
              create_table_field(mapper, field, table, mapper.pool.dbtype))
    end

    foreign_keys = [create_reference_field(r) for r in values(table.relations)]
    append!(create_table_fields, foreign_keys)

    create_table_fields = join(create_table_fields, ", ")
    if_not_exists_str = if_not_exists ? "IF NOT EXISTS" : ""
    return clean_sql("""CREATE TABLE $if_not_exists_str $(table.name) ($create_table_fields)""")
end

function create_table(mapper::DBMapper, dbtype::Type{Relational}, T::Type{<:Model}; if_not_exists::Bool=true)
    sql = create_table_query(mapper, T, if_not_exists=if_not_exists)
    @info sql
    conn = get_connection(mapper.pool)
    result = DBInterface.execute(conn, sql)
    release_connection(mapper.pool, conn)
    return result
end

insert_query(table::Table, column_names::Array, dbtype) = insert_query(table, column_names)

function insert_query(table::Table, column_names::Array)
    values_placeholder = join(repeat(['?'], length(column_names)), ",")
    column_names = join(column_names, ",")
    sql = """
INSERT INTO $(table.name) ($column_names)
VALUES ($values_placeholder)
    """
    return clean_sql(sql)
end

"""
    function insert!(mapper::DBMapper, dbtype::Type{Relational}, elem::T) where T

Insert the element in the database. Update the id of the element
"""
function insert!(mapper::DBMapper, dbtype::Type{Relational}, elem::T) where T<:Model
    (column_names, values) = struct_field_values(mapper, elem)
    table = mapper.tables[T]
    sql = insert_query(table, column_names, mapper.pool.dbtype)
    @info sql
    conn = get_connection(mapper.pool)
    stmt = DBInterface.prepare(conn, sql)
    result = DBInterface.execute(stmt, values)
    id = DBInterface.lastrowid(result)
    release_connection(mapper.pool, conn)
    if table.primary_key.has_auto_value
        setid!(elem, mapper, id)
    end
    return elem
end

update_query(table::Table, column_names::Array, dbtype) = update_query(table, column_names)

function update_query(table::Table, column_names::Array)
    column_names = join(map(cname->"$cname=?", column_names), ",")
    sql = """
    UPDATE $(table.name)
    SET $column_names
    WHERE
    $(idfield(table)) = ?
    """
    return clean_sql(sql)
end

"""
    function insert!(mapper::DBMapper, dbtype::Type{Relational}, elem::T) where T

Insert the element in the database. Update the id of the element
"""
function update!(mapper::DBMapper, dbtype::Type{Relational}, elem::T; fields::Array{Symbol}=Symbol[]) where T<:Model
    id = getid(elem, mapper)
    if id === nothing
        throw("Id not present, cannot update")
    end
    (column_names, values) = struct_field_values(mapper, elem; ignore_primary_key=true, fields=fields)
    table = mapper.tables[T]
    sql = update_query(table, column_names, mapper.pool.dbtype)
    push!(values, id)
    @info sql
    conn = get_connection(mapper.pool)
    stmt = DBInterface.prepare(conn, sql)
    result = DBInterface.execute(stmt, values)
    release_connection(mapper.pool, conn)
    return elem
end

"""
    function totuple(mapper::DBMapper, table::Table, dbtype::DataType, db_results) :: Array{Array{Pair}}

Return an array of array of tuples (field=>value). The values are converted to the julia types
"""
function totuple(mapper::DBMapper, table::Table, dbtype::DataType, db_results) :: Array{Array{Pair}}
    results = []
    for row in db_results
        r = []
        db_data = Dict(field=>getindex(row, field)
                      for field in propertynames(row))
        for field in fieldlist(table)
            if haskey(db_data, field.name)
                push!(r, field.struct_field=>unmarshal(mapper, dbtype, field.type, db_data[field.name]))
            end
        end
        push!(results, r)
    end
    return results
end

function totuple(results)
    [(;(prop=>getindex(row, prop) for prop in propertynames(row))...) for row in results]
end

function select_query(mapper::DBMapper, T::Type{<:Model}; select_one=true, fields::Array{Symbol}=Symbol[], kwargs...)
    dbtype = mapper.pool.dbtype
    table = mapper.tables[T]
    if isempty(fields)
        fields = column_names(mapper, T)
    end
    cnames = join(fields, ", ")
    conditions = join(["$field=?" for (field, value) in kwargs], " AND ")
    sql = """
    SELECT $cnames
    FROM $(table.name)
    """
    if !isempty(conditions)
        sql *= """
        WHERE $(conditions)"""
    end
    if select_one
        sql *= " LIMIT 1"
    end
    return clean_sql(sql)
end

function select(mapper::DBMapper, T::Type{<:Model}; select_one=true, fields::Array{Symbol}=Symbol[], kwargs...)
    dbtype = mapper.pool.dbtype
    table = mapper.tables[T]
    sql = select_query(mapper, T; select_one=select_one, fields=fields, kwargs...)
    values = [v[2] for v in kwargs]
    @info sql
    conn = get_connection(mapper.pool)
    stmt = DBInterface.prepare(conn, sql)
    result = totuple(mapper, table, dbtype, DBInterface.execute(stmt, values))
    release_connection(mapper.pool, conn)
    return result
end

function select_one(mapper::DBMapper, ::Type{Relational}, T::Type{<:Model}; kwargs...)
    result = select(mapper, T; select_one=true, kwargs...)
    if isempty(result)
        return nothing
    else
        return T(;result[1]...)
    end
end

function select_all(mapper::DBMapper, ::Type{Relational}, T::Type{<:Model}; fields::Array{Symbol} =Symbol[], kwargs...) :: Array{T}
    result = select(mapper, T; select_one=false, fields=fields, kwargs...)
    if isempty(result)
        return T[]
    else
        return [T(;r...) for r in result]
    end
end

clean_table_query(table::Table,  dbtype) = clean_table_query(table)

function clean_table_query(table::Table)
    return clean_sql("TRUNCATE TABLE $(table.name)")
end

function clean_table!(mapper::DBMapper, dbtype::Type{Relational}, T::Type{<:Model})
    table = mapper.tables[T]
    sql = clean_table_query(table, mapper.pool.dbtype)
    conn = get_connection(mapper.pool)
    result = DBInterface.execute(conn, sql)
    release_connection(mapper.pool, conn)
end

drop_table_query(table::Table,  dbtype) = drop_table_query(table)

function drop_table_query(table::Table)
    return clean_sql("DROP TABLE $(table.name)")
end

function drop_table!(mapper::DBMapper, dbtype::Type{Relational}, T::Type{<:Model})
    table = mapper.tables[T]
    sql = drop_table_query(table, mapper.pool.dbtype)
    @info sql
    conn = get_connection(mapper.pool)
    result = DBInterface.execute(conn, sql)
    release_connection(mapper.pool, conn)
end

function exists(mapper::DBMapper, ::Type{Relational}, T::Type{<:Model}; kwargs...) :: Bool
    table = mapper.tables[T]
    conditions = join(["$field = ?" for (field, value) in kwargs], " AND ")
    values = [v[2] for v in kwargs]
    sql = """
    SELECT COUNT(1) as count
    FROM $(table.name)
    WHERE $conditions
    """
    @info sql
    conn = get_connection(mapper.pool)
    stmt = DBInterface.prepare(conn, sql)
    result = DBInterface.execute(stmt, [value for (f, value) in kwargs])
    release_connection(mapper.pool, conn)
    result = result |> totuple
    return result[1][:count] > 0
end

function Base.delete!(mapper::DBMapper, ::Type{Relational}, T::Type{<:Model}; kwargs...)
    dbtype = mapper.pool.dbtype
    table = mapper.tables[T]
    cnames = join(column_names(mapper, T), ", ")
    conditions = join(["$field=?" for (field, value) in kwargs], " AND ")
    values = [v[2] for v in kwargs]
    sql = clean_sql("""
    DELETE FROM $(table.name)
    WHERE $(conditions)
    """)
    @info sql
    conn = get_connection(mapper.pool)
    stmt = DBInterface.prepare(conn, sql)
    DBInterface.execute(stmt, values)
    release_connection(mapper.pool, conn)
end

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
export TparamsKTR

mutable struct TparamsKTR{T} <: Tparams{T}  # specific parameters for this Krylov variant
    ζ :: T
    shifts :: Array{T, 1}
    nshifts :: Int
    function TparamsKTR{T}(; ζin ::T = 0.5) where T
        ζ = ζin
        shifts = [0.0; 10.0.^(collect(-15.0:1.0:15.0))]
        nshifts = size(shifts, 1)
        return new{T}(ζ, shifts, nshifts)
    end
end


text:
export pca, qpca, nnmf, rpca, kmeans

# principal components analysis
# minimize ||A - XY||^2
function pca(A::AbstractArray, k::Int; kwargs...)
	loss = QuadLoss()
	r = ZeroReg()
	return GLRM(A,loss,r,r,k; kwargs...)
end

# quadratically regularized principal components analysis
# minimize ||A - XY||^2 + scale*||X||^2 + scale*||Y||^2
function qpca(A::AbstractArray, k::Int; scale=1.0::Float64, kwargs...)
	loss = QuadLoss()
	r = QuadReg(scale)
	return GLRM(A,loss,r,r,k; kwargs...)
end

# nonnegative matrix factorization
# minimize_{X>=0, Y>=0} ||A - XY||^2
function nnmf(A::AbstractArray, k::Int; kwargs...)
	loss = QuadLoss()
	r = NonNegConstraint()
	GLRM(A,loss,r,r,k; kwargs...)
end

# k-means
# minimize_{columns of X are unit vectors} ||A - XY||^2
function kmeans(A::AbstractArray, k::Int; kwargs...)
	loss = QuadLoss()
	ry = ZeroReg()
	rx = UnitOneSparseConstraint()
	return GLRM(A,loss,rx,ry,k; kwargs...)
end

# robust PCA
# minimize HuberLoss(A - XY) + scale*||X||^2 + scale*||Y||^2
function rpca(A::AbstractArray, k::Int; scale=1.0::Float64, kwargs...)
	loss = HuberLoss()
	r = QuadReg(scale)
	return GLRM(A,loss,r,r,k; kwargs...)
end


text:
using ONNX, Flux, ProtoBuf
using DataFlow: Call, vertex, syntax, constant
using Test
using Base:run
# test taken from : https://github.com/onnx/onnx/tree/master/onnx/backend/test/data
# clone onnx here if onnx dir does not exist

if !("onnx" in readdir())
    # clone the package here
    println("Downloading test data....")
    Base.run(`git clone https://github.com/onnx/onnx.git`)
end

ONNX_PATH = "./onnx"

ONNX_TEST_PATH = "$ONNX_PATH/onnx/backend/test/data/node"

function read_input(folder_name)
    ar = Array{Any, 1}()
    for ele in readdir(folder_name*"/test_data_set_0")
        push!(ar, Float32.(readproto(open(folder_name*"/test_data_set_0/"*ele), ONNX.Proto.TensorProto()) |> ONNX.get_array))
    end
    return ar[1:end-1]
end

function read_output(folder_name)
    ar = Array{Any, 1}()
    for ele in readdir(folder_name*"/test_data_set_0")
        push!(ar, Float32.(readproto(open(folder_name*"/test_data_set_0/"*ele), ONNX.Proto.TensorProto()) |> ONNX.get_array))
    end
    return ar[end]
end

function read_model(folder_name)
    for ele in readdir(folder_name)
        if ele=="model.onnx"
            return readproto(open(folder_name*"/model.onnx"), ONNX.Proto.ModelProto())
        end
    end
end

function get_optype(a::ONNX.Proto.ModelProto)
    g = ONNX.convert(a.graph)
    return g.node[1].op_type
end

function get_dict(a::ONNX.Proto.ModelProto)
    g = ONNX.convert(a.graph)
    return g.node[1].attribute
end

function main_test(filename,op_expected, ip...)
        if Symbol(get_optype(read_model(filename))) == :Constant
            @test get_dict(read_model(filename))[:value] |> ONNX.get_array == op_expected

        elseif Symbol(get_optype(read_model(filename))) == :Conv

            temp = ONNX.ops[Symbol(get_optype(read_model(filename)))](get_dict(read_model(filename)),
                                                                                Symbol("ip[1]"), Symbol("ip[2]")) |> syntax

            touch("temp_conv.jl")
            open("temp_conv.jl","w") do file
                str1 = "flipkernel(x) = x[end:-1:1, end:-1:1, :, :] \n"                     # Remove when Flux directly supports it
                write(file, str1*string(temp))
            end
            model = include("temp_conv.jl")
            rm("temp_conv.jl")
            @test model == op_expected
        elseif Symbol(get_optype(read_model(filename))) == :MaxPool
            temp = ONNX.ops[Symbol(get_optype(read_model(filename)))](get_dict(read_model(filename)),
                                                                                Symbol("ip[1]")) |> syntax
            touch("temp_maxpool.jl")
            open("temp_maxpool.jl","w") do file
                write(file, string(temp))
            end

            model = include("temp_maxpool.jl")
            rm("temp_maxpool.jl")
            @test model == op_expected
        elseif Symbol(get_optype(read_model(filename))) == :AveragePool
            temp = ONNX.ops[Symbol(get_optype(read_model(filename)))](get_dict(read_model(filename)),
                                                                                Symbol("ip[1]")) |> syntax
            touch("temp_averagepool.jl")
            open("temp_averagepool.jl","w") do file
                write(file, string(temp))
            end

            model = include("temp_averagepool.jl")
            rm("temp_averagepool.jl")
            @test model ≈ op_expected atol=0.001
        elseif Symbol(get_optype(read_model(filename))) in [:GlobalAveragePool, :GlobalMaxPool]
            temp = ONNX.ops[Symbol(get_optype(read_model(filename)))](get_dict(read_model(filename)),
                                                                                Symbol("ip[1]")) |> syntax
            touch("temp_averagepool.jl")
            open("temp_averagepool.jl","w") do file
                write(file, "using Statistics \n" * string(temp))
            end
            model = include("temp_averagepool.jl")
            rm("temp_averagepool.jl")
            @test model ≈ op_expected atol=0.001
        elseif Symbol(get_optype(read_model(filename))) in [:Expand, :Concat]
            temp = ONNX.ops[Symbol(get_optype(read_model(filename)))](get_dict(read_model(filename)),
                                                                                Symbol("ip[1]"), Symbol("ip[2]")) |> syntax
            touch("temp_expand.jl")
            open("temp_expand.jl","w") do file
                write(file, string(temp))
            end

            model = include("temp_expand.jl")
            rm("temp_expand.jl")
            @test model ≈ op_expected atol=0.001
        else
                @test ONNX.ops[Symbol(get_optype(read_model(filename)))](get_dict(read_model(filename)),
                                 ip...) |> syntax |> eval ≈ op_expected atol=0.001
        end
end