text:
# include("Load_LTSpice_Net.jl")
# include("RunAnalysis.jl")
# include("SPICE2Matrix.jl")
# include("ToroidOptimizer.jl")

"""
This function takes in  the following parameters:

    LDrive, drive coil inductance (H)

    RDrive, Drive coil resistance (Ω)

    TargetZ, Amplifier/target impedance to match to  (Ω)

    DriveFreq; drive frequency (Hz)
The following are keyword arguments and should be entered after the first four (required) as;
DesignDriveFilter(LDrive,RDrive,Amp_Z,DriveFreq; [kwarg] = value). See below for an example

    CDrive = 1e6, Any series capacitance with drive coil

    NumDriveElements = 1, The number of series copies of the L,C,R for example if there are 6 pancakes of LDrive, RDrive, Cdrive, then this would be 6

    WireDiam = 2e-3, Wire diameter for toroids

    WireFillFac = 0.75, If litz wire, this is the ratio of copper to air

    PlotFTs = true, Plotting be default is true
    VSrc = 2, The source voltage for plotting
    DetermineFreq = false, IF you have a fixed capacitor (CDrive has been set) you may want to search for an ideal drive freq. in which case, set this to true.
    AddNotchFreq = nothing, This can be a single value of a notch frequency, or an array of length 2, the two values are the notch frequencies

    FilterZ = TargetZ, The characteristic impedance of the filtering section
    RDampVal = 100, [Ω] The damping resistor added in around the series L-C sections
    PerturbTxReactance = nothing, an optional element to add some purturbation to the reactance of the drive coil, which effectively enables the user to tune where the operating point ends up on the transfer function (V in → I drive)

    Example (Copy/paste):
    RDrive = 400e-3
    LDrive = 400e-6
    Amp_Z = 4
    DriveFreq = 25e3
    AddNotchFreqList = [DriveFreq*2, DriveFreq*3]

    DesignDriveFilter(LDrive,RDrive,Amp_Z,DriveFreq; AddNotchFreq = AddNotchFreqList)

"""
function DesignDriveFilter(
    LDrive,
    RDrive,
    TargetZ,
    DriveFreq;
    CDrive = 2,
    NumDriveElements = 1,
    WireDiam = 2e-3,
    WireFillFac = 0.75,
    PlotSrcR = TargetZ,
    PlotFTs = true,
    VSrc = 2,
    DetermineFreq = false,
    AddNotchFreq = nothing,
    FilterZ = TargetZ,
    RDampVal = 100,
    PerturbTxReactance = nothing,
    AxesArray = nothing,
    WriteFileName=nothing
)

    ZeroVal = 1e-9 #There are issues with zero-valued components
    ωDr = 2 * π * DriveFreq
    ZDrive =
        RDrive * NumDriveElements +
        im * ωDr * LDrive * NumDriveElements +
        NumDriveElements ./ (im * 2 * pi * DriveFreq * CDrive) #Equivalent complex impednace of the load

    Reactance_Load =
        ωDr * LDrive * NumDriveElements - NumDriveElements ./ (ωDr * CDrive)
    Reactance_Load = round(Reactance_Load; sigdigits = 3)
    display(println("The reactance of the load is: $(round(Reactance_Load;sigdigits=3)) Ω"))
    if ~DetermineFreq
        if (Reactance_Load > 0)

            SerCap,CParAct = ImpMatch_LLoad(TargetZ, ZDrive, Reactance_Load,ωDr)

            LTee_2 = ZeroVal
            LTee_2_ESR = ZeroVal
            LTee_1 = ZeroVal
            LTee_1_ESR = ZeroVal

        elseif (Reactance_Load < 0)

            LTee_2,LTee_2_ESR,CParAct = ImpMatch_CLoad(TargetZ, ZDrive, Reactance_Load,ωDr)

            LTee_1 = ZeroVal
            LTee_1_ESR = ZeroVal
            SerCap = 2

        elseif (Reactance_Load == 0)
            matchRatio = real(TargetZ) / real(ZDrive)
            Q = sqrt(matchRatio - 1)
            Xs = Q * real(ZDrive)
            LSer2 = Xs / (ωDr)
            LTee_2 = LSer2
            LTee_1 = ZeroVal
            LTee_1_ESR = ZeroVal
            CParAct = findResPair((1 + Q^(-2)) * LSer2, DriveFreq)
            LTee_2_Geom =
                ToroidOptimizer(WireDiam, LTee_2; CuFillFactor = WireFillFac)
            LTee_2_ESR = LTee_2_Geom.DCore.Resistance
            SerCap = 2
        end
    else

        matchRatio = real(TargetZ) / real(ZDrive)
        Q = √(matchRatio - 1)
        Xs = Q * real(ZDrive)

        DriveFreq = findFilterFreq(Q,real(ZDrive), LDrive * NumDriveElements, CDrive / NumDriveElements)
        ωDr = 2*π*DriveFreq
        # println("Determined Q to be $(round(Q))")
        # println("Determined drive freq. to be $(round(DriveFreq)) Hz")
        ZSer = 1*im * ωDr*LDrive * NumDriveElements -  1*im /(ωDr* CDrive / NumDriveElements) +RDrive* NumDriveElements
        YSer = 1/ZSer

        CParAct = abs.(imag(YSer))/ωDr
        LTee_2 = ZeroVal
        LTee_2_ESR = ZeroVal
        LTee_1 = ZeroVal
        LTee_1_ESR = ZeroVal
        SerCap = 1e6
    end

    (LFilt, CFilt) = Butterworth_2(FilterZ, DriveFreq)
    CFilt2 = findResPair(LFilt, DriveFreq)
    LFilt2 = findResPair(CFilt, DriveFreq)

    LFilt1_Geom =
            ToroidOptimizer(WireDiam, LFilt; CuFillFactor = WireFillFac)
    LFilt1_ESR = LFilt1_Geom.DCore.Resistance

    LFilt2_Geom =
            ToroidOptimizer(WireDiam, LFilt2; CuFillFactor = WireFillFac)
    LFilt2_ESR = LFilt2_Geom.DCore.Resistance

    CalcC(r,xl) = 1/(ωDr*(r^2/xl+xl))
    C_Par_2 = abs(CalcC(TargetZ, abs.(im * ωDr * LDrive * NumDriveElements + NumDriveElements ./ (im * ωDr * CDrive)+1/(im*ωDr*SerCap))))
    X_Par_2 =  1 ./ (im * ωDr * C_Par_2)
    println(C_Par_2)
    if ~(AddNotchFreq===nothing)
        if length(AddNotchFreq)==1
            LNotch, CNotch = makeNotchSection(AddNotchFreq[1], DriveFreq, X_Par_2)
        else
            LNotch, CNotch = makeNotchSection(AddNotchFreq[1], DriveFreq, X_Par_2*2)
        end
        LNotch_Geom =
                ToroidOptimizer(WireDiam, LNotch; CuFillFactor = WireFillFac)
        LNotch_ESR = LNotch_Geom.DCore.Resistance

        println("Adding notch")
        if length(AddNotchFreq)==2
            LNotch2, CNotch2 = makeNotchSection(AddNotchFreq[2], DriveFreq, X_Par_2*2)
            LNotch2_Geom =
                    ToroidOptimizer(WireDiam, LNotch2; CuFillFactor = WireFillFac)
            LNotch2_ESR = LNotch2_Geom.DCore.Resistance

            println("Adding second notch")
        else
            LNotch2 = 10
            CNotch2 = 1e-12
            LNotch2_ESR = 1e-3

        end

    else
        LNotch = 1e-9
        CNotch = C_Par_2

        LNotch2 = 10 #Set default values is no notch is added
        CNotch2 = 1e-12
        LNotch2_ESR = 1e-3
    end

    CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList,AxesArray = CircModel_SPICE(DriveFreq, VSrc,
    1e-3,
    RDrive,
    NumDriveElements,
    LDrive,
    CDrive,
    SerCap,
    CParAct,
    LFilt,
    LFilt2_ESR,
    LFilt1_ESR,
    CFilt2,
    CFilt,
    LFilt2,
    LNotch,
    LNotch_ESR,
    CNotch,
    LNotch2,
    LNotch2_ESR,
    CNotch2;
    PlotOn = PlotFTs,
    RDampVal = RDampVal,
    AxesArray = AxesArray)

    if ~(WriteFileName===nothing)
    writeFilterSPICE(WriteFileName,
    RDrive,
    NumDriveElements,
    LDrive,
    CDrive,
    SerCap,
    CParAct,
    LFilt,
    LFilt2_ESR,
    LFilt1_ESR,
    LFilt2,
    CFilt2,
    CFilt,
    LNotch,
    LNotch_ESR,
    CNotch,
    LNotch2,
    LNotch2_ESR,
    CNotch2,
    RDampVal,
    NumDriveElements*LDrive,
    Par(CDrive/NumDriveElements,SerCap))
    end
    DetermineComponentsTempCoeffs(SPICE_DF,InputList,1,DriveFreq,"LDrive")

    return DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList,AxesArray

end

function CircModel_SPICE(DriveFreq, VSrc,
    PlotSrcR,
    RDrive,
    NumDriveElements,
    LDrive,
    CDrive,
    SerCap,
    CParAct,
    LFilt,
    LFilt2_ESR,
    LFilt1_ESR,
    CFilt2,
    CFilt,
    LFilt2,
    LNotch,
    LNotch_ESR,
    CNotch,
    LNotch2,
    LNotch2_ESR,
    CNotch2;
    PlotOn = true,
    FreqList = 1000:10:100e3,
    ArchetypeNetFileName = nothing,
    RDampVal = nothing,
    AxesArray=nothing)

    if ArchetypeNetFileName === nothing

        ArchetypeNetFileName = joinpath(dirname(pathof(DesignDriveCoilFilt)),"Filter_Archetype_Damped_3.net")

    end
    SPICE_DF,NodeList,InputList,NumVSources = SPICE2Matrix(ArchetypeNetFileName)
    UpdateElementVal!(SPICE_DF,"RSrc",PlotSrcR)

    UpdateElementVal!(SPICE_DF,"LNotch",LNotch)
    UpdateElementESR!(SPICE_DF,"LNotch",LNotch_ESR)
    UpdateElementVal!(SPICE_DF,"CNotch",CNotch)
    UpdateElementVal!(SPICE_DF,"LNotch1",LNotch2)
    UpdateElementESR!(SPICE_DF,"LNotch1",LNotch2_ESR)
    UpdateElementVal!(SPICE_DF,"CNotch1",CNotch2)

    UpdateElementVal!(SPICE_DF,"LFilt1",LFilt)
    UpdateElementVal!(SPICE_DF,"LFilt2",LFilt)
    UpdateElementVal!(SPICE_DF,"CFilt1",CFilt2)
    UpdateElementVal!(SPICE_DF,"CFilt2",CFilt2)
    UpdateElementESR!(SPICE_DF,"LFilt1",LFilt1_ESR)

    UpdateElementVal!(SPICE_DF,"CFilt3",CFilt)
    UpdateElementVal!(SPICE_DF,"LFilt3",LFilt2)
    UpdateElementESR!(SPICE_DF,"LFilt3",LFilt2_ESR)

    UpdateElementVal!(SPICE_DF,"RDrive",RDrive)
    UpdateElementVal!(SPICE_DF,"LDrive",NumDriveElements*LDrive)
    UpdateElementVal!(SPICE_DF,"CDrive",CDrive/NumDriveElements)

    UpdateElementVal!(SPICE_DF,"CSer",SerCap)
    UpdateElementVal!(SPICE_DF,"CPar",CParAct)

    UpdateElementVal!(SPICE_DF,"LDrivePair",NumDriveElements*LDrive)
    UpdateElementESR!(SPICE_DF,"LDrivePair",LFilt1_ESR)
    UpdateElementVal!(SPICE_DF,"CDrivePair",Par(CDrive/NumDriveElements,SerCap))

    if ~(RDampVal===nothing)

        UpdateElementVal!(SPICE_DF,"RDamp2",RDampVal)
        UpdateElementVal!(SPICE_DF,"RDamp1",RDampVal)

    end

    inputs = zeros(length(InputList))
    inputs[end-(NumVSources-1):end] .= 1

    ResultNodeNames = vcat("V(".*NodeList.*")", "I(".*(InputList[end-(NumVSources-1):end]).*")")

    Results = [(inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*FreqList[i],InputList))*inputs) for i in 1:length(FreqList)]
    Results = hcat(Results...)
    ResDict = Dict(ResultNodeNames[1] => Results[1,:])
    for i in 2:length(ResultNodeNames)

         merge!(ResDict,Dict(ResultNodeNames[i] => Results[i,:]))

    end

    CurrentVec = VSrc.* plotACElCurrent(SPICE_DF,FreqList,Results,"LDrive")

    if PlotOn
        if (AxesArray===nothing)
            AxMain = PyPlot.axes([0.13, 0.15, 0.85, 0.85])

            AxesNewArray = [AxMain, AxMain] #Initialize with length 2
        else
            PyPlot.axes(AxesArray[1])
        end
        pygui(true)
        plot(FreqList,(20*log10.(abs.(CurrentVec[:]))))
        xlabel("Frequency, Hz")
        ylabel("Current in LDrive [dBAmps per Volt input]")

        StartIndex = findfirst(x-> x>(DriveFreq-1e3),FreqList)
        StopIndex = findfirst(x-> x>(DriveFreq+1e3),FreqList)
        if (AxesArray===nothing)
            AxInset = PyPlot.axes([0.68, 0.65, 0.25, 0.25]) #[left, bottom, width, height]
            AxesNewArray[2] = AxInset
        else
            PyPlot.axes(AxesArray[2])
        end

        plot(FreqList[StartIndex:StopIndex],(20*log10.(abs.(CurrentVec[StartIndex:StopIndex]))))
        plotDriveFreqDot(FreqList,CurrentVec,DriveFreq)

    end
    if (AxesArray===nothing)
        AxesArray = AxesNewArray
    end
    getElementCurrents(SPICE_DF,Results,DriveFreq)

    return CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList,AxesArray
end

function ImpMatch_LLoad(TargetZ, ZDrive, Reactance_Load,ωDr)
    println("Load is inductive")
    EquivSerL = Reactance_Load / ωDr
    println(
        "Load appears to be a inductor with a value of $(round(EquivSerL*1e6;sigdigits=3))μH and a ESR of $(real(ZDrive))Ω ",
    )

    matchRatio = real(TargetZ) / real(ZDrive)
    Q = sqrt(matchRatio - 1)
    Xs = Q * real(ZDrive) #Target reactance

    X_SerCap = Reactance_Load - Xs #Reactance of series capacitor
    SerCap = 1 / (ωDr * X_SerCap) #Series Capacitor value
    CParAct = Q / (ωDr * TargetZ)

    return SerCap,CParAct

end

function ImpMatch_CLoad(TargetZ, ZDrive, Reactance_Load,ωDr; WireDiam = 2e-3,
    WireFillFac = 0.75)
    println("Load is capacitive")
        EquivSerC = abs.(1 / (Reactance_Load * ωDr))
        println(
            "Load appears to be a capacitor with a value of $(round(EquivSerC*1e6;sigdigits=3))μF and a ESR of $(real(ZDrive))Ω ",
        )
        LSer = abs.(Reactance_Load) / (ωDr)
        matchRatio = real(TargetZ) / real(ZDrive)
        Q = sqrt(matchRatio - 1)
        Xs = Q * real(ZDrive)
        LSer2 = Xs / (ωDr)
        LTee_2 = LSer2 + LSer

        CParAct = findResPair((1 + Q^(-2)) * LSer2, ωDr*2*π)
        LTee_2_Geom =
            ToroidOptimizer(WireDiam, LTee_2; CuFillFactor = WireFillFac)
        LTee_2_ESR = LTee_2_Geom.DCore.Resistance
        println("Q = $Q")

        return LTee_2,LTee_2_ESR,CParAct

end

"""
This function determines the ideal frequency of a resonant RLC such that it can be matched to a different impedance. It assumes there is already an L, C, and R, so at some frequency it the real part of the admittance will be the matched admittance. This admittance is made fully real with a shunt reactive element

"""
function findFilterFreq(Q,RLoad,L,C)
        f = ((Q*RLoad + √(Q^2*RLoad^2+4*L/C)) / (4*π*L) , (Q*RLoad - √(Q^2*RLoad^2+4*L/C)) / (4*π*L) )

        return maximum(f)
end

"""
This function makes a filter section that notches out a frequency (NotchFreq), but is apparently real at another (DriveFreq)
    the inputs are:
        * NotchFreq: the freq that will be shunted out
        * DriveFreq: The freq that should remain real
        * RMatch: The resistance the drive coil is tuned to
        * LVal is a kwarg because either the inductor or capacitor value should be pre-defined
The outputs are:
    LVal(Henries)
    CVal(Farads)
    LTune(Henries)
"""
function makeNotchSection(NotchFreq, DriveFreq, XEquiv)
    println(XEquiv)
    LVal,CVal = findEquivLC(XEquiv,DriveFreq,NotchFreq)
    return abs(LVal),abs(CVal)
end

"""

Calculates parallel impedances.
If given two inputs, it assumes both are inputs
If given an array, it takes the parallel combo of it all

Eg.
ZList = [1 + 1*im, 1]
Par(ZList[1], ZList[2])
Par(ZList)
"""
function Par(Z1, Z2)
    return 1 / (1 / Z1 + 1 / Z2)
end
function Par(Z::Array)
    Y = 1 ./Z
    YTotal = sum(Y)
    Zeff = 1/YTotal
    return Zeff
end

"""
Parallel component lists
Input is Array{Tuple{Number,String},1}
e.g.
Freq = 25e3
CompList = [(3, "R"),(30e-6, "L"),(200e-9,"C")]
Par(CompList,Freq)

"""
function Par(C::Array{Tuple{Real,String},1},freq)

    Z = Complex.(zeros(length(C)))
    for i in 1:length(C)
        if (lowercase(C[i][2])=="r") | (lowercase(C[i][2])=="resistor")
            Z[i] = C[i][1]
        elseif (lowercase(C[i][2])=="l") | (lowercase(C[i][2])=="inductor" )
            Z[i] = Z_Ind(C[i][1],freq)
        elseif (lowercase(C[i][2])=="c") | (lowercase(C[i][2])=="capacitor" )
            Z[i] = Z_Cap(C[i][1],freq)
        elseif lowercase(C[i][2])=="z"
            Z[i] = C[i][1]
        else
            error("Unknown component")
        end

    end

    Y = 1 ./Z
    YTotal = sum(Y)
    Zeff = 1/YTotal
    return Zeff
end

"""
Converts from a polar coordinate (Mag, phase) to cartesian
if you enter a tuple or array it assumes first is mag.

"""
function Polar2Cart(Mag, Phase)
    Real = real.(Mag * exp(-im * Phase / 360 * 2 * pi))
    Imag = imag.(Mag * exp(-im * Phase / 360 * 2 * pi))
    return (Real, Imag)
end

function Polar2Cart(MagPhaseTup::Union{Tuple, Array})
    Polar2Cart(MagPhaseTup[1], MagPhaseTup[2])
end

function Z_Cap(C, f)
    return (1 ./ (im * 2 * pi .* f .* C))
end
function Z_Ind(L, f)
    return im * 2 * pi .* f .* L
end

"""
This function takes in some reactance element (E) and a frequency, and finds the pair

"""
function findResPair(E, f)
    ##Finds the resonant pair for a reactive element, E
    ## ωL = 1/(ωC) → L = 1/(ω^2C),and vice versa
    return 1 ./ ((2 * pi * f) .^ 2 * E)
end

function Butterworth_2(Z, f)
    #Sqrt(L/C) = Z ; f = 1/(2π*sqrt(LC))
    #L = Z^2*C
    #f = 1/(2π*Z C) → C = 1/(2*pi*f*Z)
    #Butterworth 2nd order = [LX = sqrt(2)Z, CX = sqrt(2)Z]
    C = 1 / (2 * pi * f) / Z * sqrt(2)
    L = 1 / (2 * pi * f) * Z * sqrt(2)
    return (L, C)
end

function PlotPhasor(
    Impedance::Tuple,
    PrevImpedance;
    Color = nothing,
    NumSteps = 1,
    Admittance = false,
    HeadWidth = 0.01,
)

    ## This function takes in an impedance, and plots it using 'phasor()'
    #optionally, it takes in a previous impedance which can be added in series or parallel
    #It then plots the line using PrevImpedance as a base
    #If Admittance is true, it means the inputs are Admittance not impedance, so parallel and series are switched
    NewImp = 0 + 0im
    dx = 0
    dy = 0
    for I = 1:NumSteps
        if I > 1
            #    Color = "gray"
        end
        if lowercase(Impedance[2]) == "series"
            if !Admittance
                NewImp = PrevImpedance + Impedance[1] / NumSteps
            else
                NewImp = Par(PrevImpedance, Impedance[1] * NumSteps)
            end
        elseif lowercase(Impedance[2]) == "parallel"
            if !Admittance
                NewImp = Par(PrevImpedance, Impedance[1] * NumSteps)
            else
                NewImp = PrevImpedance + Impedance[1] / NumSteps
            end
        else
            error("must be series or parallel")
        end

        RealParts = [real(NewImp), real(PrevImpedance)]
        ImagParts = [imag(NewImp), imag(PrevImpedance)]

        dx = -(RealParts[2] - RealParts[1])
        dy = -(ImagParts[2] - ImagParts[1])
        if Color !== nothing

            PyPlot.plot(RealParts, ImagParts, color = "gray")

            if I == 1
                phasor(PrevImpedance, color = Color)
            end
        else
            PyPlot.plot(RealParts, ImagParts)
            if I == 1
                phasor(PrevImpedance)
            end
        end
        PrevImpedance = NewImp
    end
    # arrow(real(NewImp)-dx,imag(NewImp)-dy,dx,dy,head_width=HeadWidth,fc="gray",ec="k",alpha=1.,width = 0.00001,head_length= 2*HeadWidth)

    xlabel("Real")
    ylabel("Imag.")
    return NewImp
end

"""
This function takes in an array of tuples structures such as:

    [(<Complex impedance #1> , <"series" or "Parallel">),
    (<Complex impedance #2> , <"series" or "Parallel">)
    .
    .
    .
    (<Complex impedance #n> , <"series" or "Parallel">)]

    e.g.

    ZList = [(1+1*im, "series"),
             (1+5*im, "series"),
             (5+1*im, "parallel"),]
    PlotImpedanceTransformList(ZList)
"""
function PlotImpedanceTransformList(ImpList; InitialImp = nothing, ArrHeadWidth =1)
    ColorList = ["r", "g", "b", "magenta", "teal"]
    L = zeros(length(ImpList))
    for I = 1:length(ImpList)
        L[I] = abs(ImpList[I][1])
    end
    LongestVector = maximum(L)
    HeadWidth = LongestVector / ArrHeadWidth
    ##Impedance Plotting
    if InitialImp !== nothing
        PrevImp = InitialImp
        StartVal = 1
    else
        PrevImp = ImpList[1][1]
        StartVal = 2
    end
    for I = StartVal:length(ImpList)
        PrevImp = PlotPhasor(
            ImpList[I],
            PrevImp;
            NumSteps = 500,
            Color = ColorList[I],
            HeadWidth = HeadWidth,
        )
    end
    phasor(PrevImp, color = ColorList[length(ImpList)+1])
    title("Impedance")
    ##Admittance Plotting
    figure()

    L = zeros(length(ImpList))
    for I = 1:length(ImpList)
        L[I] = abs(1 / ImpList[I][1])
    end
    LongestVector = maximum(L)
    HeadWidth = LongestVector / ArrHeadWidth
    if InitialImp !=0 nothing
        PrevImp = 1 / InitialImp
        StartVal = 1
    else
        PrevImp = 1 / ImpList[1][1]
        StartVal = 2
    end
    for I = StartVal:length(ImpList)
        PrevImp = PlotPhasor(
            (1 / ImpList[I][1], ImpList[I][2]),
            PrevImp;
            NumSteps = 500,
            Color = ColorList[I],
            Admittance = true,
            HeadWidth = HeadWidth,
        )
    end
    phasor(PrevImp, color = ColorList[length(ImpList)+1])
    title("Admittance")
end

function findMinima_Bounds(x,y;StartVal=nothing,StopVal=nothing)

    if !(StartVal===nothing)
        StartIndex = minimum(findall(x .>= StartVal))
    else
        StartIndex = 1
    end

    if !(StopVal===nothing)
        StopIndex = maximum(findall(x .<= StopVal))
    else
        StopIndex = length(x)
    end

    xSubset = x[StartIndex:StopIndex]
    ySubset = y[StartIndex:StopIndex]
    FirstDiff = diff(ySubset)
    SecondDiff = diff(FirstDiff)
    Tmp = SecondDiff .< 0
    push!(Tmp,1)
    FirstDiff[Tmp] .= 1000 #Only look at values where the second derivative is positive

    minIndex = findfirst(isequal(minimum(abs.(FirstDiff))),abs.(FirstDiff))

    Min_x = xSubset[minIndex]
    Min_y = ySubset[minIndex]
    return Min_x, Min_y

end

function findDriveDistFromTrough(fVec,CurVec,DriveFreq)
    TroughFreq,val = findMinima_Bounds(fVec,abs.(CurVec);StartVal=(DriveFreq-1e3),StopVal=(DriveFreq+1e3))
    return abs(TroughFreq-DriveFreq)
end

function plotDriveFreqDot(fVec,CurVec,DriveFreq)
    DriveFreqCurrentIndex = findfirst(fVec.>=DriveFreq)
    plot(fVec[DriveFreqCurrentIndex],20*log10.(abs.(CurVec[DriveFreqCurrentIndex])),"r*")
end

# """
# This function takes in similar inputs to DesignDriveFilter. It will iterate that function with varying matching impedances to minimize the expected drift in the current.
# """
# function DesignDriveFilter_OptimDrift(
#     LDrive,
#     RDrive,
#     TargetZ,
#     DriveFreq;
#     CDrive = 2,
#     NumDriveElements = 1,
#     WireDiam = 2e-3,
#     WireFillFac = 0.75,
#     PlotSrcR = TargetZ,
#     PlotFTs = true,
#     VSrc = 2,
#     DetermineFreq = false,
#     AddNotchFreq = nothing,
#     FilterZ = TargetZ,
#     RDampVal = FilterZ,
#     PerturbTxReactance = nothing,
#     BruteForceOpt = false,
#     MinimizeDistToTrough = true,
#     OptimIters = 10,
#     WriteFileName=nothing
# )

# AxMain = PyPlot.axes([0.13, 0.15, 0.85, 0.85])
# AxInset = PyPlot.axes([0.68, 0.65, 0.25, 0.25]) #[left, bottom, width, height]

# AxesArray = [AxMain,AxInset]

#         if MinimizeDistToTrough
#             println("minimizing distance to trough")
#         else
#             println("Minimizing drift coeff")
#         end
#         function WrapperFunkZin(Zin)
#             DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList = DesignDriveFilter(
#             LDrive,RDrive, Zin,DriveFreq;
#             CDrive = CDrive,
#             NumDriveElements = NumDriveElements,
#             WireDiam = WireDiam,
#             DetermineFreq=DetermineFreq,
#             AddNotchFreq = AddNotchFreq,
#             FilterZ = FilterZ,
#             AxesArray=AxesArray)

#             DetermineComponentsTempCoeffs(SPICE_DF,InputList,1,DriveFreq,"LDrive")
#             if CDrive>1 #Check the drift in the drive capacitor if it is used, otherwise use the CSer
#                 Drift = SPICE_DF.DriftCoeff[findfirst(isequal("CSer"),SPICE_DF.Name)]
#                 PhaseDrift = SPICE_DF.DriftCoeffPhase[findfirst(isequal("CSer"),SPICE_DF.Name)]
#             else
#                 Drift = SPICE_DF.DriftCoeff[findfirst(isequal("CDrive"),SPICE_DF.Name)]
#                 PhaseDrift = SPICE_DF.DriftCoeffPhase[findfirst(isequal("CDrive"),SPICE_DF.Name)]
#             end
#             TroughDist = findDriveDistFromTrough(FreqList,CurrentVec,DriveFreq)
#             println("$TroughDist is the distance from the drive frequency to the trough, $(round(Drift;sigdigits=3)) is the drift in magnitude, $(round(PhaseDrift;sigdigits=3)) is the drift in phase")
#             if MinimizeDistToTrough
#                 return TroughDist
#             else
#                 return abs(Drift)
#             end
#         end

#         function WrapperFunk_Per(PerturbationX)
#             DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList = DesignDriveFilter(
#             LDrive,RDrive, TargetZ,DriveFreq;
#             CDrive = CDrive,
#             NumDriveElements = NumDriveElements,
#             WireDiam = WireDiam,
#             DetermineFreq=DetermineFreq,
#             AddNotchFreq = AddNotchFreq,
#             FilterZ = TargetZ,
#             PerturbTxReactance = PerturbationX,AxesArray=AxesArray)

#             DetermineComponentsTempCoeffs(SPICE_DF,InputList,1,DriveFreq,"LDrive")
#             if CDrive>1 #Check the drift in the drive capacitor if it is used, otherwise use the CSer
#                 Drift = SPICE_DF.DriftCoeff[findfirst(isequal("CSer"),SPICE_DF.Name)]
#                 PhaseDrift = SPICE_DF.DriftCoeffPhase[findfirst(isequal("CSer"),SPICE_DF.Name)]
#             else
#                 Drift = SPICE_DF.DriftCoeff[findfirst(isequal("CDrive"),SPICE_DF.Name)]
#                 PhaseDrift = SPICE_DF.DriftCoeffPhase[findfirst(isequal("CDrive"),SPICE_DF.Name)]
#             end
#             TroughDist = findDriveDistFromTrough(FreqList,CurrentVec,DriveFreq)
#             println("$TroughDist is the distance from the drive frequency to the trough, $(round(Drift;sigdigits=3)) is the drift in magnitude, $(round(PhaseDrift;sigdigits=3)) is the drift in phase")
#             if MinimizeDistToTrough
#                 return TroughDist
#             else
#                 return abs(Drift)
#             end
#         end

#     if  (PerturbTxReactance === nothing)
#         MinZ = optimize((WrapperFunkZin),5,25,iterations = OptimIters)
#         DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList = DesignDriveFilter(
#         LDrive,RDrive, MinZ.minimizer,DriveFreq;
#         CDrive = CDrive,
#         WireDiam = WireDiam,
#         DetermineFreq=DetermineFreq,
#         AddNotchFreq = AddNotchFreq,
#         FilterZ = FilterZ,
#         WriteFileName=WriteFileName)

#         MinZVal = MinZ.minimizer
#         DetermineComponentsTempCoeffs(SPICE_DF,InputList,1,DriveFreq,"LDrive")
#         return MinZVal, DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList,AxesArray
#     else
#         println("Perturbation given")
#         if BruteForceOpt

#             PerVec = -2:.1:2
#             DriftVec = zeros(length(PerVec),1)
#             for i in 1:length(PerVec)
#                 println(i)
#                 DriftVec[i] =WrapperFunk_Per(PerVec[i])
#             end
#             MinIndex = findfirst(x-> x==minimum(DriftVec),DriftVec)
#             MinXVal = PerVec[MinIndex]
#         else

#             MinXValOptimRes = optimize((WrapperFunk_Per),-2,2,iterations = OptimIters)
#             MinXVal = MinXValOptimRes.minimizer
#         end
#         DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList = DesignDriveFilter(
#         LDrive,RDrive, TargetZ,DriveFreq;
#         CDrive = CDrive,
#         WireDiam = WireDiam,
#         DetermineFreq=DetermineFreq,
#         AddNotchFreq = AddNotchFreq,
#         FilterZ = TargetZ,
#         PerturbTxReactance = MinXVal,
#         WriteFileName=WriteFileName)

#         DetermineComponentsTempCoeffs(SPICE_DF,InputList,1,DriveFreq,"LDrive")
#         return MinXVal, DriveFreq, CurrentVec, Results, SPICE_DF,inputs,InputList,FreqList,AxesArray
#     end

# end


text:
@recipe(Coastlines) do scene
    merge(
        default_theme(scene, Lines),
        Theme(crs = (source = LonLat(), dest = LonLat()),)
    )
end

function AbstractPlotting.plot!(p::Coastlines)

    @extract p (crs,)

    new_linevec = Node(COASTLINES_LINEVEC)

    lift(crs.source, crs.dest) do source, dest
        new_linevec[] = Proj4.transform.(source, dest, COASTLINES_LINEVEC)
    end

    lines!(p, new_linevec)
end

coastlines(; kwargs...) = coastlines(1; kwargs...)
coastlines!(sc = AbstractPlotting.current_scene(); kwargs...) = plot!(sc, Coastlines, Attributes(attributes), 1)
coastlines!(; attributes...) = plot!(AbstractPlotting.current_scene(), Coastlines, Attributes(attributes), 1)

@recipe(Earth, bbox) do scene
    merge(
        default_theme(scene, Surface),
        Theme(
            crs = (
                source = LonLat(),
                dest = LonLat()
            ),
            shading = false
        )
    )
end

function AbstractPlotting.plot!(p::Earth)

    @extract p (crs,)

    lons = LinRange(-180, 179, size(EARTH_IMG)[2])
    lats = LinRange(89.5, -89.5, size(EARTH_IMG)[1])

    xs = Node([lon for lat in lats, lon in lons])
    ys = Node([lat for lat in lats, lon in lons])

    lift(crs.source, crs.dest) do source, dest

        xs.val = [lon for lat in lats, lon in lons]
        ys.val = [lat for lat in lats, lon in lons]
        Proj4.transform!(source, dest, vec(xs.val), vec(ys.val))

        xs[] = xs[]
        ys[] = ys[]
    end

    surface!(p, xs, ys, zeros(size(xs[])); color = EARTH_IMG, shading = false, show_axis = false)
end

earth(; kwargs...) = earth(1; kwargs...)
earth!(sc = AbstractPlotting.current_scene(); attributes...) = plot!(sc, Earth, Attributes(attributes), 1)
earth!(; attributes...) = plot!(AbstractPlotting.current_scene(), Earth, Attributes(attributes), 1)

ext: jl


text:
module LightDOM

import Base: convert,
	length,
	push!

export 	Node,
	TextNode,
	Element,
	text,
	count,
	namespace,
	tag,
	children,
	props

#
# Node base
#

abstract type Node end

#
# Text Node
#

struct TextNode <: Node
	text::AbstractString
end

"""
	text("This is Text)

	Creates a TextNode
"""
text(xs...) = TextNode(string(xs...))

"""
	convert a string to a Node
"""
Base.convert(::Type{Node}, s::AbstractString) = text(s)

#
# Generic Element
#

const Props = Dict{Any,Any}

"""
	Namespace and tag are part of the type (rather than a prop) to allow
	dispatch to provide overrides of methods for specific namespaces/tags
"""
struct Element{ns, tag} <: Node
	children::Vector{Union{TextNode,Element}}
	props::Props
end
Element{ns, tag}() where {ns, tag} = Element{ns,tag}([], Dict())
Element(ns::Symbol, tag::Symbol, children, props) = Element{ns, tag}(children, props)
Element(ns::Symbol, tag::Symbol, children=[]; kwargs...) = Element(ns, tag, children, Props(kwargs))
Element(tag::Symbol, children::AbstractVector=[]; kwargs...) = Element(:nothing, tag, children, Props(kwargs))
Element(tag::Symbol, textContent::AbstractString; kwargs...) = Element(:nothing, tag, text(textContent), Props(kwargs))
Element(tag::Symbol, child::Node; kwargs...) = Element(:nothing, tag, child, Props(kwargs))

Base.convert(::Type{Vector{Union{TextNode,Element}}}, x::Node) = Vector{Union{TextNode,Element}}([x])

function namespace(e::Element{ns}) where ns
	ns
end
function tag(e::Element{ns,t}) where {ns, t}
	t
end

include("show.jl")
include("collection.jl")

end # module


text:
# dbn_proposals.jl
# 2019-11-12
# David Merrell
#
# Some proposal distributions for exploring the space
# of directed graphs.

############################################
# PARENT SET PROPOSAL DISTRIBUTION
############################################
function compute_move_probs(indeg::Int64, degmax::Int64, prob_param::Float64)
    frac = (indeg/degmax)^prob_param
    comp = 1.0 - frac
    swp = 2.0 * frac * comp
    return [comp; swp; frac] ./ (frac + swp + comp)
end

function ith_nonparent(i, parents::Vector{Int64})

    n = length(parents)
    if n == 0
        return i
    end

    p = 1
    while p <= n
    if (parents[p] - p) >= i
        break
        end
    p += 1
    end
    return i + p - 1
end

"""
Return the index of x in the *complement* of a;
that is, the index of x in the array [1,2,...] setminus a.
If x is in a, then we return the index it *would* have
if it were moved into a.
"""
function searchsorted_exclude(a::Vector{Int64}, x::Int64)
    r = searchsorted(a, x)
    if r.stop < r.start
        return x - r.stop
    else
        return x - r.stop + 1
    end
end

"""
Parent set proposal distribution.

This proposal allows "add parent", "remove parent", and "swap parent" moves.

"Swap" moves preserve in-degree. This is a good thing, since our
model's posterior distribution varies pretty strongly with in-degree.
"""
@gen function smart_proposal(tr, vertex::Int64,
                             prob_param::Float64,
                             V::Int64)

    parents = copy(tr[:parent_sets => vertex => :parents])
    indeg = length(parents)
    action_probs = compute_move_probs(indeg, V, prob_param)
    action = @trace(Gen.categorical(action_probs), :action)

    if action == 1
        to_add = @trace(Gen.uniform_discrete(1, V-indeg), :to_add)
        return (parents, action, to_add)
    elseif action == 2
        to_swp_add = @trace(Gen.uniform_discrete(1, V-indeg), :to_swp_add)
        to_swp_rem = @trace(Gen.uniform_discrete(1, indeg), :to_swp_rem)
        return (parents, action, to_swp_add, to_swp_rem)
    else
        to_rem = @trace(Gen.uniform_discrete(1, indeg), :to_rem)
        return (parents, action, to_rem)
    end
end

function smart_involution(cur_tr, fwd_choices, fwd_ret, prop_args)

    parents = fwd_ret[1]
    action = fwd_ret[2]
    vertex = prop_args[1]

    update_choices = Gen.choicemap()
    bwd_choices = Gen.choicemap()

    if action == 1
        to_add = fwd_ret[3]
        add_parent = ith_nonparent(to_add, parents)
        insert_idx = searchsortedfirst(parents, add_parent)
        insert!(parents, insert_idx, add_parent)
        update_choices[:parent_sets => vertex => :parents] = parents
        bwd_choices[:action] = 3
        bwd_choices[:to_rem] = insert_idx

    elseif action == 2
        to_swp_add = fwd_ret[3]
        to_swp_rem = fwd_ret[4]
        add_parent = ith_nonparent(to_swp_add, parents)
        rem_parent = parents[to_swp_rem]

        deleteat!(parents, to_swp_rem)
        insert_idx = searchsortedfirst(parents, add_parent)
        insert!(parents, insert_idx, add_parent)

        update_choices[:parent_sets => vertex => :parents] = parents

        bwd_choices[:action] = 2
        bwd_choices[:to_swp_add] = searchsorted_exclude(parents, rem_parent)
        bwd_choices[:to_swp_rem] = insert_idx

    else
        to_rem = fwd_ret[3]
        rem_parent = parents[to_rem]
        deleteat!(parents, to_rem)
        update_choices[:parent_sets => vertex => :parents] = parents
        bwd_choices[:action] = 1
        bwd_choices[:to_add] = searchsorted_exclude(parents, rem_parent)
    end

    new_tr, weight, _, _ = Gen.update(cur_tr, Gen.get_args(cur_tr),
                                      (), update_choices)

    return new_tr, bwd_choices, weight
end

####################################
# UNIFORM GRAPH PROPOSAL
####################################

"""
In order to propose reverse-edge moves,
we need to know which of the edges are not paired.
"""
function unpaired_edges(tr, V)
    result = Set()
    for child = 1:V
        for p in tr[:parent_sets => child => :parents]
            if (child, p) in result
                delete!(result, (child,p))
            elseif p != child
                push!(result, (p,child))
            end
        end
    end
    return sort(collect(result))
end

"""
Uniform graph proposal distribution.
"""
@gen function uniform_proposal(tr, V)

    ue = unpaired_edges(tr,V)
    bound = V^2 + length(ue)
    idx = @trace(Gen.uniform_discrete(1, bound), :idx)

    if idx <= V^2
        u_idx = div(idx-1, V) + 1
        v_idx = ((idx-1) % V) + 1

        if u_idx in tr[:parent_sets => v_idx => :parents]
            return (u_idx, v_idx), bound, ue, "remove"
        else
            return (u_idx, v_idx), bound, ue, "add"
        end
    else
        return (ue[idx - V^2]), bound, ue, "reverse"
    end
end

"""
Involution for the uniform graph proposal
"""
function uniform_involution(cur_tr, fwd_choices, fwd_ret, prop_args)

    idx = fwd_ret[1]
    u_idx = idx[1]
    v_idx = idx[2]
    bound = fwd_ret[2]
    ue = fwd_ret[3]
    action = fwd_ret[4]
    V = prop_args[1]

    update_choices = Gen.choicemap()
    bwd_choices = Gen.choicemap()

    if action == "add"
        bwd_choices[:idx] = (u_idx-1)*V + v_idx

        v_parents = copy(cur_tr[:parent_sets => v_idx => :parents])
        insert_idx = searchsortedfirst(v_parents, u_idx)
        insert!(v_parents, insert_idx, u_idx)
        update_choices[:parent_sets => v_idx => :parents] = v_parents

    elseif action == "remove"
        bwd_choices[:idx] = (u_idx-1)*V + v_idx

        v_parents = copy(cur_tr[:parent_sets => v_idx => :parents])
        rem_idx = searchsortedfirst(v_parents, u_idx)
        deleteat!(v_parents, rem_idx)
        update_choices[:parent_sets => v_idx => :parents] = v_parents

    elseif action == "reverse"

        ue_idx = searchsortedfirst(ue, (u_idx, v_idx))
        deleteat!(ue, ue_idx)

        bwd_choices[:idx] = V^2 + searchsortedfirst(ue, (v_idx, u_idx))

        v_parents = copy(cur_tr[:parent_sets => v_idx => :parents])
        rem_idx = searchsortedfirst(v_parents, u_idx)
        deleteat!(v_parents, rem_idx)
        update_choices[:parent_sets => v_idx => :parents] = v_parents

        u_parents = copy(cur_tr[:parent_sets => u_idx => :parents])
        insert_idx = searchsortedfirst(u_parents, v_idx)
        insert!(u_parents, insert_idx, v_idx)
        update_choices[:parent_sets => u_idx => :parents] = u_parents
    end

    new_tr, weight, _, _ = Gen.update(cur_tr, Gen.get_args(cur_tr),
                                      (), update_choices)

    return new_tr, bwd_choices, weight
end

@gen function lambda_proposal(tr, std::Int64)

    @trace(Gen.normal(tr[:lambda], std), :lambda)

end

@gen function lambda_vec_proposal(tr, i::Int64, std::Int64)

    @trace(Gen.normal(tr[:lambda_vec => i => :lambda], std), :lambda_vec => i => :lambda)

end

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
struct Queue
    store::Vector
end

q = Queue(Vector{Float64}(undef, 3))
println(q.store)


text:
using Augmentor
using Images
using TestImages
using VisualRegressionTests
#using Plots

#ENV["MPLBACKEND"] = "Agg"
#try
#    @eval import PyPlot
#    info("Matplotlib version: $(PyPlot.matplotlib[:__version__])")
#end

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

refdir = Pkg.dir("Augmentor", "test", "refimg")
testimg = load(joinpath(refdir, "testimage.png"))

function imagetest_impl(testname, testfun)
    srand(1)
    refimgpath = joinpath(refdir, "$testname.png")
    result = test_images(VisualTest(testfun, refimgpath))
    @test success(result)
end

macro plottest(testname, expr)
    esc(quote
        imagetest_impl(string($testname), fn -> begin
            $expr
            png(fn)
        end)
    end)
end

macro imagetest(testname, expr)
    esc(quote
        imagetest_impl(string($testname), path -> begin
            save(path, $expr)
        end)
    end)
end

type FaultyOp <: Augmentor.ImageTransformation end

tests = [
    ("tst_displacementmesh.jl", "Low-level functionality for image displacement"),
    ("tst_displacementfield.jl", "Low-level functionality for image displacement"),
    ("tst_common.jl", "Utility methods"),
    ("tst_rotate.jl", "Low-level functionality for image rotation"),
    ("tst_crop.jl", "Low-level functionality for image cropping"),
    ("tst_imagetransformations.jl", "High-level abstractions for image transformations"),
    ("tst_imagesource.jl", "High-level abstractions for image sources"),
    ("tst_pipeline.jl", "Pipeline implementations"),
    ("tst_MLDataUtils.jl", "Data Access Pattern integration"),
]

for (fn, desc) in tests
    @testset "$desc ($fn)" begin
        include(fn)
    end
end


text:
module InfixFunctions

export @infix

struct InfixFunction <: Function
    name::Symbol
    operator::Function
end

InfixFunction(operator::Function) = InfixFunction(gensym(), operator)

(infix::InfixFunction)(arg₁, arg₂) = infix.operator(arg₁, arg₂)

function Base.show(io::IO, infix::InfixFunction)
    n_methods = length(methods(infix.operator))
    _methods = n_methods == 1 ? "method" : "methods"
    name = infix.name

    println(io, "$name (generic infix function with $n_methods $_methods)")
end

Base.display(infix::InfixFunction) = show(infix)

function Base.:|(arg₁, infix::InfixFunction)
    return InfixFunction(arg₂ -> infix.operator(arg₁, arg₂))
end

Base.:|(infix::InfixFunction, arg₂) = infix.operator(arg₂)

macro infix(operator::Symbol)
    return quote
        $operator::Function

        function Base.:|(arg₁, infix::typeof($operator))
            return $InfixFunction(arg₂ -> infix(arg₁, arg₂))
        end

        Base.:|(infix::typeof($operator), arg₂) = infix(arg₂)

        @info "$($operator) has been infixified"

        $operator
    end |> esc
end

end


text:
module VectorSpaces

using ..Traits: rawtype, efficientoperations, getcontent

import ..Traits: contentnames
import ...Interfaces: dimension, rank, ⊕, ⊗

export VectorSpace, NamedVectorSpace, SimpleNamedVectorSpace, ParameterSpace, ZippedNamedVectorSpace, DirectProductedNamedVectorSpace
export VectorSpaceStyle, VectorSpaceEnumerative, VectorSpaceCartesian, VectorSpaceDirectSummed, VectorSpaceDirectProducted, VectorSpaceZipped
export shape

"""
    VectorSpace{B} <: AbstractVector{B}

Abstract vector space.
"""
abstract type VectorSpace{B} <: AbstractVector{B} end
@inline Base.:(==)(vs₁::VectorSpace, vs₂::VectorSpace) = ==(efficientoperations, vs₁, vs₂)
@inline Base.isequal(vs₁::VectorSpace, vs₂::VectorSpace) = isequal(efficientoperations, vs₁, vs₂)
@inline Base.size(vs::VectorSpace) = (dimension(vs),)

"""
    VectorSpaceStyle

The style of a concrete type of vector space.
"""
abstract type VectorSpaceStyle end
@inline VectorSpaceStyle(vs::VectorSpace) = VectorSpaceStyle(typeof(vs))

@inline dimension(vs::VectorSpace) = dimension(VectorSpaceStyle(vs), vs)
@inline Base.getindex(vs::VectorSpace, i) = getindex(VectorSpaceStyle(vs), vs, i)
@inline Base.issorted(vs::VectorSpace) = issorted(VectorSpaceStyle(vs), vs)
@inline Base.findfirst(basis::B, vs::VectorSpace{B}) where B = findfirst(VectorSpaceStyle(vs), basis, vs)
@inline Base.searchsortedfirst(vs::VectorSpace{B}, basis::B) where B = searchsortedfirst(VectorSpaceStyle(vs), vs, basis)
@inline Base.in(basis::B, vs::VectorSpace{B}) where B = in(VectorSpaceStyle(vs), basis, vs)

@inline Base.issorted(::VectorSpaceStyle, vs::VectorSpace) = false
@inline function Base.findfirst(::VectorSpaceStyle, basis::B, vs::VectorSpace{B}) where B
    if issorted(vs)
        index = invoke(searchsortedfirst, Tuple{Any, Any}, vs, basis)
        return (0<index<=length(vs) && vs[index]==basis) ? index : nothing
    end
    for i = 1:length(vs)
        vs[i]==basis && return i
    end
end
@inline function Base.searchsortedfirst(::VectorSpaceStyle, vs::VectorSpace{B}, basis::B) where B
    issorted(vs) && return invoke(searchsortedfirst, Tuple{Any, Any}, vs, basis)
    for i = 1:length(vs)
        vs[i]==basis && return i
    end
    return length(vs)+1
end
@inline Base.in(::VectorSpaceStyle, basis::B, vs::VectorSpace{B}) where B = isa(findfirst(basis, vs), Integer)

"""
    VectorSpaceEnumerative <: VectorSpaceStyle

Enumerative vector space style, which indicates that the vector space has a predefined content named `:table` that contains all its bases.
"""
struct VectorSpaceEnumerative <: VectorSpaceStyle end
@inline dimension(::VectorSpaceEnumerative, vs::VectorSpace) = length(getcontent(vs, :table))
@inline Base.getindex(::VectorSpaceEnumerative, vs::VectorSpace, i::Int) = getcontent(vs, :table)[i]

"""
    VectorSpaceCartesian <: VectorSpaceStyle

Cartesian vector space style, which indicates that every basis in it could be represented by a Cartesian index.
"""
struct VectorSpaceCartesian <: VectorSpaceStyle end
function shape end
@inline dimension(::VectorSpaceCartesian, vs::VectorSpace) = mapreduce(length, *, shape(vs))
@inline Base.getindex(::VectorSpaceCartesian, vs::VectorSpace, i::Integer) = getindex(VectorSpaceCartesian(), vs, CartesianIndices(shape(vs))[i])
@inline Base.getindex(::VectorSpaceCartesian, vs::VectorSpace, i::CartesianIndex) = rawtype(eltype(vs))(i, vs)
@inline Base.issorted(::VectorSpaceCartesian, vs::VectorSpace) = true
@inline function Base.findfirst(::VectorSpaceCartesian, basis::B, vs::VectorSpace{B}) where B
    index, idx = CartesianIndex(basis, vs), CartesianIndices(shape(vs))
    index∉idx && return nothing
    return LinearIndices(shape(vs))[index-first(idx)+oneunit(typeof(index))]
end
@inline function Base.searchsortedfirst(::VectorSpaceCartesian, vs::VectorSpace{B}, basis::B) where B
    index, idx = CartesianIndex(basis, vs), CartesianIndices(shape(vs))
    index∈idx && return LinearIndices(shape(vs))[index-first(idx)+oneunit(typeof(index))]
    return searchsortedfirst(idx, index)
end
@inline Base.in(::VectorSpaceCartesian, basis::B, vs::VectorSpace{B}) where B = CartesianIndex(basis, vs)∈CartesianIndices(shape(vs))

"""
    VectorSpaceDirectSummed <: VectorSpaceStyle

Vector space style which indicates that a vector space is the direct sum of its sub-components.
"""
struct VectorSpaceDirectSummed <: VectorSpaceStyle end
@inline dimension(::VectorSpaceDirectSummed, vs::VectorSpace) = mapreduce(dimension, +, getcontent(vs, :contents))
@inline function Base.getindex(::VectorSpaceDirectSummed, vs::VectorSpace, i::Integer)
    contents = getcontent(vs, :contents)
    dimsum = cumsum(map(dimension, contents))
    m = searchsortedfirst(dimsum, i)
    n = i-(m>1 ? dimsum[m-1] : 0)
    return contents[m][n]
end

"""
    VectorSpaceDirectProducted <: VectorSpaceStyle

Vector space style which indicates that a vector space is the direct product of its sub-components.
"""
struct VectorSpaceDirectProducted <: VectorSpaceStyle end
@inline dimension(::VectorSpaceDirectProducted, vs::VectorSpace) = mapreduce(dimension, *, getcontent(vs, :contents))
@inline function Base.getindex(::VectorSpaceDirectProducted, vs::VectorSpace, i::Integer)
    contents = getcontent(vs, :contents)
    index = CartesianIndices(map(dimension, contents))[i]
    return rawtype(eltype(vs))(map(getindex, contents, Tuple(index)), vs)
end

"""
    VectorSpaceZipped <: VectorSpaceStyle

Vector space style which indicates that a vector space is the zip of its sub-components.
"""
struct VectorSpaceZipped <: VectorSpaceStyle end
@inline dimension(::VectorSpaceZipped, vs::VectorSpace) = mapreduce(dimension, min, getcontent(vs, :contents))
@inline function Base.getindex(::VectorSpaceZipped, vs::VectorSpace, i::Integer)
    return rawtype(eltype(vs))(map(content->getindex(content, i), getcontent(vs, :contents)), vs)
end

"""
    NamedVectorSpace{B} <: VectorSpace{B}

Abstract named vector space.
"""
abstract type NamedVectorSpace{B} <: VectorSpace{B} end
@inline Base.keys(vs::NamedVectorSpace) = keys(typeof(vs))
@inline Base.:(==)(vs₁::NamedVectorSpace, vs₂::NamedVectorSpace) = keys(vs₁)==keys(vs₂) && invoke(==, Tuple{VectorSpace, VectorSpace}, vs₁, vs₂)
@inline Base.isequal(vs₁::NamedVectorSpace, vs₂::NamedVectorSpace) = isequal(keys(vs₁), keys(vs₂)) && invoke(isequal, Tuple{VectorSpace, VectorSpace}, vs₁, vs₂)
@inline Base.pairs(vs::NamedVectorSpace) = NamedVectorSpacePairIteration(vs)
struct NamedVectorSpacePairIteration{VS<:NamedVectorSpace} <: AbstractVector{NamedTuple}
    vs::VS
end
@inline Base.size(ps::NamedVectorSpacePairIteration) = size(ps.vs)
@inline Base.eltype(::Type{<:NamedVectorSpacePairIteration{VS}}) where {VS<:NamedVectorSpace} = NamedTuple{keys(VS), eltype(VS)}
@inline Base.getindex(ps::NamedVectorSpacePairIteration{VS}, i::Integer) where {VS<:NamedVectorSpace} = NamedTuple{keys(VS)}(ps.vs[i])

"""
    SimpleNamedVectorSpace{N, B} <: NamedVectorSpace{B}

Abstract simple named vector space.

The subtypes must have the following predefined content:
* `:content::T`: the content of the named vector space.
"""
abstract type SimpleNamedVectorSpace{N, B} <: NamedVectorSpace{B} end
@inline contentnames(::Type{<:SimpleNamedVectorSpace}) = (:content,)
@inline dimension(vs::SimpleNamedVectorSpace) = length(getcontent(vs, :content))
@inline Base.getindex(vs::SimpleNamedVectorSpace, i) = getcontent(vs, :content)[i]
@inline Base.in(basis::B, vs::SimpleNamedVectorSpace{N, B}) where {N, B} = in(basis, getcontent(vs, :content))
@inline Base.keys(::Type{<:SimpleNamedVectorSpace{N}}) where N = (N,)
@inline Base.eltype(::Type{<:NamedVectorSpacePairIteration{VS}}) where {VS<:SimpleNamedVectorSpace} = NamedTuple{keys(VS), Tuple{eltype(VS)}}
@inline Base.getindex(ps::NamedVectorSpacePairIteration{VS}, i::Integer) where {VS<:SimpleNamedVectorSpace} = NamedTuple{keys(VS)}((ps.vs[i],))

"""
    ParameterSpace{N, T, B} <: SimpleNamedVectorSpace{N, B}

Parameter space.
"""
struct ParameterSpace{N, T, B} <: SimpleNamedVectorSpace{N, B}
    content::T
    function ParameterSpace{N}(content) where N
        @assert isa(N, Symbol) "ParameterSpace error: $(N) is not a Symbol."
        new{N, typeof(content), eltype(content)}(content)
    end
end

"""
    ZippedNamedVectorSpace{T<:Tuple{Vararg{NamedVectorSpace}}} <: NamedVectorSpace{Any}

Zipped named vector space.
"""
struct ZippedNamedVectorSpace{T<:Tuple{Vararg{NamedVectorSpace}}} <: NamedVectorSpace{Any}
    contents::T
end
@inline ZippedNamedVectorSpace(contents::NamedVectorSpace...) = ZippedNamedVectorSpace(contents)
@inline VectorSpaceStyle(::Type{<:ZippedNamedVectorSpace}) = VectorSpaceZipped()
@inline @generated Base.eltype(::Type{<:ZippedNamedVectorSpace{T}}) where {T<:Tuple{Vararg{NamedVectorSpace}}} = Tuple{map(eltype, fieldtypes(T))...}
@inline @generated Base.keys(::Type{<:ZippedNamedVectorSpace{T}}) where {T<:Tuple{Vararg{NamedVectorSpace}}} = map(first, map(keys, fieldtypes(T)))
@inline Tuple(basis::Tuple, ::ZippedNamedVectorSpace) = basis

"""
    DirectProductedNamedVectorSpace{T<:Tuple{Vararg{NamedVectorSpace}}} <: NamedVectorSpace{Any}

Zipped named vector space.
"""
struct DirectProductedNamedVectorSpace{T<:Tuple{Vararg{NamedVectorSpace}}} <: NamedVectorSpace{Any}
    contents::T
end
@inline DirectProductedNamedVectorSpace(contents::NamedVectorSpace...) = DirectProductedNamedVectorSpace(contents)
@inline VectorSpaceStyle(::Type{<:DirectProductedNamedVectorSpace}) = VectorSpaceDirectProducted()
@inline @generated Base.eltype(::Type{<:DirectProductedNamedVectorSpace{T}}) where {T<:Tuple{Vararg{NamedVectorSpace}}} = Tuple{map(eltype, fieldtypes(T))...}
@inline @generated Base.keys(::Type{<:DirectProductedNamedVectorSpace{T}}) where {T<:Tuple{Vararg{NamedVectorSpace}}} = map(first, map(keys, fieldtypes(T)))
@inline Tuple(basis::Tuple, ::DirectProductedNamedVectorSpace) = basis

"""
    ⊕(vs₁::SimpleNamedVectorSpace, vs₂::SimpleNamedVectorSpace) -> ZipNamedVectorSpace
    ⊕(vs₁::SimpleNamedVectorSpace, vs₂::ZippedNamedVectorSpace) -> ZipNamedVectorSpace
    ⊕(vs₁::ZippedNamedVectorSpace, vs₂::SimpleNamedVectorSpace) -> ZipNamedVectorSpace
    ⊕(vs₁::ZippedNamedVectorSpace, vs₂::ZippedNamedVectorSpace) -> ZipNamedVectorSpace

The zip of named vector spaces.
"""
@inline ⊕(vs₁::SimpleNamedVectorSpace, vs₂::SimpleNamedVectorSpace) = ZippedNamedVectorSpace(vs₁, vs₂)
@inline ⊕(vs₁::SimpleNamedVectorSpace, vs₂::ZippedNamedVectorSpace) = ZippedNamedVectorSpace(vs₁, vs₂.contents...)
@inline ⊕(vs₁::ZippedNamedVectorSpace, vs₂::SimpleNamedVectorSpace) = ZippedNamedVectorSpace(vs₁.contents..., vs₂)
@inline ⊕(vs₁::ZippedNamedVectorSpace, vs₂::ZippedNamedVectorSpace) = ZippedNamedVectorSpace(vs₁.contents..., vs₂.contents...)

"""
    ⊗(vs₁::SimpleNamedVectorSpace, vs₂::SimpleNamedVectorSpace) -> DirectProductedNamedVectorSpace
    ⊗(vs₁::SimpleNamedVectorSpace, vs₂::DirectProductedNamedVectorSpace) -> DirectProductedNamedVectorSpace
    ⊗(vs₁::DirectProductedNamedVectorSpace, vs₂::SimpleNamedVectorSpace) -> DirectProductedNamedVectorSpace
    ⊗(vs₁::DirectProductedNamedVectorSpace, vs₂::DirectProductedNamedVectorSpace) -> DirectProductedNamedVectorSpace

The direct product of named vector spaces.
"""
@inline ⊗(vs₁::SimpleNamedVectorSpace, vs₂::SimpleNamedVectorSpace) = DirectProductedNamedVectorSpace(vs₁, vs₂)
@inline ⊗(vs₁::SimpleNamedVectorSpace, vs₂::DirectProductedNamedVectorSpace) = DirectProductedNamedVectorSpace(vs₁, vs₂.contents...)
@inline ⊗(vs₁::DirectProductedNamedVectorSpace, vs₂::SimpleNamedVectorSpace) = DirectProductedNamedVectorSpace(vs₁.contents..., vs₂)
@inline ⊗(vs₁::DirectProductedNamedVectorSpace, vs₂::DirectProductedNamedVectorSpace) = DirectProductedNamedVectorSpace(vs₁.contents..., vs₂.contents...)

end # module


text:
# Constraint templates relating to network components or quantities not introduced by FlexPlan

## Power balance

"Power balance including candidate storage"
function constraint_power_balance_acne_dcne_strg(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    bus = _PM.ref(pm, nw, :bus, i)
    bus_arcs = _PM.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = _PM.ref(pm, nw, :ne_bus_arcs, i)
    bus_arcs_dc = _PM.ref(pm, nw, :bus_arcs_dc, i)
    bus_gens = _PM.ref(pm, nw, :bus_gens, i)
    bus_convs_ac = _PM.ref(pm, nw, :bus_convs_ac, i)
    bus_convs_ac_ne = _PM.ref(pm, nw, :bus_convs_ac_ne, i)
    bus_loads = _PM.ref(pm, nw, :bus_loads, i)
    bus_shunts = _PM.ref(pm, nw, :bus_shunts, i)
    bus_storage = _PM.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = _PM.ref(pm, nw, :bus_storage_ne, i)

    pd = Dict(k => _PM.ref(pm, nw, :load, k, "pd") for k in bus_loads)
    qd = Dict(k => _PM.ref(pm, nw, :load, k, "qd") for k in bus_loads)

    gs = Dict(k => _PM.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => _PM.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)
    constraint_power_balance_acne_dcne_strg(pm, nw, i, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
end

"Power balance (without DC equipment) including candidate storage"
function constraint_power_balance_acne_strg(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    bus = _PM.ref(pm, nw, :bus, i)
    bus_arcs = _PM.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = _PM.ref(pm, nw, :ne_bus_arcs, i)
    bus_gens = _PM.ref(pm, nw, :bus_gens, i)
    bus_loads = _PM.ref(pm, nw, :bus_loads, i)
    bus_shunts = _PM.ref(pm, nw, :bus_shunts, i)
    bus_storage = _PM.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = _PM.ref(pm, nw, :bus_storage_ne, i)

    pd = Dict(k => _PM.ref(pm, nw, :load, k, "pd") for k in bus_loads)
    qd = Dict(k => _PM.ref(pm, nw, :load, k, "qd") for k in bus_loads)

    gs = Dict(k => _PM.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => _PM.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)
    constraint_power_balance_acne_strg(pm, nw, i, bus_arcs, bus_arcs_ne, bus_gens, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
end

"Power balance including candidate storage & flexible demand"
function constraint_power_balance_acne_dcne_flex(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    bus = _PM.ref(pm, nw, :bus, i)
    bus_arcs = _PM.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = _PM.ref(pm, nw, :ne_bus_arcs, i)
    bus_arcs_dc = _PM.ref(pm, nw, :bus_arcs_dc, i)
    bus_gens = _PM.ref(pm, nw, :bus_gens, i)
    bus_convs_ac = _PM.ref(pm, nw, :bus_convs_ac, i)
    bus_convs_ac_ne = _PM.ref(pm, nw, :bus_convs_ac_ne, i)
    bus_loads = _PM.ref(pm, nw, :bus_loads, i)
    bus_shunts = _PM.ref(pm, nw, :bus_shunts, i)
    bus_storage = _PM.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = _PM.ref(pm, nw, :bus_storage_ne, i)

    gs = Dict(k => _PM.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => _PM.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)
    constraint_power_balance_acne_dcne_flex(pm, nw, i, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, gs, bs)
end

"Power balance (without DC equipment) including candidate storage & flexible demand"
function constraint_power_balance_acne_flex(pm::_PM.AbstractPowerModel, i::Int; nw::Int=_PM.nw_id_default)
    bus = _PM.ref(pm, nw, :bus, i)
    bus_arcs = _PM.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = _PM.ref(pm, nw, :ne_bus_arcs, i)
    bus_gens = _PM.ref(pm, nw, :bus_gens, i)
    bus_loads = _PM.ref(pm, nw, :bus_loads, i)
    bus_shunts = _PM.ref(pm, nw, :bus_shunts, i)
    bus_storage = _PM.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = _PM.ref(pm, nw, :bus_storage_ne, i)

    gs = Dict(k => _PM.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => _PM.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)
    constraint_power_balance_acne_flex(pm, nw, i, bus_arcs, bus_arcs_ne, bus_gens, bus_loads, bus_shunts, bus_storage, bus_storage_ne, gs, bs)
end

## AC candidate branches

"Activate a candidate AC branch depending on the investment decisions in the candidate's horizon."
function constraint_ne_branch_activation(pm::_PM.AbstractPowerModel, i::Int, prev_nws::Vector{Int}, nw::Int)
    investment_horizon = [nw]
    lifetime = _PM.ref(pm, nw, :ne_branch, i, "lifetime")
    for n in Iterators.reverse(prev_nws[max(end-lifetime+2,1):end])
        i in _PM.ids(pm, n, :ne_branch) ? push!(investment_horizon, n) : break
    end
    constraint_ne_branch_activation(pm, nw, i, investment_horizon)
end

## DC candidate branches

"Activate a candidate DC branch depending on the investment decisions in the candidate's horizon."
function constraint_ne_branchdc_activation(pm::_PM.AbstractPowerModel, i::Int, prev_nws::Vector{Int}, nw::Int)
    investment_horizon = [nw]
    lifetime = _PM.ref(pm, nw, :branchdc_ne, i, "lifetime")
    for n in Iterators.reverse(prev_nws[max(end-lifetime+2,1):end])
        i in _PM.ids(pm, n, :branchdc_ne) ? push!(investment_horizon, n) : break
    end
    constraint_ne_branchdc_activation(pm, nw, i, investment_horizon)
end

## Candidate converters

"Activate a candidate AC/DC converter depending on the investment decisions in the candidate's horizon."
function constraint_ne_converter_activation(pm::_PM.AbstractPowerModel, i::Int, prev_nws::Vector{Int}, nw::Int)
    investment_horizon = [nw]
    lifetime = _PM.ref(pm, nw, :convdc_ne, i, "lifetime")
    for n in Iterators.reverse(prev_nws[max(end-lifetime+2,1):end])
        i in _PM.ids(pm, n, :convdc_ne) ? push!(investment_horizon, n) : break
    end
    constraint_ne_converter_activation(pm, nw, i, investment_horizon)
end


text:
n = 10
m = n*(n-1)/2
c = ones(Int, n)
g = CompleteGraph(n)
@test  modularity(g, c) == 0
#
g = Graph(n)
@test modularity(g, c) == 0


text:
using MaterialModels
using Test
using Rotations
using StaticArrays
using JLD2

include("test_utils.jl")
include("test_linear_elastic.jl")
include("test_plastic.jl")
include("test_crystal_visco_plastic.jl")
include("test_crystal_visco_plastic_red.jl")
include("test_wrappers.jl")
include("test_neohook.jl")
include("test_yeoh.jl")
include("test_stvenant.jl")
include("test_straintraits.jl")


text:
"""
`X = DataUtils.expand_poly(x; degree = 5)`

Performs a polynomial basis expansion of the given `degree`
for the vector `x`.
The return value `X` is a matrix of size `(degree, length(x))`.

Note: all the features of `X` are centered and rescaled.
"""
function expand_poly(x::AbstractVector{T}; degree::Int = 5) where T<:Number
    n = length(x)
    x_vec = convert(Vector{Float64}, copy(collect(x)))
    X = zeros(Float64, (degree, n))
    @inbounds for i = 1:n
        for d = 1:degree
            X[d, i] += x_vec[i]^(d)
        end
    end
    X
end


text:
# Struct representing a point abstraction.
abstract type Element end

include("Heaps.jl")

# mutable struct Candidate{S<:String, C<:Array{Float64,1}, W<:Float64} <: Element
#     id::S
#     coords::C
#     weight::W
# end
# Candidate() = Candidate("", Array{Float64,1}(), 0.)
# Candidate(s::String) = Candidate(s, Array{Float64,1}(), 0.)

# Struct that represents the point, complete with heap of nearby points.
mutable struct Point <: Element
    # Point id
    id::String
    # Point coordinates
    coords::Vector{Float64}
    # Heap of nearest neighbors
    nn::NN
    # If true, the point is marked as active
    active::Bool
end

Point() = Point("", Vector{Float64}(), NN(), true)
Point(k::Int64) = Point("", Vector{Float64}(), NN(k), true)
Point(s::String) = Point(s, Vector{Float64}(), NN(), true)
Point(s::String, coords::Vector{Float64}) = Point(s, coords, NN(), true)
Point(s::String, coords::Vector{Float64}, k::Int64) = Point(s, coords, NN(k), true)

isActive(el::Element) = return el.active
getDimension(this::Element) = return length(this.coords)

# toCandidate(p::Point) = Candidate(p.id, p.coords, p.weight)

# Element compare by weight
Base.isless(a::Element,b::Element) = isless(a.nn.weight,b.nn.weight)

function Base.show(io::IO, obj::Element)
    println(string(obj.id), " | ", string(obj.coords), " | ", string(obj.nn.weight))
end


text:
function full_factorial(arity, disallow)
    param_cnt = length(arity)
    test_cnt = prod(arity)
    test_set = zeros(Int, param_cnt, test_cnt)
    state = ones(Int, param_cnt)
    keep_cnt = 0
    for gen_idx in 1:test_cnt
        if !disallow(state)
            keep_cnt += 1
            test_set[:, keep_cnt] = state
        end
        next_multiplicative!(state, arity)
    end
    # Test the invariant that the state rolls over to all ones.
    @assert state == ones(Int, param_cnt)
    test_set[:, 1:keep_cnt]
end


text:
module SleepstudyTest

println()
@info "Sleepstudy example"

using GLMCopula, RDatasets, Test, LinearAlgebra

# Dataframe with columns: Reaction (Float64), Days (Int32), Subject (Categorical)
# 180 rows
sleepstudy = dataset("lme4", "sleepstudy")
groups = unique(sleepstudy[!, :Subject])
n, p, m = length(groups), 2, 1
gcs = Vector{GaussianCopulaLMMObs{Float64}}(undef, n)
for (i, grp) in enumerate(groups)
    gidx = sleepstudy[!, :Subject] .== grp
    ni = count(gidx)
    yi = Float64.(sleepstudy[gidx, :Reaction])
    Ui = Float64.(sleepstudy[gidx, :Days])
    Xi = [ones(ni, 1) Ui]
    Zi = Xi
    gcs[i] = GaussianCopulaLMMObs(yi, Xi, Zi)
end
gcm = GaussianCopulaLMMModel(gcs);

# initialize β and τ from least squares solution
@info "Initial point:"
init_β!(gcm)
@show gcm.β
@show gcm.τ[1]
gcm.Σ .= diagm(ones(2))
@show gcm.Σ
@show loglikelihood!(gcm, true, false)
@show gcm.∇β
@show gcm.∇τ
@show gcm.∇Σ

# # initialize from LMM estimate
# @info "Initial point:"
# gcm.β .= [251.405, 10.4673]
# gcm.τ[1] = inv(654.94145)
# gcm.Σ .= [565.51069 10.87591; 10.87591 32.68212]
# @show loglikelihood!(gcm, true, false)
# @show gcm.∇β
# @show gcm.∇τ
# @show gcm.∇Σ

# fit model using NLP
@info "MLE:"
# @time GLMCopula.fit!(gcm, NLopt.NLoptSolver(algorithm=:LD_SLSQP, maxeval=4000))
@time GLMCopula.fit!(gcm, IpoptSolver(print_level=5))
@show gcm.β
@show gcm.τ
@show gcm.Σ
@show loglikelihood!(gcm, true, false)
# @test loglikelihood!(gcm, true, false) ≈ -163.35545251
@show gcm.∇β
@show gcm.∇τ
@show gcm.∇Σ

end


text:
export Sequence
export next!

mutable struct Sequence
  id::Int32

  Sequence() = new(0)
end

function next!(seq::Sequence)
  seq.id += 1
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
using DataStructures

###############################################################################
# dep graph when MassActionJump uses pairs to represent (species,stoich)
###############################################################################

# map from species to Set of reactions depending on that species
function spec_to_dep_rxs_map(numspec, ma_jumps::MassActionJump)

    numrxs = get_num_majumps(ma_jumps)

    # map from a species to reactions that depend on it
    spec_to_dep_rxs = [Set{Int}() for n = 1:numspec]
    for rx in 1:numrxs
        for (spec,stoch) in ma_jumps.reactant_stoch[rx]
            push!(spec_to_dep_rxs[spec], rx)
        end
    end

    spec_to_dep_rxs
end

# dependency graph is a map from a reaction to a vector of reactions
# that should depend on species it changes
function make_dependency_graph(numspec, ma_jumps::MassActionJump)

    numrxs          = get_num_majumps(ma_jumps)
    spec_to_dep_rxs = spec_to_dep_rxs_map(numspec, ma_jumps)

    # create map from rx to reactions depending on it
    dep_sets = [SortedSet{Int}() for n = 1:numrxs]
    for rx in 1:numrxs

        # rx changes spec, hence rxs depending on spec depend on rx
        for (spec,stoch) in ma_jumps.net_stoch[rx]
            for dependent_rx in spec_to_dep_rxs[spec]
                push!(dep_sets[rx], dependent_rx)
            end
        end
    end

    # convert to Vectors of Vectors
    dep_graph = Vector{Vector{Int}}(numrxs)
    for rx = 1:numrxs
        dep_graph[rx] = [dep for dep in dep_sets[rx]]
    end

    dep_graph
end

# # STATIC MATRICES:

# # map from species to Set of reactions depending on that species
# function spec_to_dep_rxs_map(numspec, ma_jumps::MassActionJump)

#     numrxs = get_num_majumps(ma_jumps)

#     # map from a species to reactions that depend on it
#     spec_to_dep_rxs = [Set{Int}() for n = 1:numspec]
#     for rx in 1:numrxs
#         for (spec,stoch) in ma_jumps.majdat[rx].reactant_stoch
#             if spec > 0
#                 push!(spec_to_dep_rxs[spec], rx)
#             end
#         end
#     end

#     spec_to_dep_rxs
# end

# # dependency graph is a map from a reaction to a vector of reactions
# # that should depend on species it changes
# function make_dependency_graph(numspec, ma_jumps::MassActionJump)

#     numrxs          = get_num_majumps(ma_jumps)
#     spec_to_dep_rxs = spec_to_dep_rxs_map(numspec, ma_jumps)

#     # create map from rx to reactions depending on it
#     dep_sets = [SortedSet{Int}() for n = 1:numrxs]
#     for rx in 1:numrxs

#         # rx changes spec, hence rxs depending on spec depend on rx
#         for (spec,stoch) in ma_jumps.majdat[rx].net_stoch
#             if spec > 0
#                 for dependent_rx in spec_to_dep_rxs[spec]
#                     push!(dep_sets[rx], dependent_rx)
#                 end
#             end
#         end
#     end

#     # convert to Vectors of Vectors
#     dep_graph = Vector{Vector{Int}}(numrxs)
#     for rx = 1:numrxs
#         dep_graph[rx] = [dep for dep in dep_sets[rx]]
#     end

#     dep_graph
# end

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
module Symbolic

using OffsetArrays

using ..TypeSynonyms

struct VarLinForm{T}
    form ::Dict{Int, T}
end

VarLinForm{T}() where T = VarLinForm(Dict{Int, T}())

Base.getindex(vlf::VarLinForm{T}, key::Int) where T = get(vlf.form, key, zero(T))

Base.setindex!(vlf::VarLinForm{T}, val::T, key::Int) where T = setindex!(vlf.form, val, key)

variables(vlf::VarLinForm) = keys(vlf.form)

function add!(dest::VarLinForm{T}, source::VarLinForm{T}, factor=1) where {T}
    for var in variables(source)
        dest[var] += source[var] * factor
    end
end

function add_conjugated!(dest::VarLinForm{T}, source::VarLinForm{T}, factor=1) where {T}
    for var in variables(source)
        dest[var] += conj(source[var]) * factor
    end
end

function mul!(form::VarLinForm, factor)
    for var in variables(form)
        form[var] *= factor
    end
end

###############################################################################

struct PolynomialForm{N <: Number}
    terms ::OffsetVector{ VarLinForm{N} }
end

function PolynomialForm{N}(index_range::UnitRange) where {N <: Number}
    PolynomialForm(
        OffsetVector{VarLinForm{N}}([VarLinForm{N}() for i in index_range], index_range)
    )
end

Base.getindex(poly::PolynomialForm, key::Int) = poly.terms[key]

firstindex(poly::PolynomialForm) = firstindex(poly.terms)
lastindex(poly::PolynomialForm) = lastindex(poly.terms)
eachpower(poly::PolynomialForm) = eachindex(poly.terms)

function add!(dest::PolynomialForm{N}, source::PolynomialForm{N}, factor=1) where {N}
    for power in eachpower(source)
        add!(dest[power], source[power],  factor)
    end
end

function mul!(poly::PolynomialForm, factor)
    for power in eachpower(poly)
        mul!(poly[power], factor)
    end
end

"""
    mul_by_power!(poly::PolynomialForm, factor)

Transform z^n to (z^n * factor^n) for each polynomial term.
Powers indices must encompass 0.
"""
function mul_by_power!(poly::PolynomialForm, factor)
    bottom = firstindex(poly)
    top  = lastindex(poly)
    if (bottom <= 0) && (top >= 0)
        min_bound = minimum(-bottom, top)
        factor_pow = factor

        # Transform central symmetrical part of polynomial
        for power in 1 : min_bound
            mul!(poly[power], factor_pow)
            mul!(poly[-power], 1/factor_pow)
            factor_pow *= factor
        end

        # Transform margins if any
        for power in min_bound+1 : top
            mul!(poly[power], factor_pow)
            factor_pow *= factor
        end
        for power in min_bound+1 : -bottom
            mul!(poly[-power], 1/factor_pow)
            factor_pow *= factor
        end
    else
        error("Indices don't encompass 0. This situation doesn't treated here")
    end
end

"""
    conjugate(poly::PolynomialForm)

Create conjugated form of a given polynomial form.
"""
function conjugate(poly::PolynomialForm{N}) where {N}
    res = PolynomialForm{N}(-eachpower(poly))
    for i in eachpower(res)
        add_conjugated!(res[i], poly[-i])
    end
    return res
end

"""
    z_conj_diff(poly::PolynomialForm)

For given polynomial form ϕ(z) creates z bar Φ(z).
"""
function z_conj_diff(poly::PolynomialForm{N}) where {N}
    bottom = firstindex(poly)
    top    = lastindex(poly)
    res = PolynomialForm{N}(-top+2 : -bottom+2)

    for n in eachpower(poly)
        add_conjugated!(res[-n+2], poly[n], n)
    end
    return res
end

function matrix_form(poly::PolynomialForm{N}) where {N}
    bottom = firstindex(poly)
    top    = lastindex(poly)
    min_var_n, max_var_n = extrema(variables(poly[bottom]))
    for n in bottom+1 : top
        min_n, max_n = extrema(variables(poly[n]))
        (min_n < min_var_n) && (min_var_n = min_n)
        (max_n < max_var_n) && (max_var_n = max_n)
    end

    matrix = OffsetArray{N, 2}(undef, bottom:top, min_var_n:max_var_n)
    fill!(matrix, zero(N))

    for pow in eachpower(poly)
        for var in variables(poly[pow])
            matrix[pow, var] = poly[pow][var]
        end
    end

    return matrix
end

#######################################

struct NormedPolynomial{N <: Number, FL <: AbstractFloat}
    poly ::PolynomialForm{N}
    r_norm ::FL
end

function add!(dest::NormedPolynomial, source::NormedPolynomial, factor=1) end

function empty!(normed_poly::NormedPolynomial) end

function conjugate(normed_poly::NormedPolynomial, r_contour) end

function re_conjugate!(normed_poly::NormedPolynomial, r_old, r_new) end

function z_conj_diff(normed_poly::NormedPolynomial, r_contour) end

function matrix_form(normed_poly::NormedPolynomial) end

#=

###############################################################################

abstract type GeneralFunction end

struct FunctionalTerm{F <: GeneralFunction}
    func   ::F
    factor ::Variable
end

"""
  (z/r)^n

power -  n,

Norming factor ''r'' is stored elsewhere.
"""
mutable struct ZPowN <: GeneralFunction
    power  ::Int64
end

function differentiate!(f::ZPowN, r_norm::Float64)
    # [(z-a)^n/r^n]' = (n/r) z^{n-1}/r^{n-1}
    multiplier = f.power / r_norm
    f.power -= 1
    return multiplier
end

###############################################################################

struct Polynomial <: GeneralFunction
    terms ::Dict{Int64, ComplexF64}
end

###############################################################################

"""
  r^{n+2} ℘^{(n)}(z)/(n+1)!

deriv -  n.

Norming factor ''r'' is stored elsewhere.
"""
mutable struct WpDerivative <: GeneralFunction
    deriv::Int64
end

function differentiate!(f::WpDerivative, r_norm::Float64)
    # [r_k^{n+2} ℘^{(n)}(z-ζ_k)/(n+1)!]' =
    # [(n+2) / r_k] *
    # r_k^{n+3} ℘^{(n+1)}(z-ζ_k)/(n+2)!
    multiplier = (f.deriv+2) / r_norm
    f.deriv += 1
    return multiplier
end

###############################################################################

"""
  r^{n+2} Q^{(n)}(z)/(n+1)!

deriv -  n.

Norming factor ''r'' is stored elsewhere.
"""
mutable struct QSpecial <: GeneralFunction
    deriv::Int64
end

function differentiate!(f::QSpecial, r_norm::Float64)
    # [r_k^{n+2} Q^{(n)}(z-ζ_k)/(n+1)!]' =
    # [(n+2) / r_k] *
    # r_k^{n+3} Q^{(n+1)}(z-ζ_k)/(n+2)!
    multiplier = (f.deriv+2) / r_norm
    f.deriv += 1
    return multiplier
end

###############################################################################

mutable struct Variable
    number      ::Int64
    is_computed ::Bool
    value       ::Float64
end

mutable struct Term{F <: GeneralFunction}
    variable    ::Variable
    coefficient ::ComplexF64
    functional  ::F
end

function differentiate!(term::Term)
    multiplier = differentiate!(term.functional)
    term.coefficient *= multiplier
    nothing
end

###############################################################################

abstract type FunctionalForm end

struct LinearForm{F<:GeneralFunction} <: FunctionalForm
    terms          ::Vector{Term{F}}
    shift_rational ::RationalComplex
    shift_raw      ::ComplexF64
    norm_factor    ::Float64
    conjugated     ::Bool
    mul_by_z       ::Bool
end

function conjugate!(form::LinearForm)
    for term in form.terms
        term.coefficient = conj(term.coefficient)
    end
    form.conjugated = !form.conjugated
    nothing
end

function mul!(form::LinearForm, factor::Float64)
    for term in form.terms
        term.coefficient *= factor
    end
    nothing
end

struct SolutionForm <:FunctionalForm
    terms ::Vector{LinearForm}
end

function differentiate!(form <:FunctionalForm)
    map(differentiate!, form.terms)
    nothing
end

function conjugate!(form <:FunctionalForm)
    map(conjugate!, form.terms)
    nothing
end

function mul!(form <:FunctionalForm, factor::Float64)
    f!(t) = mul!(t, factor)
    map(f!, form.terms)
    nothing
end
=#
end # module Symbolic


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2021 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# 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.
# ----------------------------------------------------------------------------------- #

# calculate the u-variables -
function calculate_transcription_control_array(t::Float64, x::Array{Float64,1},
    parameter_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    number_of_transcription_processes = parameter_dictionary["number_of_transcription_processes"]
    model_parameter_array = parameter_dictionary["model_parameter_array"]
    model_parameter_index_map = parameter_dictionary["model_parameter_symbol_index_map"]
    u_array = Array{Float64,1}()

    # local helper functions -
    f(x,K,n) = (x^n)/(K^n+x^n)
    u(A,R) = A/(1+A+R)

    # alias the state vector -
    gene_venus = x[1]
	mRNA_venus = x[2]
	P_venus = x[3]

    # alias the system species -
    system_array = parameter_dictionary["system_concentration_array"]
    RNAP = system_array[1]
	RIBOSOME = system_array[2]
	P_σ70 = system_array[3]

    # alias the system parameters -
    W_gene_venus = model_parameter_array[model_parameter_index_map[:W_gene_venus]]
	W_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:W_gene_venus_P_σ70]]
	K_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:K_gene_venus_P_σ70]]
	n_gene_venus_P_σ70 = model_parameter_array[model_parameter_index_map[:n_gene_venus_P_σ70]]
	K_gene_venus = model_parameter_array[model_parameter_index_map[:K_gene_venus]]
	𝛕_gene_venus = model_parameter_array[model_parameter_index_map[:𝛕_gene_venus]]
	K_P_venus = model_parameter_array[model_parameter_index_map[:K_P_venus]]
	𝛕_P_venus = model_parameter_array[model_parameter_index_map[:𝛕_P_venus]]
	𝛳_mRNA_venus = model_parameter_array[model_parameter_index_map[:𝛳_mRNA_venus]]
	𝛳_P_venus = model_parameter_array[model_parameter_index_map[:𝛳_P_venus]]

    # == CONTROL LOGIC BELOW ================================================================= #

	# - gene_venus ---------------------------------------------------------------------------------------------
	# gene_venus activation -
	W = [
			W_gene_venus	;
			W_gene_venus_P_σ70	;
	]
	gene_venus_activator_array = Array{Float64,1}()
	push!(gene_venus_activator_array, 1.0)
	P_σ70_RNAP = RNAP*f(P_σ70,K_gene_venus_P_σ70,n_gene_venus_P_σ70)
	push!(gene_venus_activator_array, P_σ70_RNAP)
	A = sum(W.*gene_venus_activator_array)

	# gene_venus repression -
	R = 0.0
	push!(u_array, u(A,R))

    # == CONTROL LOGIC ABOVE ================================================================= #

    # return -
    return u_array
end

# calculate the w-variables -
function calculate_translation_control_array(t::Float64, x::Array{Float64,1},
    parameter_dictionary::Dict{String,Any})::Array{Float64,1}

    # initialize -
    number_of_translation_processes = parameter_dictionary["number_of_translation_processes"]

    # encode: encode the w-variables -
    # translation decay: translation capacity decays over time -
	correction_term = (x[4]*(100.0^-1))
	w = ones(number_of_translation_processes).*(correction_term)

    # return -
    return w
end


text:
module ExportPublicTest

    using Test
    include("Bar.jl")
    using .Bar
    include("SimpleMath.jl")
    using .SimpleMath
    include("DemoModule.jl")
    using .DemoModule

    @testset "Bar.jl" begin
        @test foo() == 1
        @test bar() == 2
        @test barBar() == 3
        @test fooFoo() == 4
        @test_throws UndefVarError _secret_implementation_detail()
    end

    @testset "SimpleMath.jl" begin
        @test add_two(5, 5) == 10
        @test add_two(10, 10) == 20
        @test_throws UndefVarError _implementation_detail()
        @test add_squared(1, 1) == 2
        @test add_squared(2, 2) == 8
        @test add_squared(4, 4) == 32
        @test_throws UndefVarError _square_it()
        @test my_pi == 22/7
        @test_throws UndefVarError _secret_pi
    end

    @testset "DemoModule.jl" begin
        @test add(5, 5) == 10
        @test add(6, 6) == 12
        @test_throws UndefVarError _private_helper()
    end
end


text:
function asArray(f::Function, A::GBVecOrMat{T}; dropzeros=false, freeunpacked=false) where {T}
    if gbget(A, SPARSITY_STATUS) != GBDENSE
        X = similar(A)
        if X isa GBVector
            X[:] = zero(T)
        else
            X[:,:] = zero(T)
        end
        #I don't like this, it defeats the purpose of this method, which is to make no copies.
        # But somehow maintaining the input A in its original form is key to the to_vec implementation
        # for ChainRules. Temporarily it's fine, it's no worse than it originally was.
        # TODO: fix this issue with the ChainRules code.
        A = eadd(X, A)
    end
    array = _unpackdensematrix!(A)
    result = try
        f(array, A)
    finally
        if freeunpacked
            ccall(:jl_free, Cvoid, (Ptr{T},), pointer(array))
        else
            _packdensematrix!(A, array)
            if dropzeros
                select!(nonzeros, A)
            end
        end
    end
    return result
end

function asCSCVectors(f::Function, A::GBMatrix{T}; freeunpacked=false) where {T}
    colptr, rowidx, values =  _unpackcscmatrix!(A)
    result = try
        f(colptr, rowidx, values, A)
    finally
        println("I'm finally")
        if freeunpacked
            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(colptr))
            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(rowidx))
            ccall(:jl_free, Cvoid, (Ptr{T},), pointer(values))
        else
            println("I repacked!")
            _packcscmatrix!(A, colptr, rowidx, values)
        end
    end
    return result
end

function asCSRVectors(f::Function, A::GBMatrix{T}; freeunpacked=false) where {T}
    rowptr, colidx, values =  _unpackcsrmatrix!(A)
    result = try
        f(rowptr, colidx, values, A)
    finally
        if freeunpacked
            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(rowptr))
            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(colidx))
            ccall(:jl_free, Cvoid, (Ptr{T},), pointer(values))
        else
            _packcsrmatrix!(A, rowptr, colidx, values)
        end
    end
    return result
end

function asSparseMatrixCSC(f::Function, A::GBMatrix{T}; freeunpacked=false) where {T}
    colptr, rowidx, values =  _unpackcscmatrix!(A)
    array = SparseMatrixCSC{T, libgb.GrB_Index}(size(A, 1), size(A, 2), colptr, rowidx, values)
    result = try
        f(array, A)
    finally
        if freeunpacked
            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(colptr))
            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(rowidx))
            ccall(:jl_free, Cvoid, (Ptr{T},), pointer(values))
        else
            _packcscmatrix!(A, colptr, rowidx, values)
        end
    end
    return result
end

function asSparseVector(f::Function, A::GBVector{T}; freeunpacked=false) where {T}
    colptr, rowidx, values =  _unpackcscmatrix!(A)
    vector = SparseVector{T, libgb.GrB_Index}(size(A, 1), rowidx, values)
    result = try
        f(vector, A)
    finally
        if freeunpacked
            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(colptr))
            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(rowidx))
            ccall(:jl_free, Cvoid, (Ptr{T},), pointer(values))
        else
            _packcscmatrix!(A, colptr, rowidx, values)
        end
    end
    return result
end

function Base.Matrix(A::GBMatrix)
    return asArray(A) do arr, _
        return copy(arr)
    end
end

function Matrix!(A::GBMatrix)
    return asArray(A; freeunpacked=true) do arr, _
        return copy(arr)
    end
end

function Base.Vector(v::GBVector)
    return asArray(v) do vec, _
        return copy(vec)
    end
end

function Vector!(v::GBVector)
    return asArray(v; freeunpacked=true) do vec, _
        return copy(vec)
    end
end

function SparseArrays.SparseMatrixCSC(A::GBMatrix)
    return asArray(A) do arr, _
        return copy(arr)
    end
end

function SparseArrays.SparseVector(v::GBVector)
    return asArray(v) do arr, _
        return copy(arr)
    end
end


text:
export StructTree

using AbstractTrees
using Random
using ProgressMeter

const AT = AbstractTrees

struct StructTree{X}
    x::X
end

AT.children(t::StructTree{X}) where {X} =
    Tuple(f => StructTree(getfield(t.x, f)) for f in fieldnames(X))
AT.children(
    t::StructTree{T},
) where {
    T<:Union{AbstractArray,AbstractDict,MersenneTwister,ProgressMeter.Progress,Function},
} = ()
AT.children(t::Pair{Symbol,<:StructTree}) = children(last(t))
AT.children(t::StructTree{UnionAll}) = ()

AT.printnode(io::IO, t::StructTree{<:Union{Number,Symbol}}) = print(io, t.x)
AT.printnode(io::IO, t::StructTree{UnionAll}) = print(io, t.x)
AT.printnode(io::IO, t::StructTree{T}) where {T} = print(io, T.name)
AT.printnode(io::IO, t::StructTree{<:AbstractArray}) where {T} = summary(io, t.x)

function AT.printnode(io::IO, t::StructTree{String})
    s = t.x
    i = findfirst('\n', s)
    if isnothing(i)
        if length(s) > 79
            print(io, "\"$(s[1:79])...\"")
        else
            print(io, "\"$s\"")
        end
    else
        if i > 79
            print(io, "\"$(s[1:79])...\"")
        else
            print(io, "\"$(s[1:i-1])...\"")
        end
    end
end

function AT.printnode(io::IO, t::Pair{Symbol,<:StructTree})
    print(io, first(t), " => ")
    AT.printnode(io, last(t))
end

function AT.printnode(io::IO, t::Pair{Symbol,<:StructTree{<:Tuple}})
    print(io, first(t))
end


text:
include("plotHelpers.jl")
Egrid, PAgrid = logrange(10,1000,21), 6:4:90

# Plots for 4/29/21
@time @load "results/210429_lat_storage.jld2" equatorial_fluxes_042921 themis_hilat_042921 themis_midlat_042921 themis_lolat_042921
elfin_measurements_042921, elfin_error_042921 = extract_idl_csv("042921_time.csv", "042921_prec.csv",
                                                "042921_precerror.csv", "042921_precerror_time.csv", "ebins.csv", # csvs containing ELFIN measurements
                                                DateTime(2021,4,29,3,14,45),DateTime(2021,4,29,3,15,0)); # time to sample from ELFIN measurements

x1 = themis_lolat_042921.precipitating_fluxes_mean.-themis_lolat_042921.precipitating_fluxes_minus .+ 1.
x2 = themis_lolat_042921.precipitating_fluxes_mean.+themis_lolat_042921.precipitating_fluxes_plus .+ 1.
y1 = themis_midlat_042921.precipitating_fluxes_mean.-themis_midlat_042921.precipitating_fluxes_minus .+ 1.
y2 = themis_midlat_042921.precipitating_fluxes_mean.+themis_midlat_042921.precipitating_fluxes_plus .+ 1.
z1 = themis_hilat_042921.precipitating_fluxes_mean.-themis_hilat_042921.precipitating_fluxes_minus .+ 1.
z2 = themis_hilat_042921.precipitating_fluxes_mean.+themis_hilat_042921.precipitating_fluxes_plus .+ 1.
# elfin_error_042921[8] += 2000.
e1 = elfin_measurements_042921[2].-elfin_error_042921
e2 = elfin_measurements_042921[2].+elfin_error_042921
sum_of_means = @. themis_hilat_042921.precipitating_fluxes_mean + themis_midlat_042921.precipitating_fluxes_mean + themis_lolat_042921.precipitating_fluxes_mean
error_sum_minus = @. sqrt(themis_hilat_042921.precipitating_fluxes_minus^2 + themis_midlat_042921.precipitating_fluxes_minus^2 + themis_lolat_042921.precipitating_fluxes_minus^2)
error_sum_plus = @. sqrt(themis_hilat_042921.precipitating_fluxes_plus^2 + themis_midlat_042921.precipitating_fluxes_plus^2 + themis_lolat_042921.precipitating_fluxes_plus^2)
s1 = sum_of_means.-error_sum_minus .+ 1.
s2 = sum_of_means.+error_sum_plus .+ 1.

plot(Egrid, equatorial_fluxes_042921,label = L"\mathrm{equatorial\ flux\ (THEMIS\ A\ measured)}", color = bipride_pink, linewidth=2, markershape=:circle);
plot!(Egrid,themis_lolat_042921.precipitating_fluxes_mean, fillrange=x1, fillalpha = 0.2, color = c4, label=false)
plot!(Egrid,themis_lolat_042921.precipitating_fluxes_mean, fillrange=x2, fillalpha = 0.2, color = c4, label=false)
plot!(Egrid, themis_lolat_042921.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\lambda_2 = 20\degree\mathrm{,\ occurence\ rate} = 30\%)",
        # yerror=(themis_lolat_042921.precipitating_fluxes_minus, themis_lolat_042921.precipitating_fluxes_plus),
        color = c4, alpha = .7, markerstrokecolor = c4, markerstrokewidth = 2, markershape=:o, markersize = 2);
plot!(Egrid,themis_midlat_042921.precipitating_fluxes_mean, fillrange=y1, fillalpha = 0.2, color = c3, label=false)
plot!(Egrid,themis_midlat_042921.precipitating_fluxes_mean, fillrange=y2, fillalpha = 0.2, color = c3, label=false)
plot!(Egrid, themis_midlat_042921.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\lambda_2 = 30\degree\mathrm{,\ occurence\ rate} = 4\%)",
        # yerror=(themis_midlat_042921.precipitating_fluxes_minus, themis_midlat_042921.precipitating_fluxes_plus),
        color = c3, alpha = .7, markerstrokecolor = c3, markerstrokewidth = 2, markershape=:o, markersize = 2);
plot!(Egrid,themis_hilat_042921.precipitating_fluxes_mean, fillrange=z1, fillalpha = 0.2, color = bipride_orange, label=false)
plot!(Egrid,themis_hilat_042921.precipitating_fluxes_mean, fillrange=z2, fillalpha = 0.2, color = bipride_orange, label=false)
plot!(Egrid, themis_hilat_042921.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\lambda_2 = 40\degree\mathrm{,\ occurence\ rate} = 2\%)",
        # yerror=(themis_hilat_042921.precipitating_fluxes_minus, themis_hilat_042921.precipitating_fluxes_plus),
        color = bipride_orange, alpha = .7, markerstrokecolor = bipride_orange, markerstrokewidth = 2, markershape=:o, markersize = 2);
plot!(Egrid,sum_of_means, fillrange=s1, fillalpha = 0.2, color = :black, label=false)
plot!(Egrid,sum_of_means, fillrange=s2, fillalpha = 0.2, color = :black, label=false)
plot!(Egrid, sum_of_means, label=L"\mathrm{precipitating\ flux\ (sum\ of\ simulations)}",
        # yerror=(error_sum_minus, error_sum_plus),
        color = :black, linewidth=3, markershape=:o, markersize = 3);
plot!(elfin_measurements_042921[1][1:end-5], elfin_measurements_042921[2][1:end-5], fillrange = e1[1:end-5], fillalpha = 0.2, color = c5, label=false)
plot!(elfin_measurements_042921[1][1:end-5], elfin_measurements_042921[2][1:end-5], fillrange = e2[1:end-5], fillalpha = 0.2, color = c5, label=false)
plot!(elfin_measurements_042921, yerror=elfin_error_042921.+10, label = L"\mathrm{precipitating\ flux\ (ELFIN\ B\ measured)}",
        color = c5, marker = stroke(3,c5), linewidth=4, markersize = 3);
plot!(ylim =(1e3,1e8), xlim=(50,800), yscale=:log10, legend=(0.45,0.7))
plot!(xlabel=L"\mathrm{Energy\ (keV)}", ylabel=L"\mathrm{Flux\ (1/cm^{2}/s/sr/MeV)}", title=L"\mathrm{Precipitating\ Flux\ Comparison\ of\ \delta\lambda_2\ on\ 4/29/21}",
xtickfontsize=12, ytickfontsize=12, xguidefontsize=16, yguidefontsize=16, legendfontsize=10, titlefontsize=16);
plot1 = plot!(dpi = 500,size=(800,450), margin=3mm, bottom_margin=4mm)
savefig(plot1, "images/042921_flux_lat_comparison.pdf")
savefig(plot1, "images/042921_flux_lat_comparison.png")

@time @load "results/210429_pkt_storage.jld2" equatorial_fluxes_042921 themis_1pkt_042921 themis_30pkt_042921 themis_100pkt_042921;

np1 = themis_1pkt_042921.precipitating_fluxes_mean.-themis_1pkt_042921.precipitating_fluxes_minus .+ 1.
np2 = themis_1pkt_042921.precipitating_fluxes_mean.+themis_1pkt_042921.precipitating_fluxes_plus .+ 1.
op1 = themis_30pkt_042921.precipitating_fluxes_mean.-themis_30pkt_042921.precipitating_fluxes_minus .+ 1.
op2 = themis_30pkt_042921.precipitating_fluxes_mean.+themis_30pkt_042921.precipitating_fluxes_plus .+ 1.
pp1 = themis_100pkt_042921.precipitating_fluxes_mean.-themis_100pkt_042921.precipitating_fluxes_minus .+ 1.
pp2 = themis_100pkt_042921.precipitating_fluxes_mean.+themis_100pkt_042921.precipitating_fluxes_plus .+ 1.

plot(Egrid, equatorial_fluxes_042921,label = L"\mathrm{equatorial\ flux\ (THEMIS\ A\ measured)}", color = bipride_pink, linewidth=2, markershape=:circle);
plot!(Egrid,themis_1pkt_042921.precipitating_fluxes_mean, fillrange=np1, fillalpha = 0.2, color = c4, label=false)
plot!(Egrid,themis_1pkt_042921.precipitating_fluxes_mean, fillrange=np2, fillalpha = 0.2, color = c4, label=false)
plot!(Egrid, themis_1pkt_042921.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\phi = 1)",
        # yerror=(themis_nopkt_042921.precipitating_fluxes_minus, themis_nopkt_042921.precipitating_fluxes_plus),
        color = c4, alpha = .7, markerstrokecolor = c4, markerstrokewidth = 2, markershape=:o, markersize = 2);
plot!(Egrid,themis_30pkt_042921.precipitating_fluxes_mean, fillrange=op1, fillalpha = 0.2, color = c3, label=false)
plot!(Egrid,themis_30pkt_042921.precipitating_fluxes_mean, fillrange=op2, fillalpha = 0.2, color = c3, label=false)
plot!(Egrid, themis_30pkt_042921.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\phi = 30)",
        # yerror=(themis_100pkt_042921.precipitating_fluxes_minus, themis_100pkt_042921.precipitating_fluxes_plus),
        color = c3, alpha = .7, markerstrokecolor = c3, markerstrokewidth = 2, markershape=:o, markersize = 2);
plot!(Egrid,themis_100pkt_042921.precipitating_fluxes_mean, fillrange=pp1, fillalpha = 0.2, color = bipride_orange, label=false)
plot!(Egrid,themis_100pkt_042921.precipitating_fluxes_mean, fillrange=pp2, fillalpha = 0.2, color = bipride_orange, label=false)
plot!(Egrid, themis_100pkt_042921.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\phi = 100)",
        # yerror=(themis_100pkt_042921.precipitating_fluxes_minus, themis_100pkt_042921.precipitating_fluxes_plus),
        color = bipride_orange, alpha = .7, markerstrokecolor = bipride_orange, markerstrokewidth = 2, markershape=:o, markersize = 2);
# plot!(elfin_measurements_042921[1][1:end-5], elfin_measurements_042921[2][1:end-5], fillrange = e1[1:end-5], fillalpha = 0.2, color = c5, label=false)
# plot!(elfin_measurements_042921[1][1:end-5], elfin_measurements_042921[2][1:end-5], fillrange = e2[1:end-5], fillalpha = 0.2, color = c5, label=false)
# plot!(elfin_measurements_042921, yerror=elfin_error_042921.+10, label = L"\mathrm{precipitating\ flux\ (ELFIN\ B\ measured)}",
#         color = c5, marker = stroke(3,c5), linewidth=4, markersize = 3);
plot!(ylim =(1e3,1e8), xlim=(50,800), yscale=:log10, legend=(0.45,0.7));
plot!(xlabel=L"\mathrm{Energy\ (keV)}", ylabel=L"\mathrm{Flux\ (1/cm^{2}/s/sr/MeV)}", title=L"\mathrm{Precipitating\ Flux\ Comparison\ of\ \delta\phi\ on\ 4/29/21}",
xtickfontsize=12, ytickfontsize=12, xguidefontsize=16, yguidefontsize=16, legendfontsize=10, titlefontsize=16);
plot2 = plot!(dpi = 500,size=(800,450), margin=3mm, bottom_margin=4mm)
# savefig(plot2, "images/042921_flux_pkt_comparison.pdf")
# savefig(plot2, "images/042921_flux_pkt_comparison.png")

include("plotHelpers.jl")
Egrid, PAgrid = logrange(10,1000,21), 6:4:90

# Plots for 9/22/20
@load "results/200922_lat_storage.jld2" equatorial_fluxes_092220 mms_lowerlat_092220 mms_lolat_092220 mms_midlat_092220 mms_hilat_092220;
@time elfin_measurements_092220, elfin_error_092220 = extract_idl_csv("092220_time.csv", "092220_prec.csv",
                                                "092220_precerror.csv", "092220_precerror_time.csv", "ebins.csv", # csvs containing ELFIN measurements
                                                DateTime(2020,9,22,9,16,38), DateTime(2020,9,22,9,16,45)) # time to sample from ELFIN measurements
# @time elfin_measurements_092220 = extract_idl_csv("092220_time.csv", "092220_prec.csv", "ebins.csv", # csvs containing ELFIN measurements
#                                                 DateTime(2020,9,22,9,16,25), DateTime(2020,9,22,9,16,34)) # time to sample from ELFIN measurements
# @time elfin_measurements_092220 = extract_idl_csv("092220_time.csv", "092220_prec.csv", "ebins.csv", # csvs containing ELFIN measurements
#                                                 DateTime(2020,9,22,9,15,48), DateTime(2020,9,22,9,16,3)) # time to sample from ELFIN measurements
# @time elfin_measurements_all_092220 = extract_idl_csv("092220_time.csv", "092220_prec.csv", "ebins.csv", # csvs containing ELFIN measurements
#                                                 DateTime(2020,9,22,9,15,48), DateTime(2020,9,22,9,16,45)) # time to sample from ELFIN measurements

x1 = mms_lowerlat_092220.precipitating_fluxes_mean.-mms_lowerlat_092220.precipitating_fluxes_minus .+ 1.
x2 = mms_lowerlat_092220.precipitating_fluxes_mean.+mms_lowerlat_092220.precipitating_fluxes_plus .+ 1.
y1 = mms_lolat_092220.precipitating_fluxes_mean.-mms_lolat_092220.precipitating_fluxes_minus .+ 1.
y2 = mms_lolat_092220.precipitating_fluxes_mean.+mms_lolat_092220.precipitating_fluxes_plus .+ 1.
z1 = mms_hilat_092220.precipitating_fluxes_mean.-mms_hilat_092220.precipitating_fluxes_minus .+ 1.
z2 = mms_hilat_092220.precipitating_fluxes_mean.+mms_hilat_092220.precipitating_fluxes_plus .+ 1.
e1 = elfin_measurements_092220[2].-elfin_error_092220
e2 = elfin_measurements_092220[2].+elfin_error_092220

plot(Egrid, equatorial_fluxes_092220,label = L"\mathrm{equatorial\ flux\ (MMS\ 1\ measured)}", color = bipride_pink, linewidth=2, markershape=:circle);
plot!(Egrid,mms_lowerlat_092220.precipitating_fluxes_mean, fillrange=x1, fillalpha = 0.2, color = c4, label=false)
plot!(Egrid,mms_lowerlat_092220.precipitating_fluxes_mean, fillrange=x2, fillalpha = 0.2, color = c4, label=false)
plot!(Egrid, mms_lowerlat_092220.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\lambda_2 = 15\degree\mathrm{,\ occurence\ rate} = 6\%)",
        # yerror=(mms_lowerlat_092220.precipitating_fluxes_minus, mms_lowerlat_092220.precipitating_fluxes_plus),
        color = c4, alpha = .7, markerstrokecolor = c4, markerstrokewidth = 2, markershape=:o, markersize = 2);
plot!(Egrid,mms_lolat_092220.precipitating_fluxes_mean, fillrange=y1, fillalpha = 0.2, color = c3, label=false)
plot!(Egrid,mms_lolat_092220.precipitating_fluxes_mean, fillrange=y2, fillalpha = 0.2, color = c3, label=false)
plot!(Egrid, mms_lolat_092220.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\lambda_2 = 20\degree\mathrm{,\ occurence\ rate} = 6\%)",
        # yerror=(mms_lolat_092220.precipitating_fluxes_minus, mms_lolat_092220.precipitating_fluxes_plus),
        color = c3, alpha = .7, markerstrokecolor = c3, markerstrokewidth = 2, markershape=:o, markersize = 2);
plot!(Egrid,mms_hilat_092220.precipitating_fluxes_mean, fillrange=z1, fillalpha = 0.2, color = bipride_orange, label=false)
plot!(Egrid,mms_hilat_092220.precipitating_fluxes_mean, fillrange=z2, fillalpha = 0.2, color = bipride_orange, label=false)
plot!(Egrid, mms_hilat_092220.precipitating_fluxes_mean,
        label=L"\mathrm{precipitating\ flux\ }(\mathrm{simulated\ with\ }\delta\lambda_2 = 30\degree\mathrm{,\ occurence\ rate} = 2\%)",
        # yerror=(mms_hilat_092220.precipitating_fluxes_minus, mms_hilat_092220.precipitating_fluxes_plus),
        color = bipride_orange, alpha = .7, markerstrokecolor = bipride_orange, markerstrokewidth = 2, markershape=:o, markersize = 2);
plot!(elfin_measurements_092220[1][1:6], elfin_measurements_092220[2][1:6], fillrange = e1[1:end-5], fillalpha = 0.2, color = c5, label=false)
plot!(elfin_measurements_092220[1][1:6], elfin_measurements_092220[2][1:6], fillrange = e2[1:end-5], fillalpha = 0.2, color = c5, label=false)
plot!(elfin_measurements_092220[1][1:6], elfin_measurements_092220[2][1:6],
        yerror=elfin_error_092220[1:6].+10,
        label = L"\mathrm{precipitating\ flux\ (ELFIN\ A\ measured)}",
        color = c5, marker = stroke(3,c5), linewidth=4, markersize = 3);
plot!(ylim =(1e1,1e10), xlim=(50,800), yscale=:log10, legend=:topright)
plot!(xlabel=L"\mathrm{Energy\ (keV)}", ylabel=L"\mathrm{Flux\ (1/cm^{2}/s/sr/MeV)}", title=L"\mathrm{Precipitating\ Flux\ Comparison\ of\ \delta\lambda_2\ on\ 9/22/20}",
xtickfontsize=12, ytickfontsize=12, xguidefontsize=16, yguidefontsize=16, legendfontsize=10, titlefontsize=16);
plot3 = plot!(dpi = 500,size=(800,450), margin=3mm, bottom_margin=4mm)
savefig(plot3, "images/092220_flux_comparison.pdf")
savefig(plot3, "images/092220_flux_comparison.png")

max_stars_repo_name: ethantsai/nlwhistlers

max_issues_repo_name: ethantsai/nlwhistlers

max_forks_repo_name: ethantsai/nlwhistlers