text:
## Main function to generate pseudo-Data

function GenPseudoDat(model_def, pseudo_def)

    # First make a folder where to save the results
    cudi = pwd(); #@__DIR__;
    if !isdir(string(cudi,"\\Results"))
        mkdir(string(cudi,"\\Results"))
    end
    if !isdir(string(cudi, "\\Results\\", model_def["NameF"],"_",today()))
        mkdir(string(cudi, "\\Results\\", model_def["NameF"],"_",today()))
    end

    # Check which of the 2 types of simul_def have been given and make a choice on how to proceed
    if haskey(pseudo_def, "switchT")
        pseudo_def = checkStructPseudoDat(model_def, pseudo_def)
    else
        pseudo_def = extractPseudoDatCSV(model_def, pseudo_def)
    end

    simul_def = defSimulStruct();

    for (i,j) in simul_def
        simul_def[i] = pseudo_def[i];
    end
    # Make empty structure where the simulations will be
    simuls = Dict();

    # Get ODEs function
    OEDsFun = Symbol(string(join(model_def["NameF"]),"_SolveAll"));

    # loop over the different experiments to simulate
    for i in 1:simul_def["Nexp"]
        # Define simulation parameters
        global ts = collect(0.0:round(simul_def["finalTime"][i]))'; # 24 hour experiment
        global sp = [convert(Int,v) for v in (round.(simul_def["switchT"][i])')]; # 8 step experiment for now
        global ivss = simul_def["y0"][i]; # Experimental values for Y0 (IPTG, aTc, RFP, GFP)

        if model_def["Y0eqs"] != [] # ON inputs
            global pre = simul_def["preInd"][i];
        else
            global pre = [];
        end
        global samps = convert.(Int, round.(simul_def["tsamps"][i]));
        global theta = simul_def["theta"];

        # Inputs case
        if model_def["nInp"] != 0
            r = convert.(Int, 1:model_def["nInp"]:(length(simul_def["uInd"][i])));
            global inputs = zeros(convert.(Int,length(simul_def["uInd"][i])));
            for j in 1:convert.(Int,length(simul_def["uInd"][i])/model_def["nInp"])
                for k in 0:(model_def["nInp"]-1)
                    inputs[r[j]+k] = simul_def["uInd"][i][j,(k+1)];
                end
            end
        else
            r=0
        inputs=[];
        end

        simul = @eval $OEDsFun(ts, theta, sp, inputs, ivss, samps, pre);

        simuls[string("Exp_", i)] = simul;
    end

    # Now extract the entries to which the pseudodata will be generated
    pseudoDme = Dict();
    pseudoDsd = Dict();
    simObs = Dict();
    for i in 1:simul_def["Nexp"]
        pseudoDme[string("PDExp_", i)] = zeros(size(simuls[string("Exp_", i)])[1], length(pseudo_def["Obs"]), size(simuls[string("Exp_", 1)])[3]);
        pseudoDsd[string("PDExp_", i)] = zeros(size(simuls[string("Exp_", i)])[1], length(pseudo_def["Obs"]), size(simuls[string("Exp_", 1)])[3]);
        simObs[string("PDExp_", i)] = zeros(size(simuls[string("Exp_", i)])[1], length(pseudo_def["Obs"]), size(simuls[string("Exp_", 1)])[3]);
        if typeof(pseudo_def["Obs"])==Array{Int,1}
            simObs[string("PDExp_", i)] = simuls[string("Exp_", i)][:,(pseudo_def["Obs"]),:];
        else
            coun = 1;
            for j in 1:length(pseudo_def["Obs"])
               if !(occursin.("+", pseudo_def["Obs"][j])) && !(occursin.("-", pseudo_def["Obs"][j])) &&
               !(occursin.("/", pseudo_def["Obs"][j])) && !(occursin.("*", pseudo_def["Obs"][j])) &&
                !(occursin.("^", pseudo_def["Obs"][j]))
                    tind = findall(x->x==pseudo_def["Obs"][j], model_def["stName"])[1];
                    simObs[string("PDExp_", i)][:,coun,:] = simuls[string("Exp_", i)][:,tind,:];
                    coun +=1;
                else #This is the complicated stuff where we need to evaluate the function written as a string
                    if (occursin(".+", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], ".+"=>" .+ ")
                    elseif (occursin("+", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], "+"=>" .+ ")
                    end

                    if (occursin(".-", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], ".-"=>" .- ")
                    elseif (occursin("-", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], "-"=>" .- ")
                    end

                    if (occursin("./", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], "./"=>" ./ ")
                    elseif (occursin("/", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], "/"=>" ./ ")
                    end

                    if (occursin(".*", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], ".*"=>" .* ")
                    elseif (occursin("*", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], "*"=>" .* ")
                    end

                    if (occursin(".^", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], ".^"=>" .^ ")
                    elseif (occursin("^", pseudo_def["Obs"][j]))
                        pseudo_def["Obs"][j] = replace(pseudo_def["Obs"][j], "^"=>" .^ ")
                    end

                    tmp1 = pseudo_def["Obs"][j];
                    for k in 1:model_def["nStat"]
                        if occursin(model_def["stName"][k], pseudo_def["Obs"][j])
                            tmp1 = replace(tmp1, model_def["stName"][k]=>string(simuls[string("Exp_",i)][:,k,:]))
                        end
                    end
                    simObs[string("PDExp_", i)][:,coun,:] = eval(Meta.parse(tmp1));
                    coun +=1;
                end
            end
        end
    end

    # Generate pseudo-data with noise

    for j in 1:simul_def["Nexp"]
        for stat in 1:length(pseudo_def["Obs"])
            if length(size(pseudo_def["theta"]))==2
                for p in 1:convert(Int, length(pseudo_def["theta"])/model_def["nPar"])
                    if pseudo_def["NoiseType"] == "hetero"
                        pseudoDsd[string("PDExp_", j)][:,stat,p] = simObs[string("PDExp_", j)][:,stat,p] .* pseudo_def["Noise"][stat];
                    elseif pseudo_def["NoiseType"] == "homo"
                        pseudoDsd[string("PDExp_", j)][:,stat,p] = ones(length(simObs[string("PDExp_", j)][:,stat,p])) .* pseudo_def["Noise"][stat];
                    end
                    pseudoDme[string("PDExp_", j)][:,stat,p] = simObs[string("PDExp_", j)][:,stat,p] .+ rand(MvNormal(zeros(size(simObs[string("PDExp_", j)][:,stat,p])[1]),pseudoDsd[string("PDExp_", j)][:,stat,p][:]))
                    if pseudoDsd[string("PDExp_", j)][1,stat,p] == 0 # This checks for the case wher ethe user gives y0 and this is 0. This will output a mean and standard deviation of 0 which will generate issues in lekelihood computation (MLE and Stan inference)
                        pseudoDsd[string("PDExp_", j)][1,stat,p] = 1e-20;
                    end
                    if pseudoDme[string("PDExp_", j)][1,stat,p] == 0
                        pseudoDme[string("PDExp_", j)][1,stat,p] = 1e-10;
                    end
                    for r in 1:size(pseudoDme[string("PDExp_", j)])[1]
                        if pseudoDme[string("PDExp_", j)][r,stat,p] == 0 # If you have homoscedastic noise and high, this can generate negative values which do not make much sense. A proxi to a truncation of the distribution I suppsose.
                            pseudoDme[string("PDExp_", j)][r,stat,p] = 1e-10;
                        end
                    end
                end
            else
                    p=1;
                    if pseudo_def["NoiseType"] == "hetero"
                        pseudoDsd[string("PDExp_", j)][:,stat,p] = simObs[string("PDExp_", j)][:,stat,p] .* pseudo_def["Noise"][stat];
                    elseif pseudo_def["NoiseType"] == "homo"
                        pseudoDsd[string("PDExp_", j)][:,stat,p] = ones(length(simObs[string("PDExp_", j)][:,stat,p])) .* pseudo_def["Noise"][stat];
                    end
                    pseudoDme[string("PDExp_", j)][:,stat,p] = simObs[string("PDExp_", j)][:,stat,p] .+ rand(MvNormal(zeros(size(simObs[string("PDExp_", j)][:,stat,p])[1]),pseudoDsd[string("PDExp_", j)][:,stat,p][:]))
                    if pseudoDsd[string("PDExp_", j)][1,stat,p] == 0
                        pseudoDsd[string("PDExp_", j)][1,stat,p] = 1e-20;
                    end
                    if pseudoDme[string("PDExp_", j)][1,stat,p] == 0
                        pseudoDme[string("PDExp_", j)][1,stat,p] = 1e-10;
                    end

                    for r in 1:size(pseudoDme[string("PDExp_", j)])[1]
                        if pseudoDme[string("PDExp_", j)][r,stat,p] == 0 # If you have homoscedastic noise and high, this can generate negative values which do not make much sense. A proxi to a truncation of the distribution I suppsose.
                            pseudoDme[string("PDExp_", j)][r,stat,p] = 1e-10;
                        end
                    end
            end
        end
    end
    # Save path
    pseudo_def["savepath"] = string(cudi, "\\Results\\", model_def["NameF"],"_",today());
    pseudo_def["savename"] = string(model_def["NameF"],"_",today(), "_PseudoDataResults_",pseudo_def["flag"],".jld");

    # Put all results into one same structure and save
    pseudo_res = Dict()
    pseudo_res["Sims"] = simuls;
    pseudo_res["SimsObs"] = simObs;
    pseudo_res["PData"] = pseudoDme;
    pseudo_res["PError"] = pseudoDsd;

    save(string(pseudo_def["savepath"], "\\", pseudo_def["savename"]), "PseudoData", pseudo_res, "model_def", model_def, "pseudo_def", pseudo_def)

    println("")
    println("----------------------------------------- RESULTS -----------------------------------------")
    println("Pseudo-Data results are saved in the directory: ")
    println(string("                 ", pseudo_def["savepath"]))
    println(string("Under the name ",pseudo_def["savename"]))
    println("--------------------------------------------------------------------------------------")
    println("")

    if simul_def["plot"] == true
        plotPseudoDatODE(pseudo_res,model_def,pseudo_def)
        println("")
        println("----------------------------------------- PLOTS -----------------------------------------")
        println("Pseudo-Data PLOTS are saved in the directory: ")
        println(string("                 ", pseudo_def["savepath"]))
        println(string("Under the name PlotPseudoDat_Exp(i)_",pseudo_def["flag"],".png"))
        println(string("If more than one sample from the parameters are used, the plot will be the average between all traces."))
        println("--------------------------------------------------------------------------------------")
        println("")
    end

    PDatCSVGen(pseudo_res,model_def,pseudo_def);
    println("")
    println("----------------------------------------- CSVs -----------------------------------------")
    println("Pseudo-Data CSVs with results have been generated in the directory: ")
    println(string("                 ", pseudo_def["savepath"], "\\PseudoDataFiles"))
    println(string("Under the names: "))
    println(string("        Simulations: ", model_def["NameF"],"_EXP(i)", "_",
                    pseudo_def["flag"],"_Simulations.csv"))
    println(string("        Observables: ", model_def["NameF"],"_EXP(i)", "_",
                    pseudo_def["flag"],"_Observables.csv"))
    println(string("        Event Inputs (only if the model has inputs): ", model_def["NameF"],"_EXP(i)", "_",
                    pseudo_def["flag"],"_Event_Inputs.csv"))
    println("--------------------------------------------------------------------------------------")
    println("")

    return pseudo_res, model_def, pseudo_def

end

## Pseudodata general plot function

function plotPseudoDatODE(pseudo_res,model_def,pseudo_def)
    cudi = pwd(); # @__DIR__;
    for i in 1:pseudo_def["Nexp"]
        # Title for each subplot
        tit = "";
        yl1 = "";
        for k in 1:length(pseudo_def["Obs"])
            if typeof(pseudo_def["Obs"]) == Array{Int,1};
                tit = hcat(tit, string(model_def["stName"][pseudo_def["Obs"][k]]));
            else
                tit = hcat(tit, string(pseudo_def["Obs"][k]));
            end

            yl1 = hcat(yl1, "y");
        end
        tit = tit[:,2:end];
        yl1 = yl1[:,2:end];

        titu = "";
        yl2 = "";
        if model_def["nInp"] != 0
            for k in 1:model_def["nInp"]
                titu = hcat(titu, string(model_def["inpName"][k]));
                yl2 = hcat(yl2, "u");
            end
            titu = titu[:,2:end];
            yl2 = yl2[:,2:end];
        end

        tuu = hcat(tit, titu);
        yuu = hcat(yl1, yl2);

        # Elements for the inducers
        tu = Array{Array{Any,1}}(undef, length(pseudo_def["Obs"])+model_def["nInp"])
        [tu[k] = [] for k in 1:length(pseudo_def["Obs"])]
        su = Array{Array{Any,1}}(undef, length(pseudo_def["Obs"])+model_def["nInp"])
        [su[k] = [] for k in 1:length(pseudo_def["Obs"])]
        for k in 1:model_def["nInp"]
            tu[k+length(pseudo_def["Obs"])] = round.(pseudo_def["switchT"][i]);
            su[k+length(pseudo_def["Obs"])] = vcat(pseudo_def["uInd"][i][:,(k)], pseudo_def["uInd"][i][end,(k)])
        end

        t = round.(pseudo_def["tsamps"][i]);
        m = pseudo_res["PData"][string("PDExp_",i)];
        s = pseudo_res["PError"][string("PDExp_",i)];

        if length(size(pseudo_def["theta"])) == 1
            pl = plot(t, m[:,:,1],yerror = s[:,:,1], layout=length(pseudo_def["Obs"])+model_def["nInp"], label = "", title = tit, xlabel = "time", ylabel = yuu,
            color = "blue", size = [2000, 1200])

            pl = plot!(tu, su, layout=model_def["nStat"]+model_def["nInp"], label = "", linetype=:step, title = tuu);
        else
            m2 = mean(m, dims = 3);
            s2 = std(m, dims = 3);
            pl = plot(t, m2[:,:,1],yerror = s2[:,:,1], layout=length(pseudo_def["Obs"])+model_def["nInp"], label = "Average", title = tit, xlabel = "time", ylabel = yuu,
            color = "blue", size = [2000, 1200])

            pl = plot!(tu, su, layout=model_def["nStat"]+model_def["nInp"], label = "", linetype=:step, title = tuu);
        end

        savefig(string(cudi, "\\Results\\", model_def["NameF"],"_",today(), "\\PlotPseudoDat_Exp", i,"_", pseudo_def["flag"], ".png"));

    end

end

## Simulation Data structure Dict
function defPseudoDatStruct()

    # The structure of things user should give
    pseudo_def = Dict()

    pseudo_def["Nexp"] = []; # Integer indicating the number of experiments to be simulated
    pseudo_def["finalTime"] = []; # -> Vector of final times for each simulation (initial time will allways be asumed as 0, so please consider that)
    pseudo_def["switchT"] = []; # Array with the switching times of the inducer in the simulation (time 0 and final time need to be considered)
    pseudo_def["y0"] = []; # Array of Y0s for the simulations for each experiment. If you are computing the steady state this vector might not be used, however you still need to introduce it
    pseudo_def["preInd"] = []; # Level of the inducer in the ON. This entry can be empty if no Steady State equations are given and only the Y0 given by the user is considered.
    pseudo_def["uInd"] = []; # Array with the values for the inducer for each experiment at each step
    pseudo_def["theta"] = []; # Vector/Matrix with the parameter samples or directory and file location of CSV file with them. If a matrix is given, pseudo-data will be generated for each one of the theta samples
    pseudo_def["tsamps"] = []; # Array of Sampling times vectors
    pseudo_def["plot"] = []; # Bollean or yes/no string to save the resulting simulations in the results directory (false will be considered as default)
    pseudo_def["flag"] = []; # String to attach a unique flag to the generated files so it is not overwritten. If empty, nothing will be added.

    pseudo_def["Obs"] = []; # States that are observable. This is either a vecotr of strings (that has the same order asmodel_def["stName"] and model_def["eqns"]), a vector of integers insidating which states are observable or the string "All" if all states are observables. This can also be an experssion combining states (Only +,-,*,\ and ^ will be considered)
    pseudo_def["NoiseType"] = []; # String indicating if the desired noise is homoscedastic or heteroscedastic. The allowed strings are homo, homoscedastic or hetero, heteroscedastic (case independent). If empty, default will be heteroscedastic. Note that for high homoscedastic noise data could become negative (which to me does not make much sense), for now this will be set to a value close to 0.
    pseudo_def["Noise"] = []; # Percentage of heteroscedastic noise (introduced as value from 0 to 1). If empty 10% heteroscedastic noise will be assumed. In the homoscedastic case, if empty a standard deviation of 1 will be assumes. This has to be a vector of noise values for each observable.

    return(pseudo_def)

end

## Check pseudo-data structure function
# Function to check the content of the fields the user has given
function checkStructPseudoDat(model_def, pseudo_def)

    # Check taht all the dictionary entries are correct
    entries = ["Nexp", "finalTime", "switchT", "y0", "preInd", "uInd", "theta", "tsamps", "plot", "flag", "Obs", "NoiseType", "Noise"]
    if symdiff(entries,keys(pseudo_def))!=[] && symdiff(entries,keys(pseudo_def)) != ["savepath", "savename"]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the entries of the dictionary, there is something wrong...")
        println(symdiff(entries,keys(pseudo_def)))
        return
    end

    # Check one to be sure there is no empty entry
    nee = ["Nexp", "finalTime", "switchT", "y0", "uInd", "theta", "tsamps", "Obs"]; # No empty entries
    for i in 1:length(nee)
        if pseudo_def[nee[i]] == []
            println("-------------------------- Process STOPPED!!! --------------------------")
            println(string("Plese, check the field ", nee[i], ", this cannot be empty"))
            return
        end
    end

    # Checks to see if the content of the fields is the expected (also consier if single value entries have been given as a number or a vector)
    if ((typeof(pseudo_def["Nexp"][1])!=Int) || length(pseudo_def["Nexp"])!=1) && (typeof(pseudo_def["Nexp"])!=Int)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field Nexp! This should be an integer. ")
        return
    elseif (typeof(pseudo_def["finalTime"]) != Array{Int,1}) && (typeof(pseudo_def["finalTime"]) != Array{Float64,1}) && (typeof(pseudo_def["finalTime"]) != Array{Int32,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field finalTime! This should be a vector of final times. ")
        return
    elseif (typeof(pseudo_def["switchT"]) != Array{Array{Int,1},1}) && (typeof(pseudo_def["switchT"]) != Array{Array{Float64,1},1}) && (typeof(pseudo_def["switchT"]) != Array{Array{Float32,1},1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field switchT! This should be an array of arrays of numbers. ")
        return
    elseif (typeof(pseudo_def["y0"]) != Array{Array{Int,1},1}) && (typeof(pseudo_def["y0"]) != Array{Array{Float64,1},1}) && (typeof(pseudo_def["y0"]) != Array{Array{Float32,1},1}) &&
        (typeof(pseudo_def["y0"]) != Array{Array{Int,N} where N,1}) && (typeof(pseudo_def["y0"]) != Array{Array{Float64,N} where N,1}) && (typeof(pseudo_def["y0"]) != Array{Array{Float32,N} where N,1}) &&
        (typeof(pseudo_def["y0"]) != Array{Array{Int,2},1}) && (typeof(pseudo_def["y0"]) != Array{Array{Float64,2},1}) && (typeof(pseudo_def["y0"]) != Array{Array{Float32,2},1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field y0! This should be an array of arrays of numbers. ")
        return
    elseif (typeof(pseudo_def["preInd"]) != Array{Array{Int,1},1}) && (typeof(pseudo_def["preInd"]) != Array{Array{Float64,1},1}) && (typeof(pseudo_def["preInd"]) != Array{Array{Float32,1},1}) && pseudo_def["preInd"] != []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field preInd! This should be an array of arrays of numbers. ")
        return
    # elseif (typeof(pseudo_def["uInd"]) != Array{Array{Int,1},1}) && (typeof(pseudo_def["uInd"]) != Array{Array{Float64,1},1}) && (typeof(pseudo_def["uInd"]) != Array{Array{Float32,1},1})
    #     println("-------------------------- Process STOPPED!!! --------------------------")
    #     println("Please, check the field uInd! This should be an array of arrays of numbers. ")
    #     return
    elseif (typeof(pseudo_def["theta"]) != Array{Float64,1}) && (typeof(pseudo_def["theta"]) != Array{Float32,1}) &&
        (typeof(pseudo_def["theta"]) != Array{Float64,2}) && (typeof(pseudo_def["theta"]) != Array{Float32,2}) &&
        (typeof(pseudo_def["theta"]) != Array{String,1}) && (typeof(pseudo_def["theta"]) != String)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field theta! This should be a vector, matrix or string. ")
        return
    elseif (typeof(pseudo_def["tsamps"]) != Array{Array{Int,1},1}) && (typeof(pseudo_def["tsamps"]) != Array{Array{Float64,1},1}) && (typeof(pseudo_def["tsamps"]) != Array{Array{Float32,1},1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field tsamps! This should be an array of arrays of numbers. ")
        return
    elseif (pseudo_def["plot"] != [])
        try
            if (typeof(pseudo_def["plot"]) != Bool) && (typeof(pseudo_def["plot"]) != Array{Bool,1}) &&
            (typeof(pseudo_def["plot"]) != String) && (typeof(pseudo_def["plot"]) != Array{String,1})
                println("-------------------------- Process STOPPED!!! --------------------------")
                println("Please, check the field plot! This should be true, fals, \"Yes\", \"yes\", \"No\", \"no\" or an empty array ")
                return
            end
        catch
             if (typeof(pseudo_def["plot"]) != Array{Bool,1}) &&
            (typeof(pseudo_def["plot"]) != Array{String,1})
                println("-------------------------- Process STOPPED!!! --------------------------")
                println("Please, check the field plot! This should be true, fals, \"Yes\", \"yes\", \"No\", \"no\" or an empty array ")
                return
            end
        end
    elseif (typeof(pseudo_def["flag"]) != Array{String,1}) && (typeof(pseudo_def["flag"]) != String) && ((pseudo_def["flag"]) != [])
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field flag! This should be a vector, matrix or string. ")
        return
    elseif (typeof(pseudo_def["Obs"]) != Array{String,1}) && (typeof(pseudo_def["Obs"]) != Array{Int,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field Obs! This should be a vector of strings or integers. ")
        return
    elseif (typeof(pseudo_def["NoiseType"]) != Array{String,1}) && (typeof(pseudo_def["NoiseType"]) != String) && ((pseudo_def["NoiseType"]) != [])

    elseif (typeof(pseudo_def["Noise"]) != Array{Float64,1}) && (typeof(pseudo_def["Noise"]) != Array{Float32,1})  &&
        ((pseudo_def["NoiseType"]) != []) && (typeof(pseudo_def["Noise"]) != Array{Int,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field Noise! This should be a vector of error percentages from 0 to 1. ")
        return
    end

    # Extract necessary elements to ease generalisation

    if typeof(pseudo_def["Nexp"]) == Array{Int,1}
        pseudo_def["Nexp"] = pseudo_def["Nexp"][1];
    end

    if pseudo_def["plot"]==[true] || pseudo_def["plot"]==["Yes"] || pseudo_def["plot"]==["yes"] || pseudo_def["plot"]=="Yes" || pseudo_def["plot"]=="yes"
        pseudo_def["plot"]=true
    elseif pseudo_def["plot"]==[false] || pseudo_def["plot"]==["No"] || pseudo_def["plot"]==["no"] || pseudo_def["plot"]=="No" || pseudo_def["plot"]=="no" || pseudo_def["plot"]==[]
        pseudo_def["plot"]=false
    end

    if (typeof(pseudo_def["theta"]) == Array{String,1})
        pseudo_def["theta"] = pseudo_def["theta"][1];
    end

    if typeof(pseudo_def["flag"]) == Array{String,1}
        pseudo_def["flag"] = pseudo_def["flag"][1];
    elseif typeof(pseudo_def["flag"]) == String
        pseudo_def["flag"] = pseudo_def["flag"];
    elseif ((pseudo_def["flag"]) == [])
        pseudo_def["flag"] = "";
    end

    if (typeof(pseudo_def["NoiseType"]) == Array{String,1})
        pseudo_def["NoiseType"] = pseudo_def["NoiseType"][1];
    end

    # Check that all the contents make sense
    if (pseudo_def["Nexp"] != length(pseudo_def["finalTime"])) || (pseudo_def["Nexp"] != length(pseudo_def["switchT"])) ||
        (pseudo_def["Nexp"] != length(pseudo_def["y0"])) || (pseudo_def["Nexp"] != length(pseudo_def["uInd"])) ||
        (pseudo_def["Nexp"] != length(pseudo_def["tsamps"]))
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, some field does not match the number of experiments selected. Check that the number of")
        println("entries for finalTime, switchT, y0, uInd or tsamps matches the number of experiments in Nexp.")
        return
    end

    if (typeof(pseudo_def["theta"]) == Array{Float64,1}) || (typeof(pseudo_def["theta"]) == Array{Float32,1})
        if length(pseudo_def["theta"]) != model_def["nPar"]
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Number of parameters introduced does not match the specified")
            return
        end
    elseif (typeof(pseudo_def["theta"]) == Array{Float64,2}) || (typeof(pseudo_def["theta"]) == Array{Float32,2})
        if size(pseudo_def["theta"])[1] != model_def["nPar"] && size(pseudo_def["theta"])[2] != model_def["nPar"]
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Number of parameters introduced does not match the specified")
            return
        end
    elseif (typeof(pseudo_def["theta"]) == String)
        if pseudo_def["theta"][end-3:end] != ".csv"
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, add the .csv termination to the theta file! And no spaces after!")
            return
        end
        if isfile(pseudo_def["theta"])
            pseudo_def["theta"] = Matrix(CSV.read(pseudo_def["theta"]));
        else
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Sorry, but the file path you introduced does not exist!")
            return
        end
    end

    for i in 1:pseudo_def["Nexp"]
        if pseudo_def["tsamps"][i][end] > pseudo_def["finalTime"][i]
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, check finalTime. You have selected a sampling point past it.")
            return
        end
        if model_def["nInp"] != 0
            if length(pseudo_def["uInd"][i]) != (length(pseudo_def["switchT"][i])-1)
                println("-------------------------- Process STOPPED!!! --------------------------")
                println("Please, check uInd and switchT. Number of steps does not match the number of values for the inputs.")
                return
            end
        end

        if length(pseudo_def["y0"][i])/model_def["nStat"] == 1
            if length(pseudo_def["y0"][i]) != model_def["nStat"]
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but Y0 does not match the number or states in experiment ", i))
                return
            end
        else
            if size(pseudo_def["y0"][i])[1] != model_def["nStat"] && size(pseudo_def["y0"][i])[2] != model_def["nStat"]
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but Y0 does not match the number or states in experiment ", i))
                return
            end
            if length(pseudo_def["y0"][i])/model_def["nStat"] != length(pseudo_def["theta"])/model_def["nPar"]
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but Y0 does not match the dimensions of theta "))
                return
            end
            if size(pseudo_def["y0"][i])[1] == size(pseudo_def["y0"][i])[2]
                println("-------------------------- WARNING --------------------------")
                println(string("Sorry, but the number of rows and columns of the y0 matrix in experiment ",i," is the same, so the checks on "))
                println("correct orientation will not work. Please make sure that the dimensions follow: ")
                println("y0[samples, states]")
            end
        end
    end

    if (model_def["Y0eqs"] != []) && pseudo_def["preInd"]==[]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("You specified computation of Y0 as steady-state but you have not introduced an inducer value for it!")
        return
    end

    if pseudo_def["NoiseType"] == []
        pseudo_def["NoiseType"] = "hetero"
    end

    pseudo_def["NoiseType"] = lowercase(pseudo_def["NoiseType"]); # To avoid issues with uper cases and that

    if (pseudo_def["NoiseType"] == "homoscedastic")
        pseudo_def["NoiseType"] = "homo";
    elseif (pseudo_def["NoiseType"] == "heteroscedastic")
        pseudo_def["NoiseType"] = "hetero";
    end

    if pseudo_def["Noise"] == []
        pseudo_def["Noise"] = ones(length(pseudo_def["Obs"])) .* 0.1;
    end

    if (pseudo_def["NoiseType"] == "hetero") && (maximum(pseudo_def["Noise"]) > 1)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, but the noise percentage has to be between 0 and 1!")
        return
    end

    if (minimum(pseudo_def["Noise"]) < 0)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, but the noise percentage or standard deviation is negative and that CANNOT be!")
        return
    end

    if length(pseudo_def["Obs"]) != length(pseudo_def["Noise"])
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, but you need to assign an error value for each observable selected!")
        return
    end

    # Check warning in case the matrix introduced is symmetric
    if length(size(pseudo_def["theta"])) == 2
        if size(pseudo_def["theta"])[1] == size(pseudo_def["theta"])[2]
            println("-------------------------- WARNING --------------------------")
            println("Sorry, but the number of rows and columns of the theta matrix is the same, so the checks on ")
            println("correct orientation will not work. Please make sure that the dimensions follow: ")
            println("theta[samples, parameters]")
        end
    end

    # Check that no step is no smaller than 2 unit of time
    for i in 1:pseudo_def["Nexp"]
        for j in 1:length(pseudo_def["uInd"][i])
            if (pseudo_def["switchT"][i][j+1]-pseudo_def["switchT"][i][j])<=4
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but 2 of the steps in experiment ", i, " are too close. This package cannot "))
                println("handle this for now. We are working on it!")
                return
            end
        end
    end

    # Check that the inputs for the experiment are given as collumns
    try
        if model_def["nInp"]>1
            for i in 1:pseudo_def["Nexp"]
                if size(pseudo_def["uInd"][i])[2] != model_def["nInp"];
                    pseudo_def["uInd"][i] = pseudo_def["uInd"][i]';
                end
            end
        end
    catch
        println("-------------------------- Process STOPPED!!! --------------------------")
        println(string("Sorry, but there is some issue with the contents of the field uInd."))
        return
    end

    if pseudo_def["Obs"] == Array{Int,1}
        if (maximum(pseudo_def["Obs"]) > model_def["nStat"]) || (length(pseudo_def["Obs"]) > model_def["nStat"])
            println("-------------------------- Process STOPPED!!! --------------------------")
            println(string("Sorry, but there is some issue with the contents of the field Obs. It seems that "))
            println("     you have selected more observables than states or an index higher than the number of observables")
            return
        end
    end

    if typeof(pseudo_def["Obs"]) == Array{String,1}
        if convert(Bool,sum(occursin.("+", pseudo_def["Obs"]))) || convert(Bool,sum(occursin.("-", pseudo_def["Obs"]))) ||
           convert(Bool,sum(occursin.("/", pseudo_def["Obs"]))) ||
           convert(Bool,sum(occursin.("*", pseudo_def["Obs"]))) || convert(Bool,sum(occursin.("^", pseudo_def["Obs"])))
            nothing
        else
            if sum(occursin.(model_def["stName"], pseudo_def["Obs"])) == 0
                println("-------------------------- Process STOPPED!!! --------------------------")
                println(string("Sorry, but there is some issue with the contents of the field Obs."))
                println("It seems that the observable(s) selected do not match any state")
                return
            else
                pseudo_def["Obs"] = sort([findall(x->x==pseudo_def["Obs"][i], model_def["stName"])[1] for i in 1:length(pseudo_def["Obs"])])
            end
        end

    elseif typeof(pseudo_def["Obs"]) == Array{Int,1}
        pseudo_def["Obs"] = sort(pseudo_def["Obs"]);
    end

    return(pseudo_def)
end

## CSV introduction description
# function fileStructInfo()
#
#     println("ALL THE FILES MUST BE CSV FILES!!!")
#     println("")
#     println("---------------------------------------------------------------------------------------------------------")
#     println("")
#     println("The observables file entry should have the following structure: ")
#     println("Column 1: Sampling Times for the simulations")
#     println("Column 2 + Number of states: Y0 value for each each state (in order) in each column. ")
#     println("           Value might need to be repeated across all the column, but only the first row will be considered")
#     println("")
#     println("---------------------------------------------------------------------------------------------------------")
#     println("")
#     println("The event inputs file entry should have the following structure: ")
#     println("Column 1: Each Switching time (final time not considered)")
#     println("Column 2: Final time (can be repeated for as many entrances as column 1 has)")
#     println("Column 3 + number of inducers: Repeated column for the value of the inducers in the ON. ")
#     println("              If no pre-inducer will be used in the simulations, any number could be used")
#     println("Column 4 + number of inducers: Value of each inducer at each switching time.")
#     println("")
#     println("This is an Example of an experiment with 3 steps and 2 inducers: ")
#     println(" |Switching|  |FinalTime|  |IPTGpre|  |aTcPre|  |  IPTG  |  |  aTc  |")
#     println("      0           1350         1          0          1          0    ")
#     println("     400          1350         1          0         0.3         50   ")
#     println("     700          1350         1          0         0.7         35   ")
#     println("     950          1350         1          0         0.1        100   ")
#     println("")
#     println("---------------------------------------------------------------------------------------------------------")
#     println("")
#     println("For the main directory, remember that in Julia you have to type \\ twice!")
#
# end

## CSV introduction structure
function defPseudoDatStructFiles()

    fileStructInfo()

    # The structure of things user should give
    pseudo_def = Dict()

    pseudo_def["ObservablesFile"] = []; # File containing information about the states Y0 and sampling times
    pseudo_def["EventInputsFile"] = []; # -> File containing all the information about the stimuli
    pseudo_def["theta"] = []; # Vector/Matrix with the parameter samples or directory and file location of CSV file with them
    pseudo_def["MainDir"] = []; # Directory path where all files are located (to avoid repetition). This can be empty.
    pseudo_def["plot"] = [];
    pseudo_def["flag"] = [];

    pseudo_def["Obs"] = []; # States that are observable. This is either a vecotr of strings (that has the same order asmodel_def["stName"] and model_def["eqns"]), a vector of integers insidating which states are observable or the string "All" if all states are observables. This can also be an experssion combining states (Only +,-,*,\ and ^ will be considered)
    pseudo_def["NoiseType"] = []; # String indicating if the desired noise is homoscedastic or heteroscedastic. The allowed strings are homo, homoscedastic or hetero, heteroscedastic (case independent). If empty, default will be heteroscedastic. Note that for high homoscedastic noise data could become negative (which to me does not make much sense), for now this will be set to a value close to 0.
    pseudo_def["Noise"] = []; # Percentage of heteroscedastic noise (introduced as value from 0 to 1). If empty 10% heteroscedastic noise will be assumed. In the homoscedastic case, if empty a standard deviation of 1 will be assumes. This has to be a vector of noise values for each observable.

    return(pseudo_def)

end

## Extract data from CSV
# Function to extract contents of CSV files if user sellects that
function extractPseudoDatCSV(model_def, pseudo_def)

    # Check taht all the dictionary entries are correct
    entries = ["ObservablesFile", "EventInputsFile", "theta", "MainDir", "plot", "flag", "Obs", "Noise"]
    if symdiff(entries,keys(pseudo_def))!=[]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the entries of the dictionary, there is something wrong...")
        println(symdiff(entries,keys(pseudo_def)))
        return
    end

    # Checks on the different strings introduced
    if pseudo_def["MainDir"] == [] || pseudo_def["MainDir"] == ""
        pseudo_def["MainDir"] = "";
    elseif typeof(pseudo_def["MainDir"])==Array{String,1}
        pseudo_def["MainDir"] = pseudo_def["MainDir"][1]
    elseif typeof(pseudo_def["MainDir"]) == String
        nothing
    else
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, there is something wrong in the type of the field MainDir")
        return
    end

    if typeof(pseudo_def["ObservablesFile"]) != Array{String,1}
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, there is something wrong in the type of the field ObservablesFile. This should be a Array{String,1}")
        return
    end

    if length(pseudo_def["ObservablesFile"]) != length(pseudo_def["EventInputsFile"])
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, there number of entries for ObservablesFile and EventInputsFile should be the same")
        return
    end

    for i in 1:length(pseudo_def["ObservablesFile"])
        if pseudo_def["ObservablesFile"][i][end-3:end] != ".csv"
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, add the .csv termination to the ObservablesFile file! And no spaces after!")
            return
        end
    end

    if typeof(pseudo_def["EventInputsFile"]) != Array{String,1}
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, there is something wrong in the type of the field EventInputsFile. This should be a Array{String,1}")
        return
    end

    for i in 1:length(pseudo_def["EventInputsFile"])
        if pseudo_def["EventInputsFile"][i][end-3:end] != ".csv"
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, add the .csv termination to the EventInputsFile file! And no spaces after!")
            return
        end
    end

    simstu = defPseudoDatStruct();

    simstu["theta"] = pseudo_def["theta"];

    tmpInp = Array{Any}(undef, 1,length(pseudo_def["EventInputsFile"]))
    for i in 1:length(pseudo_def["EventInputsFile"])
        if pseudo_def["MainDir"] != [] && pseudo_def["MainDir"] != ""
            tmpInp[i] = Matrix(CSV.read(string(pseudo_def["MainDir"], "\\\\",pseudo_def["EventInputsFile"][i])));
        else
            tmpInp[i] = Matrix(CSV.read(pseudo_def["EventInputsFile"][i]));
        end
    end

    tmpObs = Array{Any}(undef, 1,length(pseudo_def["ObservablesFile"]))
    for i in 1:length(pseudo_def["ObservablesFile"])
        if pseudo_def["MainDir"] != [] && pseudo_def["MainDir"] != ""
            tmpObs[i] = Matrix(CSV.read(string(pseudo_def["MainDir"], "\\\\",pseudo_def["ObservablesFile"][i])));
        else
            tmpObs[i] = Matrix(CSV.read(pseudo_def["ObservablesFile"][i]));
        end
    end

    try
        simstu["Nexp"] = length(pseudo_def["EventInputsFile"]);
        simstu["finalTime"] = [tmpInp[i][1,2] for i in 1:length(pseudo_def["EventInputsFile"])];
        simstu["switchT"] = [vcat(tmpInp[i][:,1], tmpInp[i][1,2]) for i in 1:length(pseudo_def["EventInputsFile"])];
        simstu["y0"] = [tmpObs[i][1,2:(1+model_def["nStat"])] for i in 1:length(pseudo_def["ObservablesFile"])];
        simstu["preInd"] = [tmpInp[i][1,3:(2+model_def["nInp"])] for i in 1:length(pseudo_def["EventInputsFile"])];
        simstu["uInd"] = [tmpInp[i][:,(3+model_def["nInp"]):(3+model_def["nInp"]+(model_def["nInp"]-1))] for i in 1:length(pseudo_def["EventInputsFile"])];
        simstu["theta"] = pseudo_def["theta"];
        simstu["tsamps"] = [tmpObs[i][:,1] for i in 1:length(pseudo_def["ObservablesFile"])];
        simstu["plot"] = pseudo_def["plot"];
        simstu["flag"] = pseudo_def["flag"];
        simstu["Obs"] = pseudo_def["Obs"];
        simstu["NoiseType"] = pseudo_def["NoiseType"];
        simstu["Noise"] = pseudo_def["Noise"];
    catch
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, but there seems to be some issue with the internal structure of the files you have given...")
        return
    end

    simstu = checkStructPseudoDat(model_def, simstu);

    return(simstu)

end

function PDatCSVGen(pseudo_res,model_def,pseudo_def)
    ## Generate CSV Files
    if !isdir(string(pseudo_def["savepath"], "\\PseudoDataFiles"))
        mkdir(string(pseudo_def["savepath"], "\\PseudoDataFiles"))
    end

    for expp in 1:pseudo_def["Nexp"]

        # Simulations Matrix
        a,b,c = size(pseudo_res["Sims"][string("Exp_", expp)]);
        Simulations = zeros(a, (b*c)+1);
        Simulations[:,1] = pseudo_def["tsamps"][expp];

        sims_head = ["time"];

        coun = 2;
        for i in 1:c
            Simulations[:,coun:(coun+b-1)] = pseudo_res["Sims"][string("Exp_", expp)][:,:,i];
            coun = coun+b;
            for j in 1:b
                sims_head = vcat(sims_head, string(model_def["stName"][j], "_Theta_", i))
            end
        end

        df1 = DataFrame[Simulations][1]
        rename!(df1, [sims_head[i] for i in 1:(b*c)+1])
        CSV.write(string(pseudo_def["savepath"], "\\PseudoDataFiles\\",model_def["NameF"],"_EXP",expp, "_",
                pseudo_def["flag"],"_Simulations.csv"), df1);

        # Observables matrix
        a,b,c = size(pseudo_res["PData"][string("PDExp_", expp)]);
        Observables = zeros(a, (3*b)*c);
        obs_head = Array{String,1}(undef,(3*b)*c);

        r = 1:3:(3*b)*c
        for i in 1:length(r)
            Observables[:,r[i]] = pseudo_def["tsamps"][expp];
            obs_head[r[i]] = string("time",i);
        end

        co = 2;
        for i in 1:c
            for j in 1:b
                Observables[:,co] = pseudo_res["PData"][string("PDExp_", expp)][:,j,i];
                Observables[:,co+1] = pseudo_res["PError"][string("PDExp_", expp)][:,j,i]
                obs_head[co] = string(pseudo_def["Obs"][j],"_theta", i, "_Mean");
                obs_head[co+1] = string(pseudo_def["Obs"][j],"_theta", i, "_Std");
                co +=3;
            end
        end

        df2 = DataFrame[Observables][1]
        rename!(df2, [obs_head[i] for i in 1:(3*b)*c])
        CSV.write(string(pseudo_def["savepath"], "\\PseudoDataFiles\\",model_def["NameF"],"_EXP",expp, "_",
                pseudo_def["flag"],"_Observables.csv"), df2);

        # Event Inputs matrix
        if model_def["nInp"] != 0
            EvnInputs = zeros(length(pseudo_def["uInd"][expp]),2+(model_def["nInp"]*2));
            evins_head = Array{String,1}(undef,2+(model_def["nInp"]*2));
            evins_head[1] = "Switchingtimes";
            evins_head[2] = "FinalTime";

            EvnInputs[:,1] = pseudo_def["switchT"][expp][1:end-1];
            EvnInputs[:,2] = repeat([pseudo_def["switchT"][expp][end]], outer = [length(pseudo_def["uInd"][expp])]);

            for i in 1:model_def["nInp"]
                if pseudo_def["preInd"] != []
                    EvnInputs[:,(2+i)] = repeat([pseudo_def["preInd"][expp][i]], outer = [length(pseudo_def["uInd"][expp])]);
                else
                    EvnInputs[:,(2+i)] = repeat([0], outer = [length(pseudo_def["uInd"][expp])]);
                end
                evins_head[(2+i)] = string(model_def["inpName"][i],"_Pre");
                EvnInputs[:,(2+model_def["nInp"])+i] = pseudo_def["uInd"][expp];
                evins_head[(2+model_def["nInp"])+i] = string(model_def["inpName"][i]);
            end

            df3 = DataFrame[EvnInputs][1]
            rename!(df3, [evins_head[i] for i in 1:2+(model_def["nInp"]*2)])
            CSV.write(string(pseudo_def["savepath"], "\\PseudoDataFiles\\",model_def["NameF"],"_EXP",expp, "_",
                    pseudo_def["flag"],"_Events_Inputs.csv"), df3);
        end

    end

end

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


text:
module AOC2021

const data_dir = normpath(joinpath(@__FILE__, "..", "..", "data"))

const PKGDIR = pkgdir(AOC2021)
src_path(day::AbstractString) = joinpath(PKGDIR, "src", "$day.jl")

const DAYS = ["day" * lpad(day, 2, "0") for day = 1:25]
const SRC_PATHS = let
    paths = src_path.(DAYS)
    existing_paths = filter(isfile, paths)
end

for path in SRC_PATHS
    include(path)
end

export solve

function solve(day::Int64, part::Int64)
    module_string = "Day" * lpad(day, 2, "0")
    part_string = string(part)
    module_and_part_string = module_string * ".part" * part_string * "()"
    eval(Meta.parse(module_and_part_string))
end

end # module

ext: jl


text:
## Generic ##

Tracker.dual(x::Bool, p) = x
Base.prevfloat(r::TrackedReal) = track(prevfloat, r)
@grad function prevfloat(r::Real)
    prevfloat(data(r)), Δ -> Δ
end
Base.nextfloat(r::TrackedReal) = track(nextfloat, r)
@grad function nextfloat(r::Real)
    nextfloat(data(r)), Δ -> Δ
end

for f = [:hcat, :vcat]
    for c = [
        [:TrackedReal],
        [:AbstractVecOrMat, :TrackedReal],
        [:TrackedVecOrMat, :TrackedReal],
    ]
        cnames = map(_ -> gensym(), c)
        @eval begin
            function Base.$f(
                $([:($x::$c) for (x, c) in zip(cnames, c)]...),
                x::Union{TrackedArray,TrackedReal},
                xs::Union{AbstractArray,Number}...,
            )
                return track($f, $(cnames...), x, xs...)
            end
        end
    end
    @eval begin
        @grad function $f(x::Real)
            $f(data(x)), (Δ) -> (Δ[1],)
        end
        @grad function $f(x1::Real, x2::Real)
            $f(data(x1), data(x2)), (Δ) -> (Δ[1], Δ[2])
        end
        @grad function $f(x1::AbstractVector, x2::Real)
            $f(data(x1), data(x2)), (Δ) -> (Δ[1:length(x1)], Δ[length(x1)+1])
        end
    end
end

function Base.copy(
    A::TrackedArray{T, 2, <:Adjoint{T, <:AbstractTriangular{T, <:AbstractMatrix{T}}}},
) where {T <: Real}
    return track(copy, A)
end
@grad function Base.copy(
    A::TrackedArray{T, 2, <:Adjoint{T, <:AbstractTriangular{T, <:AbstractMatrix{T}}}},
) where {T <: Real}
    return copy(data(A)), ∇ -> (copy(∇),)
end

Base.:*(A::TrackedMatrix, B::AbstractTriangular) = track(*, A, B)
Base.:*(A::AbstractTriangular{T}, B::TrackedVector) where {T} = track(*, A, B)
Base.:*(A::AbstractTriangular{T}, B::TrackedMatrix) where {T} = track(*, A, B)
Base.:*(A::Adjoint{T, <:AbstractTriangular{T}}, B::TrackedMatrix) where {T} = track(*, A, B)
Base.:*(A::Adjoint{T, <:AbstractTriangular{T}}, B::TrackedVector) where {T} = track(*, A, B)

function Base.fill(
    value::TrackedReal,
    dims::Vararg{Union{Integer, AbstractUnitRange}},
)
    return track(fill, value, dims...)
end
@grad function Base.fill(value::Real, dims...)
    return fill(data(value), dims...), function(Δ)
        size(Δ) ≢  dims && error("Dimension mismatch")
        return (sum(Δ), map(_->nothing, dims)...)
    end
end

## StatsFuns ##

logsumexp(x::TrackedArray; dims=:) = _logsumexp(x, dims)
_logsumexp(x::TrackedArray, dims=:) = track(_logsumexp, x, dims)
@grad function _logsumexp(x::TrackedArray, dims)
    lse = logsumexp(data(x), dims = dims)
    return lse, Δ -> (Δ .* exp.(x .- lse), nothing)
end

## Linear algebra ##

# Work around https://github.com/FluxML/Tracker.jl/pull/9#issuecomment-480051767

upper(A::AbstractMatrix) = UpperTriangular(A)
lower(A::AbstractMatrix) = LowerTriangular(A)
function upper(C::Cholesky)
    if C.uplo == 'U'
        return upper(C.factors)
    else
        return copy(lower(C.factors)')
    end
end
function lower(C::Cholesky)
    if C.uplo == 'U'
        return copy(upper(C.factors)')
    else
        return lower(C.factors)
    end
end

LinearAlgebra.LowerTriangular(A::TrackedMatrix) = lower(A)
lower(A::TrackedMatrix) = track(lower, A)
@grad lower(A) = lower(Tracker.data(A)), ∇ -> (lower(∇),)

LinearAlgebra.UpperTriangular(A::TrackedMatrix) = upper(A)
upper(A::TrackedMatrix) = track(upper, A)
@grad upper(A) = upper(Tracker.data(A)), ∇ -> (upper(∇),)

function LinearAlgebra.cholesky(A::TrackedMatrix; check=true)
    factors_info = turing_chol(A, check)
    factors = factors_info[1]
    info = data(factors_info[2])
    return Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
end
function LinearAlgebra.cholesky(A::Symmetric{<:Any, <:TrackedMatrix}; check=true)
    uplo = A.uplo == 'U' ? (:U) : (:L)
    factors_info = symm_turing_chol(parent(A), check, uplo)
    factors = factors_info[1]
    info = data(factors_info[2])
    return Cholesky{eltype(factors), typeof(factors)}(factors, 'U', info)
end

function turing_chol(A::AbstractMatrix, check)
    chol = cholesky(A, check=check)
    (chol.factors, chol.info)
end
function symm_turing_chol(A::AbstractMatrix, check, uplo)
    chol = cholesky(Symmetric(A, uplo), check=check)
    (chol.factors, chol.info)
end

turing_chol(A::TrackedMatrix, check) = track(turing_chol, A, check)
@grad function turing_chol(A::AbstractMatrix, check)
    C, back = ZygoteRules.pullback(data(A), check) do A, check
        cholesky(A, check=check)
    end
    return (C.factors, C.info), Δ->back((factors=data(Δ[1]),))
end

symm_turing_chol(A::TrackedMatrix, check, uplo) = track(symm_turing_chol, A, check, uplo)
@grad function symm_turing_chol(A::AbstractMatrix, check, uplo)
    C, back = ZygoteRules.pullback(data(A), check, uplo) do A, check, uplo
        cholesky(Symmetric(A, uplo), check=check)
    end
    return (C.factors, C.info), Δ->back((factors=data(Δ[1]),))
end

# Specialised logdet for cholesky to target the triangle directly.
logdet_chol_tri(U::AbstractMatrix) = 2 * sum(log, U[diagind(U)])
logdet_chol_tri(U::TrackedMatrix) = track(logdet_chol_tri, U)
@grad function logdet_chol_tri(U::AbstractMatrix)
    U_data = data(U)
    return logdet_chol_tri(U_data), Δ->(Matrix(Diagonal(2 .* Δ ./ diag(U_data))),)
end

function LinearAlgebra.logdet(C::Cholesky{<:TrackedReal, <:TrackedMatrix})
    return logdet_chol_tri(C.U)
end

# Tracker's implementation of ldiv isn't good. We'll use Zygote's instead.

zygote_ldiv(A::AbstractMatrix, B::AbstractVecOrMat) = A \ B
function zygote_ldiv(A::TrackedMatrix, B::TrackedVecOrMat)
    return track(zygote_ldiv, A, B)
end
function zygote_ldiv(A::TrackedMatrix, B::AbstractVecOrMat)
    return track(zygote_ldiv, A, B)
end
zygote_ldiv(A::AbstractMatrix, B::TrackedVecOrMat) =  track(zygote_ldiv, A, B)
@grad function zygote_ldiv(A, B)
    Y, back = ZygoteRules.pullback(\, data(A), data(B))
    return Y, Δ->back(data(Δ))
end

function Base.:\(a::Cholesky{<:TrackedReal, <:TrackedArray}, b::AbstractVecOrMat)
    return (a.U \ (a.U' \ b))
end

# SpecialFunctions

SpecialFunctions.logabsgamma(x::TrackedReal) = track(logabsgamma, x)
@grad function SpecialFunctions.logabsgamma(x::Real)
    return logabsgamma(data(x)), Δ -> (digamma(data(x)) * Δ[1],)
end
ZygoteRules.@adjoint function SpecialFunctions.logabsgamma(x::Real)
    return logabsgamma(x), Δ -> (digamma(x) * Δ[1],)
end

# Zygote fill has issues with non-numbers

ZygoteRules.@adjoint function fill(x::T, dims...) where {T}
    return ZygoteRules.pullback(x, dims...) do x, dims...
        return reshape([x for i in 1:prod(dims)], dims)
    end
end

# isprobvec

function Distributions.isprobvec(p::TrackedArray{<:Real})
    pdata = Tracker.data(p)
    all(x -> x ≥ zero(x), pdata) && isapprox(sum(pdata), one(eltype(pdata)), atol = 1e-6)
end

# Some array functions - workaround https://github.com/FluxML/Tracker.jl/issues/4

import Base: +, -, *, /, \
import LinearAlgebra: dot

for f in (:+, :-, :*, :/, :\, :dot), (T1, T2) in [
    (:TrackedArray, :AbstractArray),
    (:TrackedMatrix, :AbstractMatrix),
    (:TrackedMatrix, :AbstractVector),
    (:TrackedVector, :AbstractMatrix),
]
    @eval begin
        function $f(a::$T1{T}, b::$T2{<:TrackedReal}) where {T <: Real}
            return $f(convert(AbstractArray{TrackedReal{T}}, a), b)
        end
        function $f(a::$T2{<:TrackedReal}, b::$T1{T}) where {T <: Real}
            return $f(a, convert(AbstractArray{TrackedReal{T}}, b))
        end
    end
end


text:
# ---
# title: 883. Projection Area of 3D Shapes
# id: problem883
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/projection-area-of-3d-shapes/description/>
# hidden: true
# ---
#
# On a `N * N` grid, we place some `1 * 1 * 1 `cubes that are axis-aligned with
# the x, y, and z axes.
#
# Each value `v = grid[i][j]` represents a tower of `v` cubes placed on top of
# grid cell `(i, j)`.
#
# Now we view the  _projection_  of these cubes onto the xy, yz, and zx planes.
#
# A projection is like a shadow, that maps our 3 dimensional figure to a 2
# dimensional plane.
#
# Here, we are viewing the "shadow" when looking at the cubes from the top, the
# front, and the side.
#
# Return the total area of all three projections.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[2]]
#     Output: 5
#
#
# **Example 2:**
#
#
#
#     Input: [[1,2],[3,4]]
#     Output: 17
#     Explanation:
#     Here are the three projections ("shadows") of the shape made with each axis-aligned plane.
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/08/02/shadow.png)
#
#
# **Example 3:**
#
#
#
#     Input: [[1,0],[0,2]]
#     Output: 8
#
#
# **Example 4:**
#
#
#
#     Input: [[1,1,1],[1,0,1],[1,1,1]]
#     Output: 14
#
#
# **Example 5:**
#
#
#
#     Input: [[2,2,2],[2,1,2],[2,2,2]]
#     Output: 21
#
#
#
#
# **Note:**
#
#   * `1 <= grid.length = grid[0].length <= 50`
#   * `0 <= grid[i][j] <= 50`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
## These are the figures that appear in the paper
# "Spectra of Jacobi operators via connection coefficients" by Marcus Webb and Sheehan Olver (2018)

using ApproxFun, SpectralMeasures, Plots; Plots.gr(linewidth=3,legend=false)
using ComplexPhasePortrait, PyPlot, FileIO; pygui(false) # PyPlot is good for Wegert Plots

x = linspace(-3, 3, 300)
  Z = x .+ flipdim(x,1)'*im

mkdir("Figs")

# Example for nondecreasingness of resolvents in intervals of continuity
MyResolvent = x -> real(1 + 2*sqrt(complex(x-1.)).*sqrt(complex(x+1.))-2*x + x./(9-x.^2))
  Plots.plot([-1,1],[0,0],color=:black)
  plot!([-3,-3,-3],[-3,0,3],line=:dash,color=:black,marker=:circle,markersize=2,linewidth=1)
  plot!([3,3,3],[-3,0,3],line=:dash,color=:black,marker=:circle,markersize=2,linewidth=1)
  plot!([-1,-1],[-3,3],line=:dash,color=:black,linewidth=1)
  plot!([1,1],[-3,3],line=:dash,color=:black,linewidth=1)
  plot!(-8:0.001:8,[MyResolvent(-8:0.001:-1);NaN*(-1.001:0.001:0.999);MyResolvent(1:0.001:8)],xlabel="\\lambda",ylabel="f\\(\\lambda\\)",xlims=(-7.5,7.5),ylims=(-2,2))
  Plots.savefig("Figs/secular.pdf")

# Pertrubation of the [1,1] element
for a in [0, 0.15, 0.35, 0.5, 0.75, 1.]
  Plots.plot(spectralmeasure([a],[.5]),xlims=(-2,2),ylims=(0,2),title="\\alpha=$(a)")
  Plots.savefig("Figs/spectralmeasurealpha=$(round(Int,100*a)).pdf")

  R = principalresolvent([a],[.5])
  save("Figs/principalresolventalpha=$(round(Int,100*a)).png", portrait(R(Z),PTstepmod))
  imshow(imread("Figs/principalresolventalpha=$(round(Int,100*a)).png"),extent=[-3,3,-3,3])
  PyPlot.plot([-1.,1.],[0.,0.],linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/principalresolventalpha=$(round(Int,100*a)).png",transparent=true,bbox_inches="tight")
  clf()

  r = discresolvent([a],[.5])
  save("Figs/discresolventalpha=$(round(Int,100*a)).png", portrait(r(Z),PTstepmod))
  imshow(imread("Figs/discresolventalpha=$(round(Int,100*a)).png"),extent=[-3,3,-3,3])
  PyPlot.plot(cos(linspace(0,2pi,100)),sin(linspace(0,2pi,100)),linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/discresolventalpha=$(round(Int,100*a)).png",transparent=true,bbox_inches="tight")
  clf()
end

# Perturbation of the [1,2] element
for b in [.5, 0.707, .85, 1, 1.2, 1.5]
  Plots.plot(spectralmeasure([0.],[b]/sqrt(2)),xlims=(-2,2),ylims=(0,2),title="\\beta=$(b)")
  Plots.savefig("Figs/spectralmeasurebeta=$(round(Int,1000*b)).pdf")

  R = principalresolvent([0.],[b]/sqrt(2))
  save("Figs/principalresolventbeta=$(round(Int,1000*b)).png", portrait(R(Z),PTstepmod))
  imshow(imread("Figs/principalresolventbeta=$(round(Int,1000*b)).png"),extent=[-3,3,-3,3])
  PyPlot.plot([-1.,1.],[0.,0.],linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/principalresolventbeta=$(round(Int,1000*b)).png",transparent=true,bbox_inches="tight")
  clf()

  r = discresolvent([0.],[b]/sqrt(2))
  save("Figs/discresolventbeta=$(round(Int,1000*b)).png", portrait(r(Z),PTstepmod))
  imshow(imread("Figs/discresolventbeta=$(round(Int,1000*b)).png"),extent=[-3,3,-3,3])
  PyPlot.plot(cos(linspace(0,2pi,100)),sin(linspace(0,2pi,100)),linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/discresolventbeta=$(round(Int,1000*b)).png",transparent=true,bbox_inches="tight")
  clf()
end

# 3 by 3 Perturbation
a = [3/4,-1/4,1/2]; b=[1,3/4]
  μ = spectralmeasure(a,b)
  R = principalresolvent(a,b)
  r = discresolvent(a,b)

  Plots.plot(μ,xlims=(-3,3),ylims=(0,2))
  Plots.savefig("Figs/spectralmeasure3by3pert.pdf")

  save("Figs/principalresolvent3by3pert.png", portrait(R(Z),PTstepmod))
  imshow(imread("Figs/principalresolvent3by3pert.png"),extent=[-3,3,-3,3])
  PyPlot.plot([-1.,1.],[0.,0.],linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/principalresolvent3by3pert.png",transparent=true,bbox_inches="tight")
  clf()

  save("Figs/discresolvent3by3pert.png", portrait(r(Z),PTstepmod))
  imshow(imread("Figs/discresolvent3by3pert.png"),extent=[-3,3,-3,3])
  PyPlot.plot(cos(linspace(0,2pi,100)),sin(linspace(0,2pi,100)),linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/discresolvent3by3pert.png",transparent=true,bbox_inches="tight")
  clf()

# Legendre
for n in [1,2,3,10,30,100]
  bLeg = (1:n-1)./sqrt(4*(1:n-1).^2-1)
  Plots.plot(spectralmeasure(zeros(n),bLeg),title="Legendre approximation, n = $(n)",xlims=(-1.5,1.5),ylims=(0,1),yticks=[0,.25,.5,.75,1])
  Plots.savefig("Figs/spectralmeasureLegn=$(n).pdf")

  R = principalresolvent([0.],bLeg)
  save("Figs/principalresolventLegn=$(n).png", portrait(R(Z),PTstepmod))
  imshow(imread("Figs/principalresolventLegn=$(n).png"),extent=[-3,3,-3,3])
  PyPlot.plot([-1.,1.],[0.,0.],linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/principalresolventLegn=$(n).png",transparent=true,bbox_inches="tight")
  clf()

  r = discresolvent([0.],bLeg)
  save("Figs/discresolventLegn=$(n).png", portrait(r(Z),PTstepmod))
  imshow(imread("Figs/discresolventLegn=$(n).png"),extent=[-3,3,-3,3])
  PyPlot.plot(cos(linspace(0,2pi,100)),sin(linspace(0,2pi,100)),linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/discresolventLegn=$(n).png",transparent=true,bbox_inches="tight")
  clf()
end

# Ultraspherical polynomials
γ = 0.6
  for n in [1,2,3,10,30,100]
  bUlt = .5*sqrt(((1:n).*(2γ+(0:n-1)))./((γ+(0:n-1)).*(γ+(1:n))))
  Plots.plot(spectralmeasure(zeros(n),bUlt),title="Ultraspherical\\($(γ)\\), n = $(n)",xlims=(-1.5,1.5),ylims=(0,1),yticks=[0,.25,.5,.75,1])
  Plots.savefig("Figs/spectralmeasureUltra0pt6n=$(n).pdf")

  R = principalresolvent([0.],bUlt)
  save("Figs/principalresolventUltra0pt6n=$(n).png", portrait(R(Z),PTstepmod))
  imshow(imread("Figs/principalresolventUltra0pt6n=$(n).png"),extent=[-3,3,-3,3])
  PyPlot.plot([-1.,1.],[0.,0.],linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/principalresolventUltra0pt6n=$(n).png",transparent=true,bbox_inches="tight")
  clf()

  r = discresolvent([0.],bUlt)
  save("Figs/discresolventUltra0pt6n=$(n).png", portrait(r(Z),PTstepmod))
  imshow(imread("Figs/discresolventUltra0pt6n=$(n).png"),extent=[-3,3,-3,3])
  PyPlot.plot(cos(linspace(0,2pi,100)),sin(linspace(0,2pi,100)),linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/discresolventUltra0pt6n=$(n).png",transparent=true,bbox_inches="tight")
  clf()
end

# Jacobi polynomials
α = 0.4; β = 1.9
  for n in [1,2,3,10,30,100]
  aJac = (β.^2-α.^2)./((2.*(0:n)+α+β).*(2.*(1:n+1)+α+β))
  bJac = 2*sqrt(((1:n).*(α+(1:n)).*(β+(1:n)).*(α+β+(1:n)))./((2.*(1:n)+α+β-1).*((2.*(1:n)+α+β).^2).*(2.*(1:n)+α+β+1)))
  Plots.plot(spectralmeasure(aJac,bJac),title="Jacobi\\($(α),$(β)\\), n = $(n)",xlims=(-1.5,1.5),ylims=(0,1),yticks=[0,.25,.5,.75,1])
  Plots.savefig("Figs/spectralmeasureJacobin=$(n).pdf")

  R = principalresolvent(aJac,bJac)
  save("Figs/principalresolventJacobin=$(n).png", portrait(R(Z),PTstepmod))
  imshow(imread("Figs/principalresolventJacobin=$(n).png"),extent=[-3,3,-3,3])
  PyPlot.plot([-1.,1.],[0.,0.],linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/principalresolventJacobin=$(n).png",transparent=true,bbox_inches="tight")
  clf()

  r = discresolvent(aJac,bJac)
  save("Figs/discresolventJacobin=$(n).png", portrait(r(Z),PTstepmod))
  imshow(imread("Figs/discresolventJacobin=$(n).png"),extent=[-3,3,-3,3])
  PyPlot.plot(cos(linspace(0,2pi,100)),sin(linspace(0,2pi,100)),linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/discresolventJacobin=$(n).png",transparent=true,bbox_inches="tight")
  clf()
end

# Random
mkdir("Figs")
srand(200)
longaRand = 3*(2rand(100)-1)./(1:100).^2
for n in [1,2,3,10,50,100]
  aRand = longaRand[1:n]
  Plots.plot(spectralmeasure(aRand,Float64[]),title="Random, n = $(n)",xlims=(-1.5,1.5),ylims=(0,1),yticks=[0,.25,.5,.75,1])
  Plots.savefig("Figs/spectralmeasurerandomn=$(n).pdf")

  R = principalresolvent(aRand,Float64[])
  save("Figs/principalresolventrandomn=$(n).png", portrait(R(Z),PTstepmod))
  imshow(imread("Figs/principalresolventrandomn=$(n).png"),extent=[-3,3,-3,3])
  PyPlot.plot([-1.,1.],[0.,0.],linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/principalresolventrandomn=$(n).png",transparent=true,bbox_inches="tight")
  clf()

  r = discresolvent(aRand,Float64[])
  save("Figs/discresolventrandomn=$(n).png", portrait(r(Z),PTstepmod))
  imshow(imread("Figs/discresolventrandomn=$(n).png"),extent=[-3,3,-3,3])
  PyPlot.plot(cos(linspace(0,2pi,100)),sin(linspace(0,2pi,100)),linewidth=2,color=:black)
  PyPlot.grid("on")
  PyPlot.savefig("Figs/discresolventrandomn=$(n).png",transparent=true,bbox_inches="tight")
  clf()
end

a = zeros(100); a[1] = 1
  Plots.plot(spectralmeasure(a,Float64[]),title="\\alpha=[1,0,...]",xlims=(-1,1.5),ylims=(0,4.2),yticks=[0,1,2,3,4])
  Plots.savefig("Figs/spectralmeasuredelta1.pdf")

a = zeros(100); a[2] = 1
  Plots.plot(spectralmeasure(a,Float64[]),title="\\alpha=[0,1,0,...]",xlims=(-1,1.5),ylims=(0,4.2),yticks=[0,1,2,3,4])
    Plots.savefig("Figs/spectralmeasuredelta2.pdf")

a = zeros(100); a[5] = 1
  Plots.plot(spectralmeasure(a,Float64[]),title="\\alpha=[0,0,0,0,1,0,...]",xlims=(-1,1.5),ylims=(0,4.2),yticks=[0,1,2,3,4])
  Plots.savefig("Figs/spectralmeasuredelta5.pdf")


text:
@testset "AffineScalingLayer" begin
    include("affineScalingLayerTest.jl")
end

@testset "DoubleSymLayer" begin
    include("doubleSymLayerTest.jl")
end

@testset "NormLayer" begin
    include("normLayerTest.jl")
end

@testset "SingleLayer" begin
    include("singleLayerTest.jl")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 04af0d10-9070-11eb-2fbf-493475b83bd0
using Plots

# ╔═╡ bca6bf50-915e-11eb-11e8-1d2fc21da187
using CSV;

# ╔═╡ 1b7aa6b0-906d-11eb-1bcd-9fd48415fade
md"""
# Solving for resistance

The problem is setup as a rectangle with width 2L and height 2b. In the center there is a hole with radius a. The conductivtiy of the material is $\sigma$. a will be varied.

$ R = \frac{b}{\sigma L}\frac{L^2+a^2}{b^2-a^2} $
"""

# ╔═╡ a6e8f7a0-906e-11eb-1a0e-51f4a4b751f7
σ = 5.998e7;

# ╔═╡ d3572c80-906e-11eb-0a39-73a83bffb283
b = 30;

# ╔═╡ d9e9e9c0-906e-11eb-2507-2bf61260bf5a
L = 50;

# ╔═╡ 94eb35b0-9161-11eb-1c77-5d337997146a
simulationData = CSV.File("./data/resistance2.csv")

# ╔═╡ de4fdbf2-906e-11eb-2715-6f5bc8d68d2a
a = simulationData.a; #Use same values of a as in the simulation

# ╔═╡ 03943b70-9166-11eb-3f37-6fe974598a65
simulationData.R

# ╔═╡ 662d3360-906f-11eb-2826-d18803f2c32e
R = (b*(L^2 .+ a.^2))./(σ*L*(b^2 .- a.^2))

# ╔═╡ 3c5e1dee-9070-11eb-1369-6509d6df4f3c
plot(a, R, xlabel="a [m]", ylabel="Resistance [Ohms]", label="Analytical R")

# ╔═╡ 6fb114e0-9166-11eb-36af-a5a5faa376e1
plot(a, simulationData.R, xlabel="a [m]", ylabel="Resistance [Ohms]", label="Simulated R")

# ╔═╡ 5ff3b170-9166-11eb-25e0-1f2818ee7004
plot(a, [R simulationData.R], xlabel="a [m]", ylabel="Resistance [Ohms] (log)", label=["Analytical R" "Simulated R"], yaxis=:log)

# ╔═╡ 09dae670-9164-11eb-1ef0-f14e5b401a2d
ϵ = abs.((simulationData.R - R) ./ (simulationData.R + R))

# ╔═╡ 89572940-9164-11eb-1669-41593429cd22
plot(a, ϵ, xlabel="a [m]", ylabel="Difference in result", label="Epsilon")

# ╔═╡ Cell order:
# ╟─1b7aa6b0-906d-11eb-1bcd-9fd48415fade
# ╠═04af0d10-9070-11eb-2fbf-493475b83bd0
# ╠═bca6bf50-915e-11eb-11e8-1d2fc21da187
# ╠═a6e8f7a0-906e-11eb-1a0e-51f4a4b751f7
# ╠═d3572c80-906e-11eb-0a39-73a83bffb283
# ╠═d9e9e9c0-906e-11eb-2507-2bf61260bf5a
# ╠═94eb35b0-9161-11eb-1c77-5d337997146a
# ╠═de4fdbf2-906e-11eb-2715-6f5bc8d68d2a
# ╠═03943b70-9166-11eb-3f37-6fe974598a65
# ╠═662d3360-906f-11eb-2826-d18803f2c32e
# ╠═3c5e1dee-9070-11eb-1369-6509d6df4f3c
# ╠═6fb114e0-9166-11eb-36af-a5a5faa376e1
# ╠═5ff3b170-9166-11eb-25e0-1f2818ee7004
# ╠═09dae670-9164-11eb-1ef0-f14e5b401a2d
# ╠═89572940-9164-11eb-1669-41593429cd22

max_stars_repo_name: tskarhed/julia-notebooks

max_issues_repo_name: tskarhed/julia-notebooks

max_forks_repo_name: tskarhed/julia-notebooks


text:
module MultiStartLH

using Dates, DocStringExtensions, FileIO, JLD2, Random

export test_dir;

include("logging.jl");
include("history.jl");
include("update_rule.jl");
include("multistart.jl");

test_dir() = normpath(joinpath(@__DIR__, "..", "testfiles"));

end # module


text:
@testset "Properties" begin

    atoms = readPDB("./structure.pdb", "protein")

    @test mass(atoms) ≈ 11079.704440000156
    f = formula(select(atoms, "residue < 5"))
    @test f.formula == [("H", 25), ("C", 19), ("N", 4), ("O", 7), ("S", 1)]

    @test element.(select(atoms, "residue = 1")) ==
          ["N", "H", "H", "H", "C", "H", "C", "H", "H", "H", "C", "O"]
    @test element_name.(select(atoms, "residue = 1")) == [
        "Nitrogen",
        "Hydrogen",
        "Hydrogen",
        "Hydrogen",
        "Carbon",
        "Hydrogen",
        "Carbon",
        "Hydrogen",
        "Hydrogen",
        "Hydrogen",
        "Carbon",
        "Oxygen",
    ]
    @test atomic_number.(select(atoms, "residue = 1")) ==
          [7, 1, 1, 1, 6, 1, 6, 1, 1, 1, 6, 8]
    residues = collect(eachresidue(atoms))
    @test length(residues) == 104
    @test name(Residue(atoms, 1:12)) == "ALA"
    @test Residue(atoms, 1:12).range == 1:12

    m = maxmin(atoms)
    @test m.xmin ≈ [-14.18, -17.561, -15.369]
    @test m.xmax ≈ [18.694, 14.182, 15.909]
    @test m.xlength ≈ [32.873999999999995, 31.743000000000002, 31.278]

    s = select(atoms, "residue = 3")
    @test coor(s) ≈ [
        SVector{3,Float64}(-4.383, -11.903, -6.849),
        SVector{3,Float64}(-4.51, -11.263, -6.096),
        SVector{3,Float64}(-3.903, -11.262, -8.062),
        SVector{3,Float64}(-3.731, -12.076, -8.767),
        SVector{3,Float64}(-4.938, -10.279, -8.612),
        SVector{3,Float64}(-4.417, -9.552, -9.06),
        SVector{3,Float64}(-5.543, -9.911, -7.784),
        SVector{3,Float64}(-5.867, -10.85, -9.684),
        SVector{3,Float64}(-5.451, -10.837, -10.863),
        SVector{3,Float64}(-6.974, -11.289, -9.3),
        SVector{3,Float64}(-2.626, -10.48, -7.749),
        SVector{3,Float64}(-1.94, -10.014, -8.658)
    ]

    s2 = select(atoms, "residue = 5")
    @test distance(s, s2) ≈ 3.6750402718881863
    x1 = coor(s)
    x2 = coor(s2)
    @test distance(x1, x2) ≈ 3.6750402718881863
    @test distance(residues[3], residues[5]) ≈ 3.6750402718881863

    r = Residue(select(atoms, "residue = 3"))
    @test coor(s) == coor(r)
    @test coor(select(atoms, "residue = 3")) == coor(residues[3])
    @test PDBTools.same_residue(atoms[1], atoms[2]) == true
    @test PDBTools.same_residue(atoms[1], atoms[20]) == false

    #
    # Dispatch of closest and distance functions
    #
    r1 = select(atoms, "residue = 3")
    r2 = select(atoms, "residue = 5")

    @test all(closest(r1,r2) .≈ (11, 2, 3.6750402718881863))
    @test all(closest(coor(r1),r2) .≈ (11, 2, 3.6750402718881863))
    @test all(closest(r1,coor(r2)) .≈ (11, 2, 3.6750402718881863))
    @test all(closest(coor(r1),coor(r2)) .≈ (11, 2, 3.6750402718881863))

    @test all(closest(r1[1],coor(r2)) .≈ (1, 2, 5.121218702613667))
    @test all(closest(coor(r1[1]),coor(r2)) .≈ (1, 2, 5.121218702613667))
    @test all(closest(coor(r1[1]),r2) .≈ (1, 2, 5.121218702613667))
    @test all(closest(coor(r1[1]),coor(r2[2])) .≈ (1, 1, 5.121218702613667))

    @test all(closest(r1[1],coor(r2[2])) .≈ (1, 1, 5.121218702613667))
    @test all(closest(coor(r1[1]),r2[2]) .≈ (1, 1, 5.121218702613667))
    @test all(closest(r1[1],r2[2]) .≈ (1, 1, 5.121218702613667))

    @test distance(r1,r2) ≈ 3.6750402718881863
    @test distance(coor(r1),r2) ≈ 3.6750402718881863
    @test distance(r1,coor(r2)) ≈ 3.6750402718881863
    @test distance(coor(r1),coor(r2)) ≈ 3.6750402718881863

    @test distance(r1[1],coor(r2)) ≈ 5.121218702613667
    @test distance(coor(r1[1]),coor(r2)) ≈ 5.121218702613667
    @test distance(coor(r1[1]),r2) ≈ 5.121218702613667
    @test distance(coor(r1[1]),coor(r2[2])) ≈ 5.121218702613667

    @test distance(r1[1],coor(r2[2])) ≈ 5.121218702613667
    @test distance(coor(r1[1]),r2[2]) ≈ 5.121218702613667
    @test distance(r1[1],r2[2]) ≈ 5.121218702613667

    #
    # Resiude name functions
    #
    @test residuename("glu") == "Glutamic acid"
    @test oneletter("Glu") == "E"
    @test threeletter("E") == "GLU"

    #
    # Mass of sequences
    #
    seq = "AEG"
    @test mass(Sequence(seq)) ≈ 257.2432
    seq = ["A", "E", "G"]
    @test mass(Sequence(seq)) ≈ 257.2432
    seq = ['A', 'E', 'G']
    @test mass(Sequence(seq)) ≈ 257.2432
    seq = ["ALA", "GLU", "GLY"]
    @test mass(Sequence(seq)) ≈ 257.2432
    seq = ["Alanine", "Glutamic acid", "Glycine"]
    @test mass(Sequence(seq)) ≈ 257.2432

end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
import Distributions: rand, logpdf, pdf, estimate

function sample_chain(chain)
    parms = (Symbol.(chain.name_map.parameters)...,)
    n = size(chain, 1)*size(chain, 3)
    idx = rand(1:n)
    vals = map(x -> chain[x][idx], parms)
    return NamedTuple{parms}(vals)
end

function _posterior_predictive(m, chain, f=x -> x)
    parms = sample_chain(chain)
    return f(m(parms))
end

"""
**posterior_predictive**

Returns posterior predictive distribution and optionally applies function to samples
    on each replication
- `model`: the data generating function of a model
- `chain`: an MCMCChains chain object
- `n_samples`: the number of samples
- `f`: a function that is applied to each sample from posterior predictive distribution
Function Signature
````julia
posterior_predictive(model, chain, n_samples::Int, f=x -> x)
````
"""
function posterior_predictive(model, chain, n_samples::Int, f=x -> x)
    return map(x -> _posterior_predictive(model, chain, f), 1:n_samples)
end

find_index(actr::AbstractACTR; check_value=true, criteria...) = find_index(actr.declarative.memory; check_value, criteria...)

"""
**find_index**

Returns the index of first chunk that matches a set of criteria
- `chunks`: an array of chunks
- `criteria`: a set of keyword arguments for slot-value pairs

Function Signature
````julia
find_index(chunks::Array{<:Chunk,1}; criteria...)
````
## Example
````julia
chunks = [Chunk(animal=:dog), Chunk(animal=cat)]
find_index(chunks; animal=:dog)
````
"""
function find_index(chunks::Array{<:Chunk,1}; check_value=true, criteria...)
    for (i,c) in enumerate(chunks)
        match(c; check_value, criteria...) ? (return i) : nothing
    end
    return -100
end

"""
**find_index**

Returns the index of first chunk that matches a set of criteria
- `actr`: ACTR object
- `funs`: a set of functions
- `criteria`: a set of keyword arguments for slot-value pairs

Function Signature
````julia
find_index(actr::ACTR, funs...; criteria...)
````
## Example
````julia
chunks = [Chunk(animal=:dog), Chunk(animal=cat)]
find_index(chunks; animal=:dog)
````
"""
find_index(actr::ACTR, funs...; check_value=true, criteria...) = find_index(actr.declarative.memory, funs...; check_value, criteria...)

"""
**find_index**

Returns the index of first chunk that matches a set of criteria
- `chunks`: an array of chunks
- `funs`: a set of functions
- `criteria`: a set of keyword arguments for slot-value pairs

Function Signature
````julia
find_index(chunks::Array{<:Chunk,1}; criteria...)
````
## Example
````julia
chunks = [Chunk(animal=:dog), Chunk(animal=cat)]
find_index(chunks; animal=:dog)
````
"""
function find_index(chunks::Array{<:Chunk,1}, funs...; check_value=true, criteria...)
    for (i,c) in enumerate(chunks)
        match(c, funs...; check_value, criteria...) ? (return i) : nothing
    end
    return -100
end

"""
**find_indices**

Returns the index of first chunk that matches a set of criteria
- `actr`: an ACTR object
- `criteria`: a set of keyword arguments for slot-value pairs

Function Signature
````julia
find_indices(actr::ACTR; criteria...)
````
## Example
````julia
chunks = [Chunk(animal=:dog), Chunk(animal=:dog), Chunk(animal=cat)]
find_indices(chunks; animal=:dog)
````
"""
find_indices(actr::ACTR; check_value = true, criteria...) = find_indices(actr.declarative.memory; check_value, criteria...)

"""
**find_indices**

Returns the index of first chunk that matches a set of criteria
- `chunks`: an array of chunks
- `criteria`: a set of keyword arguments for slot-value pairs

**Function Signature**
````julia
find_indices(actr::ACTR; criteria...)
````
**Example**
````julia
chunks = [Chunk(animal=:dog), Chunk(animal=:dog), Chunk(animal=cat)]
find_indices(chunks; animal=:dog)
````
"""
function find_indices(chunks::Array{<:Chunk,1}; check_value=true, criteria...)
    idx = Int64[]
    for (i,c) in enumerate(chunks)
        match(c; check_value, criteria...) ? push!(idx, i) : nothing
    end
    return idx
end

"""
**find_indices**

Returns the index of first chunk that matches a set of criteria
- `actr`: an ACTR object
- `funs`: a set of functions
- `criteria`: a set of keyword arguments for slot-value pairs

**Function Signature**
````julia
find_indices(actr::ACTR; criteria...)
````
**Example**
````julia
chunks = [Chunk(animal=:dog), Chunk(animal=:dog), Chunk(animal=cat)]
find_indices(chunks; animal=:dog)
````
"""
find_indices(actr::ACTR, funs...; check_value=true, criteria...) = find_indices(actr.declarative.memory, funs...; check_value, criteria...)

"""
**find_indices**

Returns the index of first chunk that matches a set of criteria
- `chunks`: an array of chunks
- `funs`: a set of functions
- `criteria`: a set of keyword arguments for slot-value pairs

**Function Signature**
````julia
find_indices(actr::ACTR; criteria...)
````
** Example**
````julia
chunks = [Chunk(animal=:dog), Chunk(animal=:dog), Chunk(animal=cat)]
find_indices(chunks; animal=:dog)
````
"""
function find_indices(chunks::Array{<:Chunk,1}, funs...; check_value=true, criteria...)
    idx = Int64[]
    for (i,c) in enumerate(chunks)
        match(c, funs...; check_value, criteria...) ? push!(idx, i) : nothing
    end
    return idx
end

"""
**get_iconic_memory**

Returns array of chunks or visual objects representing iconic memory
- `actr`: an ACTR object

**Function Signature**
````julia
get_iconic_memory(actr)
````
"""
get_iconic_memory(actr) = actr.visual_location.iconic_memory

"""
**get_visicon**

Returns array of chunks or visual objects representing all visual objects
within the simulation

**Function Signature**
````julia
get_visicon(actr)
````
"""
get_visicon(actr) = actr.visicon

"""
**get_buffer**

Returns posterior predictive distribution and optionally applies function to samples
    on each replication
- `actr`: an ACTR object
- `m`: name of module as a symbol

**Function Signature**
````julia
get_buffer(actr, m)
````

**Example**
````julia
get_buffer(actr, :imaginal)
````
"""
get_buffer(actr, m) = getfield(actr, m).buffer
"""
**get_buffer**

Returns posterior predictive distribution and optionally applies function to samples
    on each replication
- `actr`: an ACTR object
- `m`: name of module as a symbol
- `v`: new value

**Function Signature**
````julia
set_buffer!(actr, m, v)
````

Example
````julia
set_buffer!(actr, :imaginal, [chunk])
````
"""
set_buffer!(actr, m, v) = getfield(actr, m).buffer = v

"""
**get_mean_activations**

Returns a vector of mean activations i.e. activation without noise
- `actr`: an ACTR object

**Function Signature**
````julia
get_mean_activations(actr)
````
"""
get_mean_activations(actr::AbstractACTR) = get_mean_activations(actr.declarative)

"""
**get_mean_activations**

Returns a vector of mean activations i.e. activation without noise
- `declarative`: an Declarative memory object

**Function Signature**
````julia
get_mean_activations(memory)
````
"""
get_mean_activations(declarative::Declarative) = get_mean_activations(declarative.memory)

"""
    get_mean_activations(chunks::Vector{<:Chunk})

Returns a vector of mean activations i.e. activation without noise

# Arguments

- `chunks`: an array of chunks
"""
function get_mean_activations(chunks::Vector{<:Chunk})
    return map(x->x.act - x.act_noise, chunks)
end

"""
    get_name(actr::AbstractACTR)

Returns the name of an ACT-R model object.

- `actr`: an ACTR model object
"""
function get_name(actr::AbstractACTR)
    return actr.name
end

get_time(actr::AbstractACTR) = get_time(actr.scheduler)
get_time(scheduler) = scheduler.time

add_time!(actr::AbstractACTR, t) = add_time!(actr.scheduler, t)
add_time!(scheduler, t) = scheduler.time += t

reset_time!(actr::AbstractACTR) = reset_time!(actr.scheduler)
reset_time!(scheduler) = scheduler.time = 0.0

set_time!(actr::AbstractACTR, t) = set_time!(actr.scheduler, t)
set_time!(scheduler, t) = scheduler.time = t

rnd_time(μ) = rand(Uniform(μ * (2/3), μ * (4/3)))


text:
module Docile

import AnsiColor, Markdown

import Base: triplequoted, writemime, push!
import Base.Meta: isexpr

export query, doctest, @query, @doc, @docstrings, @tex_mstr

# internal
macro docref(ref)
    node = ref.args[end].args
    :($(esc(node[end])) = ($(node[1].args[1] + 1), @__FILE__))
end

include("types.jl")
include("docstrings.jl")

# Initialise docstrings for this module. Can't document previous files
# since they define the docstring methods/macros.
@docstrings

include("render.jl")
include("doctest.jl")
include("query.jl")

# link the readme into the package docs
@doc { :file => "../README.md" } -> Docile

include("docs.jl")

end # module


text:
"""
    IMERGMonthly{ST<:AbstractString, DT<:TimeType} <: IMERGDataset

Object containing information on Montly IMERG datasets to be downloaded
"""
struct IMERGMonthly{ST<:AbstractString, DT<:TimeType} <: IMERGDataset
	npdID :: ST
	lname :: ST
	doi   :: ST
    dtbeg :: DT
    dtend :: DT
    sroot :: ST
    hroot :: ST
    fpref :: ST
    fsuff :: ST
end

"""
    IMERGMonthly(
        ST = String,
        DT = Date;
        dtbeg :: TimeType,
        dtend :: TimeType,
        sroot :: AbstractString = homedir(),
    ) -> npd :: IMERGMonthly{ST,DT}

Creates a `IMERGMonthly` dataset `npd` to retrieve datasets for Monthly output

Keyword Arguments
=================
- `dtbeg` : Date at which download / analysis of the dataset begins
- `dtend` : Date at which download / analysis of the dataset ends
- `sroot` : The directory in which the folder `imergmonthly` will be created for data downloads, storage and analysis, default is the home directoy called by `homedir()`

The following fields in `npd` will be fixed as below:
- `npdID` : imergmonthly
- `lname` : IMERG Monthly
- `doi`   : 10.5067/GPM/IMERG/3B-MONTH/06
- `hroot` : https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGM.06
- `fpref` : 3B-MO.MS.MRG.3IMERG
- `fsuff` : V06B.HDF5
"""
function IMERGMonthly(
    ST = String,
    DT = Date;
    dtbeg :: TimeType,
    dtend :: TimeType,
    sroot :: AbstractString = homedir(),
)

	@info "$(modulelog()) - Setting up data structure containing information on IMERG Monthly data to be downloaded"

    fol = joinpath(sroot,"imergmonthly"); if !isdir(fol); mkpath(fol) end

	dtbeg = Date(year(dtbeg),1,1)
	dtend = Date(year(dtend),12,31)

    return IMERGMonthly{ST,DT}(
		"imergmonthly", "IMERG Monthly", "10.5067/GPM/IMERG/3B-MONTH/06",
        dtbeg, dtend,
		joinpath(sroot,"imergmonthly"),
        "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGM.06",
        "3B-MO.MS.MRG.3IMERG", "V06B.HDF5",
    )

end

function show(io::IO, npd::IMERGMonthly{ST,DT}) where {ST<:AbstractString, DT<:TimeType}
    print(
		io,
		"The NASA Precipitation Dataset {$ST,$DT} has the following properties:\n",
		"    Dataset ID      (npdID) : ", npd.npdID, '\n',
		"    Logging Name    (lname) : ", npd.lname, '\n',
		"    DOI URL          (doi)  : ", npd.doi,   '\n',
		"    Data Directory  (sroot) : ", npd.sroot, '\n',
		"    Date Begin      (dtbeg) : ", npd.dtbeg, '\n',
		"    Date End        (dtend) : ", npd.dtend, '\n',
		"    Timestep                : 1 Month\n",
        "    Data Resolution         : 0.1º\n",
        "    Data Server     (hroot) : ", npd.hroot, '\n',
	)
end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
export readpar, SpectralLines

#-------------------------------------------------------------------------------
#mappin between isotopologue labels and integer indices

const ISOINDEX = Dict{Char,Int64}(
    '1'=>1,  '2'=>2,  '3'=>3,  '4'=>4,  '5'=>5,  '6'=>6,
    '7'=>7,  '8'=>8,  '9'=>9,  '0'=>10, 'A'=>11, 'B'=>12,
    'C'=>13, 'D'=>14, 'E'=>15, 'F'=>16, 'G'=>17, 'H'=>18,
    'I'=>19, 'J'=>20, 'K'=>21, 'L'=>22, 'M'=>23, 'N'=>24,
    'O'=>25, 'P'=>26, 'Q'=>27, 'R'=>28, 'S'=>29, 'T'=>30,
    'U'=>31, 'V'=>32, 'W'=>33, 'X'=>34, 'Y'=>35, 'Z'=>36
)

#-------------------------------------------------------------------------------
# struct for molecule parameters

struct MolParam
    #molecule number
    M::Int64
    #molecular formula
    formula::String
    #molecule name
    name::String
    #global isotopologue codes
    I::Vector{Int64}
    #isopologue formulae
    isoform::Vector{String}
    #AFGL isotopologue codes
    AFGL::Vector{Int64}
    #abundance fractions
    A::Vector{Float64}
    #molecular masses [kg/mole]
    μ::Vector{Float64}
    #Qref
    Qref::Vector{Float64}
    #flag, has interpolating chebyshev polynomial
    hascheb::Vector{Bool}
    #length of cheby polys
    ncheb::Vector{Int64}
    #maximum rel err of cheb polys
    maxrelerr::Vector{Float64}
    #chebyshev expansion coefficients
    cheb::Vector{Vector{Float64}}
end

#for empty structs
MolParam() = MolParam(-1, "", "", [], [], [], [], [], [], [], [], [], [])

#-------------------------------------------------------------------------------

"""
    readpar(filename; νmin=0, νmax=Inf, Scut=0, I=[], maxlines=-1, progress=true)

Read an absoption line file from the HITRAN database, which should have the ".par" extension. These files are available at [`https://hitran.org/lbl`](https://hitran.org/lbl) after registering for a free account.

# Keyword Arguments
* `νmin`: smallest line wavenumber to include
* `νmax`: largest line wavenumber to include
* `Scut`: smallest spectral line intensity
* `I`: array of isotopologue numbers to include (excludes all others)
* `maxlines`: maximum number of lines to include (includes only the most intense `maxlines` lines)
* `progress`: whether to display the progress meter

A dictionary of vectors is returned, reflecting the definitions from
1. [HITRAN website](https://hitran.org/docs/definitions-and-units`)
2. [Rothman, Laurence S., et al. "The HITRAN 2004 molecular spectroscopic database." Journal of quantitative spectroscopy and radiative transfer 96.2 (2005): 139-204.](https://www.sciencedirect.com/science/article/abs/pii/S0022407313002859)

| Key | Vector Type | Description |
| --- | :---------- | :---------- |
| `M`    | `Int16`   | [HITRAN molecular identification number](https://hitran.org/docs/molec-meta) |
| `I`    | `Char`    | HITRAN isotopologue identification symbol |
| `ν`    | `Float64` | spectral line wavenumber [cm``^{-1}``] in a vacuum |
| `S`    | `Float64` | spectral line intensity [cm``^{-1}``/(molecule``\\cdot``cm``^{-2}``)] at 296 K |
| `A`    | `Float64` | Einstein-A coefficient (s``^{-1}``) of a transition |
| `γa`   | `Float64` | air-broadened half width at half maximum (HWHM) [cm``^{-1}``/atm] at 296 K and 1 atm |
| `γs`   | `Float64` | self-broadened half width at half maximum (HWHM) [cm``^{-1}``/atm] at 296 K and 1 atm |
| `Epp`  | `Float64` | lower-state energy of the transition [cm``^{-1}``] |
| `na`   | `Float64` | coefficient of temperature dependence of air-broadened half width |
| `δa`   | `Float64` | pressure shift [cm``^{-1}``/atm] at 296 K and 1 atm of the line position with respect to vacuum transition wavenumber |
| `Vp`   | `String`  | upper-state "global" quanta |
| `Vpp`  | `String`  | lower-state "global" quanta |
| `Qp`   | `String`  | upper-state "local" quanta |
| `Qpp`  | `String`  | lower-state "local" quanta |
| `Ierr` | `String`  | uncertainty indices |
| `Iref` | `String`  | reference indices |
| `*`    | `Char`    | flag (?) |
| `gp`   | `String`  | statistical weight of upper state |
| `gpp`  | `String`  | statistical weight of lower state |
"""
function readpar(filename::String;
                 νmin::Real=0,
                 νmax::Real=Inf,
                 Scut::Real=0,
                 I::Vector=[],
                 maxlines::Int=-1,
                 progress::Bool=true)
    @assert filename[end-3:end] == ".par" "expected file with .par extension, downloaded from https://hitran.org/lbl/"
    lines = readlines(filename)
    N = length(lines)
    #hard code the format to make sure nothing slows it down ¯\_(ツ)_/¯
    par = Dict(
        "M"   =>zeros(Int16, N),
        "I"   =>Vector{Char}(undef, N),
        "ν"   =>zeros(Float64, N),
        "S"   =>zeros(Float64, N),
        "A"   =>zeros(Float64, N),
        "γa"  =>zeros(Float64, N),
        "γs"  =>zeros(Float64, N),
        "Epp" =>zeros(Float64, N),
        "na"  =>zeros(Float64, N),
        "δa"  =>zeros(Float64, N),
        "Vp"  =>Vector{String}(undef, N),
        "Vpp" =>Vector{String}(undef, N),
        "Qp"  =>Vector{String}(undef, N),
        "Qpp" =>Vector{String}(undef, N),
        "Ierr"=>Vector{String}(undef, N),
        "Iref"=>Vector{String}(undef, N),
        "*"   =>Vector{Char}(undef, N),
        "gp"  =>Vector{String}(undef, N),
        "gpp" =>Vector{String}(undef, N)
    )
    #parse the values
    if progress
        prg = Progress(length(lines), 0.1, "Parsing \"$(splitdir(filename)[end])\" ")
    end
    @inbounds for i ∈ eachindex(lines)
        #get the line in question
        line = lines[i]
        #parse all the fields
        par["M"][i]    = parse(Int16, line[1:2])
        par["I"][i]    = line[3]
        par["ν"][i]    = parse(Float64, line[4:15])
        par["S"][i]    = parse(Float64, line[16:25])
        par["A"][i]    = parse(Float64, line[26:35])
        par["γa"][i]   = parse(Float64, line[36:40])
        par["γs"][i]   = parse(Float64, line[41:45])
        par["Epp"][i]  = parse(Float64, line[46:55])
        par["na"][i]   = parse(Float64, line[56:59])
        par["δa"][i]   = parse(Float64, line[60:67])
        par["Vp"][i]   = line[68:82]
        par["Vpp"][i]  = line[83:97]
        par["Qp"][i]   = line[98:112]
        par["Qpp"][i]  = line[113:127]
        par["Ierr"][i] = line[128:133]
        par["Iref"][i] = line[134:145]
        par["*"][i]    = line[146]
        par["gp"][i]   = line[147:153]
        par["gpp"][i]  = line[154:160]
        #update progress meter
        progress && next!(prg)
    end
    #filtering
    mask = ones(Bool, N)
    mask .&= par["ν"] .>= νmin
    mask .&= par["ν"] .<= νmax
    mask .&= par["S"] .>= Scut
    if length(I) > 0
        for j = 1:N
            #isotopologue character
            c = par["I"][j]
            #isotopologue integer
            i = ISOINDEX[c]
            #check if I and its integer counterpart are not in I
            if !(c in I) & !(i in I)
                mask[j] = false
            end
        end
    end
    #check that there will be information left
    @assert any(mask) "par information has been filtered to nothing!"
    #slice all the par arrays
    for (k,v) ∈ par
        par[k] = v[mask]
    end
    #take strongest lines if needed
    if maxlines > 0
        @assert maxlines > 0 "maxlines must be a positive integer"
        #check if there are fewer lines than desired already
        if N > maxlines
            idx = reverse(sortperm(par["S"]))[1:maxlines]
            for (k,v) in par
                par[k] = v[idx]
            end
        end
    end
    #ensure sorted by wavenumber
    idx = sortperm(par["ν"])
    for (k,v) ∈ par
        par[k] = v[idx]
    end
    return par
end

"""
Organizing type for spectral line data of a single gas

| Field | Type | Description |
| ----- | :--- | :---------- |
| `name`    | `String` | gas name |
| `formula` | `String` | gas formula |
| `N`       | `Int64` | number of lines |
| `M`       | `Int16` | see [`readpar`](@ref) |
| `I`       | `Vector{Int16}` | see [`readpar`](@ref) |
| `μ`       | `Vector{Float64}` | molar mass of isotopologues [kg/mole] |
| `A`       | `Vector{Float64}` | isotopologue abundance (Earth) |
| `ν`       | `Vector{Float64}` | see [`readpar`](@ref) |
| `S`       | `Vector{Float64}` | see [`readpar`](@ref) |
| `γa`      | `Vector{Float64}` | see [`readpar`](@ref) |
| `γs`      | `Vector{Float64}` | see [`readpar`](@ref) |
| `Epp`     | `Vector{Float64}` | see [`readpar`](@ref) |
| `na`      | `Vector{Float64}` | see [`readpar`](@ref) |

# Constructors

    SpectralLines(par::Dict)

Construct a `SpectralLines` object from a dictionary of line data. That dictionary can be created with [`readpar`](@ref).

    SpectralLines(filename, νmin=0, νmax=Inf, Scut=0, I=[], maxlines=-1)

Read a `.par` file directly into a `SpectralLines` object. Keyword arguments are passed through to [`readpar`](@ref).
"""
struct SpectralLines
    #gas name
    name::String
    #gas formula
    formula::String
    #number of lines
    N::Int64
    #molecule number (https://hitran.org/docs/molec-meta/)
    M::Int16
    #ordered isotopologue integer (https://hitran.org/docs/iso-meta/)
    I::Vector{Int16}
    #molar mass of isotopologues [kg/mole]
    μ::Vector{Float64}
    #isotopologue abundance in HITRAN
    A::Vector{Float64}
    #wavenumbers of absorption lines [cm^-1]
    ν::Vector{Float64}
    #spectral line intensity [cm^-1/(molecule*cm^-2)]
    S::Vector{Float64}
    #air-broadened half-width [cm^-1/atm]
    γa::Vector{Float64}
    #self-broadened half-width [cm^-1/atm]
    γs::Vector{Float64}
    #lower state energy, E'' [cm^-1]
    Epp::Vector{Float64}
    #temperature-dependence exponent for γair [unitless]
    na::Vector{Float64}
end

function SpectralLines(par::Dict)
    #length
    N = length(par["ν"])
    #make sure there is only one molecule present
    @assert length(unique(par["M"])) == 1 "SpectralLines objects must contain only one molecule's lines"
    #get the molecule name and formula
    M = par["M"][1]
    name = MOLPARAM[M].name
    form = MOLPARAM[M].formula
    #create arrays for isotopologue index, abundance, and molar mass
    I = map(i->ISOINDEX[par["I"][i]], 1:N)
    A = map(i->MOLPARAM[M].A[I[i]], 1:N)
    μ = map(i->MOLPARAM[M].μ[I[i]], 1:N)
    #get sorting indices to make sure ν vector is in ascending order
    idx = sortperm(par["ν"])
    #create a SpectralLines structure
    SpectralLines(
        name,
        form,
        N,
        M,
        I[idx],
        μ[idx],
        A[idx],
        par["ν"][idx],
        par["S"][idx],
        par["γa"][idx],
        par["γs"][idx],
        par["Epp"][idx],
        par["na"][idx]
    )
end

SpectralLines(fn::String; kwargs...) = SpectralLines(readpar(fn; kwargs...))


text:
printstyled(stdout,"  extended methods\n", color=:light_green)

printstyled(stdout,"    getindex\n", color=:light_green)
(S,T) = mktestseis()
@test findid(T, S) == [0, 0, 4, 5]

printstyled(stdout,"    getindex + Int on SeisData ==> SeisChannel\n", color=:light_green)
C = S[1]
@test S[firstindex(S)] == C
i_targ = 3
C = S[i_targ]
test_fields_preserved(C,S,i_targ)

printstyled(stdout,"    getindex + Range on SeisData ==> SeisData\n", color=:light_green)
D = S[i_targ:i_targ+1]
test_fields_preserved(D,S,2,i_targ+1)

printstyled(stdout,"    setindex!\n", color=:light_green)
A = SeisData(3)
setindex!(A, C, 3)
test_fields_preserved(C, A, 3)
A[1:2]=D
test_fields_preserved(A, S, 2, i_targ+1)
test_fields_preserved(C, S, 3)

printstyled(stdout,"    in\n", color=:light_green)
@test ("XX.TMP01.00.BHZ" in S.id)

printstyled(stdout,"      findid\n", color=:light_green)
@test ≈(findid("CC.LON..BHZ",S), findid(S,"CC.LON..BHZ"))
@test ≈(findid(S,"CC.LON..BHZ"), 4)
@test ≈(findid(S,C), 3)

printstyled(stdout,"    isempty\n", color=:light_green)
D = SeisData()
@test isempty(D)

printstyled(stdout,"    append!\n", color=:light_green)
(S,T) = mktestseis()
append!(S, T)
sizetest(S, 9)

C = deepcopy(S[4])
deleteat!(S, 4)
sizetest(S, 8)
@test ≈(length(findall(S.name.==C.name)),1)

C = deepcopy(S[3:4])
deleteat!(S,3:4)
nt = 6
@test ≈(S.n, nt)
@test ≈(maximum([length(getfield(S,i)) for i in datafields]), nt)
@test ≈(minimum([length(getfield(C,i)) for i in datafields]), 2)
@test length(findall(S.name.==C.name[1])).*length(findall(S.id.==C.id[1])) == 0
@test length(findall(S.name.==C.name[2])).*length(findall(S.id.==C.id[2])) == 1

s = "CC.LON..BHZ"
delete!(S, s)
sizetest(S, 5)

s = r"EH"
# @test_throws BoundsError S - s
delete!(S, s, exact=false)
sizetest(S, 2)

# untested methods in SeisData
for i = 1:5
  S = randSeisData()
  @test sizeof(S) > 0

  r = S.id[1]
  U = S - r
  @test sizeof(U) < sizeof(S)
  @test S.n == U.n + 1

  T = pull(S,r)
  @test isa(T, SeisChannel)
  @test(U == S)

  S = randSeisData(5)
  C = S[1]
  @test firstindex(S) == 1
  @test S[firstindex(S)] == S[1] == C
  @test length(S) == S.n == lastindex(S)

  push!(S, SeisChannel())
  T = prune(S)
  prune!(S)
  @test S == T
  @test S.n == T.n == 5
end

H = randSeisHdr()
@test sizeof(H) > 0
clear_notes!(H)
@test length(H.notes) == 1

(S,T) = mktestseis()
U = S-T
sizetest(S,5)
sizetest(T,4)
sizetest(U,3)

(S,T) = mktestseis()
U = deepcopy(S)
deleteat!(U, 1:3)
@test (S - [1,2,3]) == U
sizetest(S,5)

(S,T) = mktestseis()
@test in("UW.SEP..EHZ",S)
U = S[3:S.n]
V = deepcopy(S)
deleteat!(S, 1:2)
deleteat!(V, [1,2])
@test S == U == V
delete!(S, "CC.", exact=false)
delete!(U,V)

(S,T) = mktestseis()
X = deepcopy(T)
U = pull(S, 5)
@test U.id == "UW.SEP..EHZ"
sizetest(S, 4)
@test findid("UW.SEP..EHZ", S.id) == 0

U = pull(S, 3:4)
sizetest(S, 2)
@test findid(U.id[1], S.id) == 0
@test findid(U.id[2], S.id) == 0

V = pull(T, [2,3,4])
sizetest(T, 1)
@test findid(V.id[1], T.id) == 0
@test findid(V.id[2], T.id) == 0
@test findid(V.id[3], T.id) == 0
V = sort(V)

deleteat!(X,1)
@test findid(V,X) == [2,3,1]
Y = sort(X)
@test V == Y

# added 2019-02-23
S = SeisData(randSeisData(5), SeisChannel(), SeisChannel(),
    SeisChannel(id="UW.SEP..EHZ", name="Darth Exploded",
    loc=GeoLoc(lat=46.1967, lon=-122.1875, el=1440.0), t=[1 0; 1024 0], x=rand(1024)))
prune!(S)
@test (S.n == 6)
J = findchan("EHZ",S)
@test (6 in J)

printstyled("  SeisChannel methods\n", color=:light_green)

id = "UW.SEP..EHZ"
name = "Darth Exploded"

Ch = randSeisChannel()
Ch.id = id
Ch.name = name
S = SeisData(Ch)

@test in(id, Ch) == true
@test isempty(Ch) == false
@test convert(SeisData, Ch) == SeisData(Ch)
@test findid(Ch, S) == 1
@test sizeof(Ch) > 0
@test lastindex(S) == 1

Ch.gain = 1.0
Ch.fs = 0.0
@test isempty(Ch) == false