text:
# Several functions in this file are adapted from Frank T. Bergman
# Date: 2019
# Availability: https://groups.google.com/forum/#!topic/sbml-discuss/inS4Lzp3Ri8 or
# https://www.dropbox.com/s/2bfpiausejp0gd0/convert_reactions.py?dl=0
###################################################################################

@ModelingToolkit.parameters t
D = ModelingToolkit.Differential(t)
export D

"""
    simulatesbml(sbmlfile::String,tspan;saveat::Real=1.,jac::Bool=true,
                 solver=OrdinaryDiffEq.Tsit5())

Simulate model from SBML file and return an ODESolution.

# Arguments
- `sbmlfile::String`: path to SBML file.
- `tspan`: The solution u(t) will be computed for tspan[1] ≤ t ≤ tspan[2].
- `saveat::Real=1.`: time intervals of u(t) output.
- `jac::Bool=true`: Jacobian matrix.
- `solver=OrdinaryDiffEq.Tsit5()`: Integration algorithm.

# Example
```julia-repl
julia> sol = simulatesbml("mymodel.sbml",(0.0,1.0))
retcode: Success
Interpolation: 1st order linear
t: 2-element Array{Float64,1}:
 0.0
 1.0
u: 2-element Array{Array{Float64,1},1}:
 [0.0, 1.0]
 [0.4, 0.6]
```
"""
function simulatesbml(sbmlfile::String,tspan;saveat::Real=1.,jac::Bool=true,
                      solver=OrdinaryDiffEq.Tsit5())
    prob = sbml2odeproblem(sbmlfile,tspan,jac=jac)
    OrdinaryDiffEq.solve(prob,solver,saveat=saveat)
end

"""
    sbml2odeproblem(sbmlfile::String,tspan;jac::Bool=true)

Convert model from SBML file to an `ODEProblem` and return it.

# Arguments
- `sbmlfile::String`: path to SBML file.
- `tspan`: The solution u(t) will be computed for tspan[1] ≤ t ≤ tspan[2].
- `jac::Bool=true`: Jacobian matrix.

# Example
```julia-repl
julia> prob = sbml2odeproblem(SBML_FILE,(0.0,10.0))
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: [0.0, 1.0]
```
"""
function sbml2odeproblem(sbmlfile::String,tspan;jac::Bool=true)
    model = SbmlInterface.getmodel(sbmlfile)
    p = getparameters(model)
    u0 = getinitialconditions(model)
    sys,u0,p = sbml2odesystem(sbmlfile)
    ModelingToolkit.ODEProblem(sys,u0,tspan,p)
end

"""
    sbml2odesystem(sbmlfile::String)

Given an `sbmlfile` return a tuple of `ODESystem`, `u0` and `p`.

# Example
```julia-repl
julia> sys,u0,p = sbml2odesystem(SBML_FILE)
(Model ##ODESystem#268 with 2 equations
States (2):
  B(t)
  A(t)
Parameters (3):
  k1
  k2
  compartment,
Pair{Num,Float64}[A(t) => 1.0, B(t) => 0.0],
Pair{Num,Float64}[k1 => 0.8, k2 => 0.6, compartment => 1.0])
```
"""
function sbml2odesystem(sbmlfile::String)
    if !isfile(sbmlfile)
        throw(DomainError("`sbmlfile` is not a file"))
    elseif !(splitext(sbmlfile)[2] in (".xml", ".sbml"))
        throw(DomainError("`sbmlfile` must have `.xml` or `.sbml` ending"))
    end
    model = getmodel(sbmlfile)
    p = getparameters(model)
    u0 = getinitialconditions(model)
    eqs = getodes(model)
    sys = ModelingToolkit.ODESystem(eqs)
    sys,u0,p
end

"""
    getmodel(sbmlfile::String)

Given an `sbmlfile` return a PyObject (libsbml.Model).

# Example
```julia-repl
julia> model = getmodel(SBML_FILE)
PyObject <Model conversion_reaction_0 "Conversion Reaction 0">
```
"""
function getmodel(sbmlfile::String)
    doc = libsbml.readSBMLFromFile(sbmlfile)
    if doc.getNumErrors(libsbml.LIBSBML_SEV_FATAL) == 1
            @warn("Encountered serious errors while reading file")
            @error(doc.getErrorLog().toString())
    end

    # clear errors
    doc.getErrorLog().clearLog()

    # perform conversions
    props = libsbml.ConversionProperties()
    props.addOption("promoteLocalParameters", true)

    if doc.convert(props) != libsbml.LIBSBML_OPERATION_SUCCESS
        @warn("The document could not be converted")
        @warn(doc.getErrorLog().toString())
    end

    props = libsbml.ConversionProperties()
    props.addOption("expandFunctionDefinitions", true)  # Todo: ask PEtab developers set this to `True` when creating `petab.problem.Problem()`

    if doc.convert(props) != libsbml.LIBSBML_OPERATION_SUCCESS
        @warn("The document could not be converted")
        @warn(doc.getErrorLog().toString())
    end

    mod = doc.getModel()
end

"""
    getodes(model)::Array

Create array of ODEs from a libsbml Model.

# Example
```julia-repl
julia> model = getodes(SBML_FILE)
2-element Array{Equation,1}:
 Differential(t)(B(t)) ~ compartment*k1*A(t) - (compartment*k2*(B(t)))
 Differential(t)(A(t)) ~ compartment*k2*B(t) - (compartment*k1*(A(t)))
```
"""
function getodes(model)::Array
    assignments = []
    for a in model.getListOfRules()
        if (a.getMath().getName() == "piecewise") # || (a.getMath().getName() == nothing)
            @error("Assignment rules are not implemented.")
        end
        push!(assignments, a.getId() => a.getMath().getName())
    end

    reactions = Dict()  # Dict of reaction => kinetic formula
    for i in 0:model.getNumReactions()-1
        reaction = model.getReaction(i)
        kinetics = reaction.getKineticLaw()
        kinetic_components = kinetics.getFormula()
        reactions[reaction.getId()] = kinetic_components
    end

    species = Dict()  # Dict of species and stoichiometry-reactionId tuple they are involved in
    for i in 0:model.getNumSpecies()-1
        specie = model.getSpecies(i)
        if (specie.getBoundaryCondition() == true) || (specie.getId() in keys(species))
            continue
        end
        species[specie.getId()] = []
    end

    for i in 0:model.getNumReactions()-1
        reaction = model.getReaction(i)
        kinetics = reaction.getKineticLaw()
        for j in 0:reaction.getNumReactants()-1
            ref = reaction.getReactant(j)
            specie = model.getSpecies(ref.getSpecies())
            products = [r.getSpecies() for r in reaction.getListOfProducts()]
            if specie.getBoundaryCondition() == true
                println("continuing...")
                continue
            end
            stoich = -ref.getStoichiometry()

            if specie.getName() in products
                stoich_p = reaction.getProduct(specie.getId()).getStoichiometry()
                stoich = stoich_p + stoich
            end
            if stoich < 0
                stoich = string(stoich)
            elseif stoich > 0
                stoich = "+"*string(stoich)
            else
                stoich = ""
            end
            if stoich != ""
                push!(species[specie.getId()], (stoich, reaction.getId()))
            end
        end

        for j in 0:reaction.getNumProducts()-1
            ref = reaction.getProduct(j)
            specie = model.getSpecies(ref.getSpecies())
            reactants = [r.getSpecies() for r in reaction.getListOfReactants()]
            if (specie.getBoundaryCondition() == true) || (specie.getName() in reactants)
                println("continuing")
                continue
            end
            push!(species[specie.getId()], ("+"*string(ref.getStoichiometry()), reaction.getId()))
        end
    end

    getparameters(model)
    getinitialconditions(model)

    # Write ODEs
    eqs = ModelingToolkit.Equation[]
    for specie in keys(species)  # For every species
        if species[specie] != nothing
            lhs = eval(Meta.parse("D($specie)"))
            rhs = "0"
            for (coef, reaction_name) in species[specie]  # For every reaction
                reactionformula = " $coef * ( $(reactions[reaction_name]) )"
                reactionformula = replace(reactionformula, "pow"=>"^")
                rhs = rhs*reactionformula
            end
        rhs = eval(Meta.parse(rhs))
        eqn = ModelingToolkit.Equation(lhs, rhs)  # Todo: kick this line out
        push!(eqs, lhs ~ rhs)
        end
    end
    eqs
end

"""
    getparameters(model)

Create array of parameters from a libsbml Model.

# Example
```julia-repl
julia> pars = getparameters(model)
3-element Array{Pair{Num,Float64},1}:
          k1 => 0.8
          k2 => 0.6
 compartment => 1.0
```
"""
function getparameters(model)
    parameters = Pair{ModelingToolkit.Num,Float64}[]
    for i in 0:model.getNumParameters()-1
        par = model.getParameter(i)
        parname = par.getId()
        parval = par.getValue()
        eval(Meta.parse("@ModelingToolkit.parameters $parname"))
        eval(Meta.parse("export $parname"))
        push!(parameters, eval(Meta.parse("$parname => $parval")))
    end
    for i in 0:model.getNumCompartments()-1
        comp = model.getCompartment(i)
        compname = comp.getId()
        compval = comp.getSize()
        eval(Meta.parse("@ModelingToolkit.parameters $compname"))
        eval(Meta.parse("export $compname"))
        push!(parameters, eval(Meta.parse("$compname => $compval")))
    end
    parameters
end

"""
    getinitialconditions(model)

Create array of initial conditions from a libsbml Model.

# Example
```julia-repl
julia> u0 = getinitialconditions(model)
2-element Array{Pair{Num,Float64},1}:
 A(t) => 1.0
 B(t) => 0.0
```
"""
function getinitialconditions(model)
    eval(Meta.parse("@ModelingToolkit.parameters t"))

    initialassignments = Dict()
    for var in model.getListOfInitialAssignments()
        varname = var.getId()
        varval = var.getMath().getName()
        if varval isa String
            varval = model.getParameter(varval).getValue()
        end
        initialassignments[varname] = varval
    end

    initialconditions = Pair{ModelingToolkit.Num,Float64}[]
    for i in 0:model.getNumSpecies()-1
        var = model.getSpecies(i)
        varname = var.getId()
        if var.isSetInitialConcentration()
            varval = var.getInitialConcentration()
        elseif var.isSetInitialAmount()
            varval = var.getInitialAmount()
        else
            varval = initialassignments[varname]
        end
        eval(Meta.parse("@ModelingToolkit.variables $varname(t)"))
        eval(Meta.parse("export $varname"))
        push!(initialconditions, eval(Meta.parse("$varname => $varval")))
    end
    initialconditions
end


text:
@testset "Parse PDB and PDBML" begin

    txt(code) = joinpath(pwd(), "data", string(uppercase(code), ".pdb"))
    xml(code) = joinpath(pwd(), "data", string(uppercase(code), ".xml"))

    @testset "2VQC: Missings" begin

        code = "2VQC"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        @test findfirst(x -> x.id.number == "4", pdb) == findfirst(x -> x.id.number == "4", pdbml)
        @test findfirst(x -> x.id.number == "73",pdb) == findfirst(x -> x.id.number == "73",pdbml)
    end

    @testset "Test downloadpdb" begin

        code = "2VQC"
        pdb = read(txt(code), PDBFile)
        filename = downloadpdb(code)
        try
            pdbml = read(filename, PDBML)

            @test findfirst(x -> x.id.number == "4", pdb) == findfirst(x -> x.id.number == "4", pdbml)
            @test findfirst(x -> x.id.number == "73",pdb) == findfirst(x -> x.id.number == "73",pdbml)
        finally
            rm(filename)
        end

        filename = downloadpdb(code,
            headers=Dict("User-Agent" => "Mozilla/5.0 (compatible; MSIE 7.01; Windows NT 5.0)"),
            format=PDBFile)
        try
            d_pdb = read(filename, PDBFile)

            @test findfirst(x -> x.id.number == "4", pdb) == findfirst(x -> x.id.number == "4", d_pdb)
            @test findfirst(x -> x.id.number == "73",pdb) == findfirst(x -> x.id.number == "73",d_pdb)
        finally
            rm(filename)
        end
    end

    @testset "1H4A: Chain A (auth) == Chain X (label)" begin

        file = string(xml("1H4A"), ".gz")
        auth = read(file, PDBML, label=false)
        label = read(file, PDBML)

        @test unique([ res.id.chain for res in auth ]) == ["X"]
        @test unique([ res.id.chain for res in label]) == ["A", "B"]
    end

    @testset "1SSX: Residues with insert codes, 15A 15B" begin

        code = "1SSX"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        for residue_list in [pdb, pdbml]

            @test findall(res -> res.id.number == "15A", residue_list) == [1]
            @test findall(res -> isresidue(res,All,All,All,"15A"), residue_list) == [1]

            @test findall(res -> res.id.number == "15B", residue_list) == [2]
            @test findall(res -> isresidue(res,All,All,All,"15B"), residue_list) == [2]
        end

        @testset "@residues" begin

            @test (@residues pdb model All chain All group All residue "141") ==
                   filter(res -> isresidue(res, All, All, All,"141"), pdb)

            # Testing the macro in let block:
            mo = "1"
            ch = "A"
            gr = All
            re = "141"
            @test (@residues pdb model mo chain ch group gr residue re) ==
                  (@residues pdb model All chain All group gr residue "141")
        end

        @testset "Occupancy != 1.0 and @atom" begin

            @test sum(map(a -> (a.atom == "HH22" ? a.occupancy : 0.0), filter(r -> r.id.number == "141", pdbml)[1].atoms)) == 1.0
            @test sum([a.occupancy for a in @atoms pdbml model "1" chain "A"  group All residue "141" atom "HH22" ] ) == 1.0

            # Testing the macro in let block:
            mo = "1"
            ch = "A"
            gr = All
            re = "141"
            at = "HH22"
            @test sum([a.occupancy for a in @atoms pdbml model mo chain ch group gr residue re atom at ] ) == 1.0
        end

        @testset "Best occupancy" begin

            atoms_141 = @atoms    pdbml model "1" chain "A"  group All residue "141" atom "HH22"
            resid_141 = @residues pdbml model "1" chain "A"  group All residue "141"

            @test bestoccupancy(atoms_141)[1].occupancy == 0.75
            @test bestoccupancy(reverse(atoms_141))[1].occupancy == 0.75
            @test bestoccupancy(PDBAtom[ atoms_141[2] ])[1].occupancy == 0.25

            @test length(resid_141[1]) == 48
            @test selectbestoccupancy(resid_141[1], collect(1:48)) == 1
            @test selectbestoccupancy(resid_141[1], [1, 2]) == 1

            @test_throws AssertionError selectbestoccupancy(resid_141[1], Int[])
            @test_throws AssertionError selectbestoccupancy(resid_141[1], collect(1:100))
        end

        @testset "@atom with All" begin

            # ATOM      2  CA  ALA A  15A     22.554  11.619   6.400  1.00  6.14           C
            @test atoms(pdb, "1", "A", "ATOM", All, r"C.+")[1].atom == "CA"
        end
    end

    @testset "read with occupancyfilter=true" begin
        # `read` only atoms with the best occupancy

        code = "1SSX"
        pdb = read(txt(code), PDBFile, occupancyfilter=true)
        pdbml = read(xml(code), PDBML, occupancyfilter=true)

        res_pdb   = @residues pdbml model "1" chain "A"  group All residue "141"
        res_pdbml = @residues pdbml model "1" chain "A"  group All residue "141"

        atm_pdbml = @atoms pdbml model "1" chain "A"  group All residue "141" atom "HH22"

        @test length( atm_pdbml ) == 1
        @test atm_pdbml[1].occupancy == 0.75
        @test length( @atoms pdb   model "1" chain "A"  group All residue "141" atom "HH22" ) == 1

        @test length(res_pdb[1])   == 24
        @test length(res_pdbml[1]) == 24
    end

    @testset "CBN: Identical PDBe ResNum for Residue 22" begin
    # <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="22" dbResName="SER">
    #   <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="1cbn" dbResNum="22" dbResName="SER" dbChainId="A"/>
    #   <crossRefDb dbSource="UniProt" dbCoordSys="UniProt" dbAccessionId="P01542" dbResNum="22" dbResName="P"/>
    #   ....
    # <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="22" dbResName="PRO">
    #   <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="1cbn" dbResNum="22" dbResName="PRO" dbChainId="A"/>
    #   <crossRefDb dbSource="UniProt" dbCoordSys="UniProt" dbAccessionId="P01542" dbResNum="22" dbResName="P"/>
    #   ...

        code = "1CBN"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        @test length( @residues pdb   model "1" chain "A" group All residue "22" ) == 2
        @test length( @residues pdbml model "1" chain "A" group All residue "22" ) == 2

        @test [r.id.name for r in  @residues pdb   model "1" chain "A" group All residue "22"] == ["SER", "PRO"]
        @test [r.id.name for r in  @residues pdbml model "1" chain "A" group All residue "22"] == ["SER", "PRO"]
    end

    @testset "1AS5: NMR" begin

        code = "1AS5"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        @test length( @residues pdbml model "1"  chain "A" group All residue All ) == 25
        @test length( @residues pdbml model "14" chain "A" group All residue All ) == 25

        @test length( @residues pdbml model All  chain "A" group All residue All ) == 25*14
    end

    @testset "1DPO: Inserted residues lack insertion letters" begin
    # Single unnamed chain in 1DPO contains insertions at postions 184 (Gly, Phe),
    # 188 (Gly, Lys), and 221 (Ala, Leu) but no insertion letters.

        code = "1DPO"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        # Single "A" chain for PDB (auth_asym_id in PDBML)
        @test unique([r.id.chain for r in pdb  ]) == [ "A" ]
        # But 'A':'H' chains for PDBML (label_asym_id)
        @test unique([r.id.chain for r in pdbml]) == [ string(chain) for chain in 'A':'H' ]

        @test [r.id.name for r in @residues pdb   model "1" chain "A" group All residue r"^184[A-Z]?$"] == ["GLY", "PHE"]
        @test [r.id.name for r in @residues pdbml model "1" chain "A" group All residue r"^184[A-Z]?$"] == ["GLY", "PHE"]

        @test [r.id.name for r in @residues pdb   model "1" chain "A" group All residue r"^188[A-Z]?$"] == ["GLY", "LYS"]
        @test [r.id.name for r in @residues pdbml model "1" chain "A" group All residue r"^188[A-Z]*"] == ["GLY", "LYS"]

        @test [r.id.name for r in @residues pdb   model "1" chain "A" group All residue r"^221[A-Z]?$"] == ["ALA", "LEU"]
        @test [r.id.name for r in @residues pdbml model "1" chain "A" group All residue r"^221[A-Z]?$"] == ["ALA", "LEU"]
    end

    @testset "1IGY: Insertions" begin
    # Insertions have more than one copy of the same amino acid in a single insertion block.
    # For example, chain B in 1IGY contains a block of four residues inserted at sequence position 82.
    # The block contains Leu-Ser-Ser-Leu.

        code = "1IGY"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        @test [r.id.name for r in @residues pdb   model "1" chain "B" group All residue r"^82[A-Z]?$" ] ==
                ["LEU", "SER", "SER", "LEU"]
        @test [r.id.name for r in @residues pdbml model "1" chain "B" group All residue r"^82[A-Z]?$" ] ==
                ["LEU", "SER", "SER", "LEU"]

        @test sum([r.id.group for r in @residues pdb model "1" chain "D" group All residue All] .== "HETATM") ==
                length(@residues pdb model "1" chain "D" group "HETATM" residue All)
        @test sum([r.id.group for r in @residues pdb model "1" chain "D" group All residue All] .== "ATOM") ==
                length(@residues pdb model "1" chain "D" group "ATOM" residue All)
    end

    @testset "1HAG" begin
    # Chain E begins with 1H, 1G, 1F, ... 1A, then 1 (in reverse alphabetic order)

        code = "1HAG"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        @test unique([res.id.chain for res in pdb  ]) == ["E", "I"]
        @test unique([res.id.chain for res in pdbml]) == ["A", "B", "C", "D", "E"]

        # The chain E of PDB is the chain A of PDBML
        @test [r.id.number for r in @residues pdb   model "1" chain "E" group All residue r"^1[A-Z]?$"] ==
                [string(1, code) for code in vcat(collect('H':-1:'A'), "")]
        @test [r.id.number for r in @residues pdbml model "1" chain "A" group All residue r"^1[A-Z]?$"] ==
                [string(1, code) for code in vcat(collect('H':-1:'A'), "")]
    end

    @testset "1NSA" begin
    # Contains a single (unnamed) protein chain with sequence 7A-95A that continues 4-308.

        code = "1NSA"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        # Single "A" chain for PDB (auth_asym_id in PDBML)
        @test unique([r.id.chain for r in pdb  ]) == ["A"]
        # But 'A':'F' chains for PDBML (label_asym_id)
        @test unique([r.id.chain for r in pdbml]) == [string(chain) for chain in 'A':'F']

        ind = findall(r -> r.id.number=="95A", pdbml)[1]
        @test pdbml[ind + 1].id.number == "4"

        ind = findall(r -> r.id.number=="95A", pdb)[1]
        @test pdb[ind + 1].id.number == "4"
    end

    @testset "1IAO" begin
    # Contains in chain B (in this order) 1S, 2S, 323P-334P, 6-94, 94A, 95-188, 1T, 2T

        code = "1IAO"
        pdb = read(txt(code), PDBFile)
        pdbml = read(xml(code), PDBML)

        pdb_B   = @residues pdb   model "1" chain "B" group All residue All
        pdbml_B = @residues pdbml model "1" chain "B" group All residue All

        for B in [pdb_B, pdbml_B]
            @test B[findall(r -> r.id.number == "2S", B)[1]   + 1].id.number == "323P"
            @test B[findall(r -> r.id.number == "334P", B)[1] + 1].id.number == "6"
            @test B[findall(r -> r.id.number == "94", B)[1]   + 1].id.number == "94A"
            @test B[findall(r -> r.id.number == "94A", B)[1]  + 1].id.number == "95"
            @test B[findall(r -> r.id.number == "188", B)[1]  + 1].id.number == "1T"
            @test B[findall(r -> r.id.number == "1T", B)[1]   + 1].id.number == "2T"
        end
    end
end

@testset "RESTful PDB Interface" begin

    @test getpdbdescription("4HHB")["resolution"] == "1.74"
    @test_throws KeyError getpdbdescription("104D")["resolution"] # NMR
end

@testset "Write PDB files" begin

    txt(code) = joinpath(pwd(), "data", string(uppercase(code), ".pdb"))

    @testset "2VQC" begin
        code = "2VQC"
        io = IOBuffer()
        pdb = read(txt(code), PDBFile)
        print(io, pdb, PDBFile)
        printed = split(String(take!(io)), '\n')

        @test length(printed) == 609 # Only ATOM, HETATM & END + 1 because the trailing \n
        @test printed[1]   == "ATOM      1  N   THR A   4       2.431  19.617   6.520  1.00 24.37           N  "
        @test printed[607] == "HETATM  607  O   HOH A2025      13.807  38.993   2.453  1.00 33.00           O  "
    end

    @testset "NMR" begin
        code = "1AS5"
        io = IOBuffer()

        pdb = read(txt(code), PDBFile)
        print(io, pdb, PDBFile)
        printed = split(String(take!(io)), '\n')

        @test sum(map(x -> startswith(x, "MODEL "), printed)) == 14 # 14 models
        @test sum(map(x -> x == "ENDMDL", printed)) == 14 # 14 models
    end

    @testset "2 Chains" begin
        code = "1IAO"
        io = IOBuffer()

        pdb = read(txt(code), PDBFile)
        print(io, pdb, PDBFile)
        printed = split(String(take!(io)), '\n')

        # MIToS only prints TER for the ATOM group if the chain changes.
        # Some modified residues are annotated as HETATM in the middle of the ATOM chain:
        # TER can not be printed from ATOM to HETATM if the chain doesn’t change.

        # Only prints TER between chain A and B
        @test sum(map(x -> startswith(x, "TER "), printed)) == 1

        @test filter!(s -> occursin(r"TER ", s), printed)[1] == "TER    1418      TRP A 178 "
    end

    @testset "read/write consistency" begin
       io = IOBuffer()
       for code in ["2VQC", "1IAO", "1NSA", "1HAG", "1IGY", "1DPO", "1AS5", "1CBN", "1SSX"]

           readed = read(txt(code), PDBFile)
           print(io, readed, PDBFile)
           readed_writed_readed = parse(String(take!(io)), PDBFile)

           @test readed_writed_readed == readed
       end
    end
end

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("onlineuvstats.jl")
include("onlinemvstats.jl")
include("effective_sample_size.jl")


text:
adv_eq_dir = joinpath(output_root, "AdvectionEquation")

print_norms = false

plot_solution(grid, aux, wave_speed, filename) = nothing
plot_solution_burgers(grid, aux, velocity_sign, filename) = nothing

using Requires
@init @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
    using .Plots

    function plot_solution(grid, aux, wave_speed, filename)
        mkpath(dirname(filename))
        domain_range = over_elems_real(grid)
        x = [grid.zc[k] for k in domain_range]
        if wave_speed == -1
            y = [aux[:ϕ_initial, k] for k in domain_range]
            plot(
                y,
                x,
                label = "initial",
                markercolor = "blue",
                linecolor = "blue",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
            y = [q[:ϕ, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "numerical",
                markercolor = "black",
                linecolor = "black",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
            y = [aux[:ϕ_error, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "error",
                markercolor = "red",
                linecolor = "red",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
            y = [aux[:ϕ_analytic, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "analytic",
                markercolor = "green",
                linecolor = "green",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
        else
            y = [q[:ϕ, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "",
                markercolor = "black",
                linecolor = "black",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
            y = [aux[:ϕ_error, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "",
                markercolor = "red",
                linecolor = "red",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
            y = [aux[:ϕ_analytic, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "",
                markercolor = "green",
                linecolor = "green",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
        end
        wave_speed == 1 && savefig(filename)
    end

    function plot_solution_burgers(grid, aux, velocity_sign, filename)
        mkpath(dirname(filename))
        domain_range = over_elems_real(grid)
        x = [grid.zc[k] for k in domain_range]
        if velocity_sign == -1
            y = [aux[:w_initial, k] for k in domain_range]
            plot(
                y,
                x,
                label = "initial",
                markercolor = "blue",
                linecolor = "blue",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
            y = [q[:w, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "numerical",
                markercolor = "black",
                linecolor = "black",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
        else
            y = [aux[:w_initial, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "",
                markercolor = "blue",
                linecolor = "blue",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
            y = [q[:w, k] for k in domain_range]
            plot!(
                y,
                x,
                label = "",
                markercolor = "black",
                linecolor = "black",
                markershapes = markershape,
                markersize = 2,
                legend = :topleft,
            )
        end
        velocity_sign == 1 && savefig(joinpath(directory, name))
    end

end

@testset "Linear advection, ∂_t ϕ + ∇•(cϕ) = 0 ∈ ∂Ω, ϕ(t=0) = Gaussian(σ, μ), ConservativeForm, explicit Euler" begin
    domain_set = DomainSet(gm = 1)
    domain_subset = DomainSubSet(gm = true)
    σ, μ = 0.1, 0.5
    δz = 0.2
    Triangle(z) = μ - δz < z < μ + δz ?
        (μ > z ? (z - (μ - δz)) / δz : ((μ + δz) - z) / δz) : 0.0
    Square(z) = μ - δz < z < μ + δz ? 1 : 0.0
    Gaussian(z) = exp(-1 / 2 * ((z - μ) / σ)^2)
    for n_elems_real in (64, 128)
        grid = UniformGrid(0.0, 1.0, n_elems_real)
        domain_range = over_elems_real(grid)
        x = [grid.zc[k] for k in domain_range]
        unknowns = ((:ϕ, domain_subset),)
        vars = (
            (:ϕ_initial, domain_subset),
            (:ϕ_error, domain_subset),
            (:ϕ_analytic, domain_subset),
        )
        q = StateVec(unknowns, grid, domain_set)
        aux = StateVec(vars, grid, domain_set)
        rhs = deepcopy(q)
        CFL = 0.1
        Δt = CFL * grid.Δz
        T = 0.25
        maxiter = Int(T / Δt)
        for scheme in
            (UpwindAdvective(), UpwindCollocated(), CenteredUnstable())
            for distribution in (Triangle, Square, Gaussian)
                for wave_speed in (-1, 1)
                    scheme_name = replace(string(scheme), "()" => "")
                    distribution_name =
                        joinpath("LinearAdvection", string(distribution))
                    directory =
                        joinpath(adv_eq_dir, distribution_name, scheme_name)
                    print_norms && print(
                        "\n",
                        directory,
                        ", ",
                        n_elems_real,
                        ", ",
                        wave_speed,
                        ", ",
                    )

                    for k in over_elems_real(grid)
                        aux[:ϕ_initial, k] = distribution(grid.zc[k])
                        q[:ϕ, k] = aux[:ϕ_initial, k]
                    end
                    amax_w = max([
                        max(aux[:ϕ_initial, k]) for k in over_elems_real(grid)
                    ]...)
                    for i in 1:maxiter
                        for k in over_elems_real(grid)
                            ϕ = q[:ϕ, Cut(k)]
                            ϕ_dual = q[:ϕ, Dual(k)]
                            w = [wave_speed, wave_speed, wave_speed]
                            w_dual = [wave_speed, wave_speed]
                            rhs[:ϕ, k] =
                                -advect_old(
                                    ϕ,
                                    ϕ_dual,
                                    w,
                                    w_dual,
                                    grid,
                                    scheme,
                                    Δt,
                                )
                        end
                        for k in over_elems(grid)
                            q[:ϕ, k] += Δt * rhs[:ϕ, k]
                        end
                        apply_Dirichlet!(q, :ϕ, grid, 0.0, Zmax())
                        apply_Dirichlet!(q, :ϕ, grid, 0.0, Zmin())
                    end
                    sol_analtyic = [
                        distribution(grid.zc[k] - wave_speed * T)
                        for k in over_elems(grid)
                    ]
                    sol_error =
                        [sol_analtyic[k] - q[:ϕ, k] for k in over_elems(grid)]
                    L2_norm = sum(sol_error .^ 2) / length(sol_error)

                    if !(scheme == CenteredUnstable())
                        @test all(abs.(sol_error) .< 100 * grid.Δz)
                        @test all(L2_norm < 100 * grid.Δz)
                    end
                    for k in over_elems(grid)
                        aux[:ϕ_error, k] = sol_error[k]
                        aux[:ϕ_analytic, k] = sol_analtyic[k]
                    end

                    print_norms && print("L2_norm(err) = ", L2_norm)

                    name = string(n_elems_real)
                    markershape = wave_speed == -1 ? :dtriangle : :utriangle
                    # Skip plotting...
                    # plot_solution(
                    #     grid,
                    #     aux,
                    #     wave_speed,
                    #     joinpath(directory, name),
                    # )
                end
            end
        end
    end
end

@testset "Non-linear Bergers, ∂_t w + ∇•(ww) = 0 ∈ ∂Ω, u(t=0) = Gaussian(σ, μ), ConservativeForm, explicit Euler" begin
    domain_set = DomainSet(gm = 1)
    domain_subset = DomainSubSet(gm = true)
    σ, μ = 0.1, 0.5
    δz = 0.2
    Triangle(z, velocity_sign) = μ - δz < z < μ + δz ?
        (
        μ > z ? velocity_sign * (z - (μ - δz)) / δz :
            velocity_sign * ((μ + δz) - z) / δz
    ) :
        0.0
    Square(z, velocity_sign) = μ - δz < z < μ + δz ? velocity_sign : 0.0
    Gaussian(z, velocity_sign) = velocity_sign * exp(-1 / 2 * ((z - μ) / σ)^2)
    for n_elems_real in (64, 128)
        grid = UniformGrid(0.0, 1.0, n_elems_real)
        domain_range = over_elems_real(grid)
        x = [grid.zc[k] for k in domain_range]
        unknowns = ((:w, domain_subset),)
        vars = (
            (:w_initial, domain_subset),
            (:w_error, domain_subset),
            (:w_analytic, domain_subset),
        )
        q = StateVec(unknowns, grid, domain_set)
        aux = StateVec(vars, grid, domain_set)
        rhs = deepcopy(q)
        CFL = 0.1
        Δt = CFL * grid.Δz
        T = 0.25
        maxiter = Int(T / Δt)
        for scheme in
            (UpwindAdvective(), UpwindCollocated(), CenteredUnstable())
            for distribution in (Triangle, Square, Gaussian)
                for velocity_sign in (-1, 1)
                    for k in over_elems_real(grid)
                        aux[:w_initial, k] =
                            distribution(grid.zc[k], velocity_sign)
                        q[:w, k] = aux[:w_initial, k]
                    end
                    amax_w = max([
                        max(aux[:w_initial, k]) for k in over_elems_real(grid)
                    ]...)
                    for i in 1:maxiter
                        for k in over_elems_real(grid)
                            w = q[:w, Cut(k)]
                            w_dual = (w[1:2] + w[2:3]) / 2
                            rhs[:w, k] =
                                -advect_old(
                                    w,
                                    w_dual,
                                    w,
                                    w_dual,
                                    grid,
                                    scheme,
                                    Δt,
                                )
                        end
                        for k in over_elems(grid)
                            q[:w, k] += Δt * rhs[:w, k]
                        end
                        apply_Dirichlet!(q, :w, grid, 0.0, Zmax())
                        apply_Dirichlet!(q, :w, grid, 0.0, Zmin())
                    end

                    scheme_name = replace(string(scheme), "()" => "")
                    distribution_name =
                        joinpath("BurgersEquation", string(distribution))
                    directory =
                        joinpath(adv_eq_dir, distribution_name, scheme_name)
                    name = string(n_elems_real)
                    markershape = velocity_sign == -1 ? :dtriangle : :utriangle
                    plot_solution_burgers(
                        grid,
                        aux,
                        velocity_sign,
                        joinpath(directory, name),
                    )
                end
            end
        end
    end

end


text:
@testset "TDVP" begin
    @testset "DAS" begin
        l = 6
        pchs = U1Charges(-2:2)
        vchs = U1Charges(-4:4)
        pds = [1 for i in pchs]
        vds = [5 for i in vchs]
        mpo = hubbardMPO(1, 1, pchs,l, T = ComplexF64)
        mps = randmps(l,ComplexF64, (U1(), (pchs, vchs), (pds, vds)))
        gs, es  = dmrg(mps, mpo)

        cmps = normalize(gs)
        mpss = tdvp(cmps,mpo,10,0.01, verbose = false, collectall = true)
        vals = [sqrt(abs2(inner(cmps, mps))) for mps in mpss[2]]
        @test all(vals .≈ 1)

        cmps = normalize(mps)
        mpss = tdvp(cmps,mpo,10,0.01im, verbose = false, collectall = true)
        vals = [sqrt(abs2(inner(cmps, mps))) for mps in mpss[2]]
        @test all(diff(vals) .< 0)
    end

    @testset "Dense" begin
        l = 6
        mpo = hubbardMPO(1, 1, U1Charges(-2:2),l, T = ComplexF64)
        ansatz = randmps(l, ComplexF64, (5,20))
        gsa, eas = dmrg(ansatz, todense(mpo))

        cmpsa = normalize(ansatz)
        mpsas = tdvp(cmpsa,todense(mpo),10,0.01, verbose = false, collectall = true)
        vals = [sqrt(abs2(inner(cmpsa, mps))) for mps in mpsas[2]]
        @test all(diff(vals) .< 0)

        cmpsa = normalize(gsa)
        mpsas = tdvp(cmpsa,todense(mpo),10,0.01, verbose = false, collectall = true)
        vals = [sqrt(abs2(inner(cmpsa, mps))) for mps in mpsas[2]]
        @test all(vals .≈ 1)
    end

    @testset "ED" begin
        l = 3
        pchs = U1Charges(-1:1)
        vchs = U1Charges(-5:5)
        pds = [1 for i in pchs]
        vds = [5 for i in vchs]
        @testset "Dense" begin
            #= hamiltonian=#
            mpo = todense(hubbardMPO(1, 1, pchs,l, T = ComplexF64))
            left, blk, right = mpo[[1,2,l]]
            @tensor h[1,2,3,4,5,6] := left[1,-1,4] * blk[2,-1,-2,5] * right[3,-2,6];
            harr = toarray(fuselegs(h,((1,2,3),(4,5,6)))[1])

            #= state =#
            mps = todense(convert(MPS,normalize(randmps(l,ComplexF64,(U1(), (pchs,vchs),(pds,vds))))))
            v1,v2,v3 = mps[[1,2,3]]
            @tensor v[1,2,3] := v1[1,-1] * v2[-1,2,-2] * v3[-2,3]

            t = 1e-0
            #= ED =#
            varr = toarray(fuselegs(v,((1,2,3),))[1])
            varr2 = exp(-1im * t * harr) * varr
            norm(varr) ≈ norm(varr2) ≈ 1
            sqrt(abs2(dot(varr,varr2)))

            #= tdvp =#
            mps2 = convert(MPS,tdvp(mps, mpo, 10, t))
            @test inner(mps2) ≈ inner(mps) ≈ 1
            res = sqrt(abs2(inner(mps2,mps)))

            v1p,v2p,v3p = mps2[[1,2,3]]
            @tensor v2[1,2,3] := v1p[1,-1] * v2p[-1,2,-2] * v3p[-2,3]
            varr2p = toarray(fuselegs(v2,((1,2,3),))[1])

            @test varr2p ≈ varr2
        end

        @testset "DAS" begin
            l = 3
            mpo = hubbardMPO(1, 1, pchs,l, T = ComplexF64)
            left, blk, right = mpo[[1,2,l]]
            @tensor h[1,2,3,4,5,6] := left[1,-1,4] * blk[2,-1,-2,5] * right[3,-2,6];
            harr = toarray(fuselegs(h,((1,2,3),(4,5,6)),InOut(1,-1))[1])

            mps  = convert(MPS,normalize(randmps(l,ComplexF64, (U1(), (pchs,vchs),(pds,vds)))))
            v1,v2,v3 = mps[[1,2,3]]
            @tensor v[1,2,3] := v1[1,-1] * v2[-1,2,-2] * v3[-2,3]

            t = 1e-0
            #= ED =#
            varr = toarray(fuselegs(v,((1,2,3),))[1])
            varr2 = exp(-1im * t * harr) * varr
            norm(varr) ≈ norm(varr2) ≈ 1
            sqrt(abs2(dot(varr,varr2)))

            #= tdvp =#
            mps2 = convert(MPS,tdvp(mps, mpo, 10, t))
            @test inner(mps2) ≈ inner(mps) ≈ 1
            res = sqrt(abs2(inner(mps2,mps)))

            v1p,v2p,v3p = mps2[[1,2,3]]
            @tensor v2[1,2,3] := v1p[1,-1] * v2p[-1,2,-2] * v3p[-2,3]
            varr2p = toarray(fuselegs(v2,((1,2,3),),InOut(1))[1])
            @test varr2p ≈ varr2
        end
    end
end


text:
using Distributions, Random, Statistics, Test
import ThorinDistributions
import DynamicPolynomials, Statistics, StatsBase

####################################
# Previous versions of the code, used now for testing purposes only.
####################################
function Bell_partial_matrix(N,K,x)
    M = zeros(eltype(x .+ 1.0),(N+1,K+1))
    for n in 0:N
        for k in 0:n
            if n == 0 && k == 0
                M[n+1,k+1] = 1
            elseif n == 0 || k == 0
                M[n+1,k+1] = 0
            else
                M[n+1,k+1] = sum([binomial(n-1,i-1)*x[i]*M[n-i+1,k] for i in 1:n-k+1])
            end
        end
    end
    return M
end
function cumulants_from_moments(μ, mu0 = 1)
    # we ned to construct kappa that are the derivatives of the log times mu
    n = length(μ)
    M = Bell_partial_matrix(n-1,n-1,μ)
    log_der_fact = [log(mu0),[factorial(big(i-1)) * (-1)^(i-1) * mu0^(-i) for i in 1:(n-1)]...]
    return M''log_der_fact
end
# According to wikipedia, we could also obtain them by recursion.
# see end of paragraph https://en.wikipedia.org/wiki/Cumulant#Cumulants_and_moments
function cum_from_mom_rec!(κ,μ,mu0=nothing)
    if isnothing(mu0)
        mu0 = μ[1]
        μ = μ[2:end]
    end
    # much faster algorithm !
    κ[1] = log(mu0)
    κ[2] = μ[1] ./ mu0
    for n in 2:(length(μ)-1)
        κ[n+1] = (μ[n] - sum(binomial(big(n-1),big(m-1))*κ[m+1]*μ[n-m] for m in 1:(n-1))) / mu0
    end
    return κ
end
function cum_from_mom_rec(μ;mu0)
    if isnothing(mu0)
        mu0 = μ[1]
        μ = μ[2:end]
    end
    κ = Array{Any}(undef, length(μ)) # needs to be Array{Any} so that the test with polynomials passes...
    cum_from_mom_rec!(κ,μ,mu0)
    return κ
end
function old_thorin_moments(D,t_star,m)
    n = length(D)
    mu = sum(D .^(0:(m+1))' .* exp.(t_star .* D),dims=1)/n
    return [log(mu[1]),(cum_from_mom_rec(mu[2:end],mu0=mu[1])[2:end] ./ factorial.(big.(0:(m-1))))...]
end
function old_resemps_thorin_moments(M,D,t_star,m)
    n = length(D)
    T = eltype(D)
    mu_mat = D .^(0:(m+1))' .* exp.(t_star .* D)
    kappa = zeros(T,(M,m+1))
    facts = [T(1),factorial.(T.(0:(m-1)))...]
    for i in 1:M
        @views kappa[i,:] = cum_from_mom_rec!(kappa[i,:], StatsBase.mean(mu_mat[StatsBase.sample(1:n,n,replace=true),:],dims=1)) ./ facts
    end
    return kappa
end
function cum_from_mom_rec_simplified(μ;mu0)
    # for test only
    if isnothing(mu0)
        mu0 = μ[1]
        μ = μ[2:end]
    end
    η = Array{Any}(undef, length(μ)+1)
    τ = Array{Any}(undef, length(μ)+1)
    η .= [mu0,(μ ./ factorial.(big.(1:length(μ))))...]
    η[2:end] ./= mu0
    τ[1] = log.(mu0)
    τ[2] = η[2]
    for k in 3:length(η)
        τ[k] = (k-1)*η[k]
        for j in 1:(k-2)
            τ[k] -= τ[j+1] * η[k-j]
        end
    end
    κ = τ .* [1,factorial.(big.(0:(length(τ)-2)))...]
    return κ[1:end-1]
end

@testset "ThorinMoments.jl tests" begin

    DynamicPolynomials.@polyvar μ[1:10];
    κ1 = cumulants_from_moments(μ,1);
    κ2 = cum_from_mom_rec(μ,mu0=1);
    κ3 = cum_from_mom_rec_simplified(μ,mu0=1);
    κ4 = cumulants_from_moments(μ,0.5);
    κ5 = cum_from_mom_rec(μ,mu0=0.5);
    κ6 = cum_from_mom_rec_simplified(μ,mu0=0.5);

    n = 100
    D = BigFloat.(rand(LogNormal(),n));

    @testset "thorin_moment_function is alright" begin
        τ = ThorinDistributions.thorin_moments(D,-1,9)
        ηs = zeros(eltype(D),(n,11))
        ThorinDistributions.ηs_from_data!(ηs,D,-1)
        η_mom = Statistics.mean(ηs,dims=1) .* factorial.(big.(0:10))'
        τ1 = cum_from_mom_rec(η_mom,mu0=nothing) ./ factorial.([big(0),big.(0:8)...])
        τ2 = cum_from_mom_rec_simplified(η_mom,mu0=nothing) ./ factorial.([big(0),big.(0:8)...])
        τ3 = cumulants_from_moments(η_mom[2:end],η_mom[1]) ./ factorial.([big(0),big.(0:8)...])

        @test τ ≈ τ1
        @test τ ≈ τ2
        @test τ ≈ τ3
    end

    @testset "firsts cumulants corretness, mu0=1" begin
        @test κ1[1] == 0
        @test κ1[2] == μ[1]
        @test κ1[3] == μ[2] - μ[1]^2
        @test κ1[4] == 2.0μ[1]^3 - 3.0μ[1]μ[2] + μ[3]
        @test κ1[5] ==  -6.0μ[1]^4 + 12.0μ[1]^2*μ[2] - 4.0μ[1]μ[3] - 3.0μ[2]^2 + μ[4]
        @test κ1[6] ==   24.0μ[1]^5 - 60.0μ[1]^3*μ[2] + 20.0μ[1]^2*μ[3] + 30.0μ[1]*μ[2]^2 - 5.0μ[1]*μ[4] - 10.0μ[2]*μ[3] + μ[5]
    end

    @testset "cumulants recursivity" begin
        @test κ1 == κ2
        @test κ4 == κ5
    end

    @testset "η/τ simplification" begin
        @test all(κ3 .≈ κ1)
        @test all(κ3 .≈ κ2)
        @test all(κ6 .≈ κ4)
        @test all(κ6 .≈ κ5)
    end

    @testset "Bell partial matrix" begin
        M1 = [1   0      0       0       0       0       0      0     0   0  0
        0   1      0       0       0       0       0      0     0   0  0
        0   2      1       0       0       0       0      0     0   0  0
        0   3      6       1       0       0       0      0     0   0  0
        0   4     24      12       1       0       0      0     0   0  0
        0   5     80      90      20       1       0      0     0   0  0
        0   6    240     540     240      30       1      0     0   0  0
        0   7    672    2835    2240     525      42      1     0   0  0
        0   8   1792   13608   17920    7000    1008     56     1   0  0
        0   9   4608   61236  129024   78750   18144   1764    72   1  0
        0  10  11520  262440  860160  787500  272160  41160  2880  90  1]
        @test Int.(Bell_partial_matrix(10,10,1:30)) == M1
    end

    @testset "Resemp is OK" begin
        M = 10
        Random.seed!(123)
        x1 = old_resemps_thorin_moments(M,D,-1,10)
        Random.seed!(123)
        x2 = ThorinDistributions.resemps_thorin_moments(M,D,-1,10)
        @test x1 ≈ x2
        Random.seed!(123)
        x3 = ThorinDistributions.resemps_thorin_moments(M,[D 2D]',-1,10)
        @test x1 ≈ x3[:,1,:]

    end
end


text:
using Flux, Plots
using Base.Iterators: repeated
using Flux: @epochs
using LinearAlgebra

gridsize = 150;
tstart=0.0
tend=20.0
myfunc(x) = sin(x);
x = collect(range(tstart,stop=tend,length=gridsize));
y= myfunc.(x);

xtest=collect(tend:0.1:2*tend)
xpred=collect(tstart:0.1:2*tend)

#Build the input data
data = []
for i in 1:length(x)
    push!(data, ([x[i]], y[i])) #have to input a list
end

Q = 20;
ann = Chain(Dense(1,Q,tanh),Dense(Q,Q,tanh),Dense(Q,1));

#Simple Mean-square loss2
function loss(x, y)
    pred=ann(x)
    loss=Flux.mse(ann(x), y)
    return loss
end

opt = ADAM()

ps=params(ann)
@epochs 5000 Flux.train!(loss,ps, data, opt)

plot(xpred, ann(xpred')', color=:red, lw=2.0, label="")
scatter!(x,myfunc.(x), color=:blue, legend=false)
scatter!(xtest,myfunc.(xtest), color=:green, lw=2.0, label="", markershape=:cross)
ylims!((-1.5,1.5))
xlims!((tstart,2*tend))
title!("Curve Fitting with Neural Networks")
xlabel!("Input")
ylabel!("Output")

savefig("sin_poor.png")

max_stars_repo_name: khannay/CurveFittingNeuralNets

max_issues_repo_name: khannay/CurveFittingNeuralNets

max_forks_repo_name: khannay/CurveFittingNeuralNets


text:
# fin_example.jl
# Script to provide example of how the simulation method runs.

using Statistics

include("fin_simulate.jl")

Random.seed!(111)

# Generate a random set of patients, with varying admission times over a horizon of 30 weeks ( = 210 days).
num_patients = 5000
num_days = 7 * 30;

day_of_entry = rand(1:(num_days-15), num_patients)
rday_of_surgery_completion = rand(1:3, num_patients)
rday_of_discharge = rday_of_surgery_completion + rand(1:5, num_patients)

# Generate the endpoint (0/1 binary variable indicating if the patient is readmitted
# or not); here roughly 5% of patients will be readmitted.
endpoint = rand( num_patients ) .>= 0.95

# Generate predictions.
# Note that the predictions below are random; here, we would read the predictions from a CSV or
# some other source which would contain the predictions that our model would make on each patient.
# Note that as defined, predictions on the readmissions (endpoint = 1) are Uniform(0,0.9)
# random variables, and on the non-readmissions (endpoint = 0) are Uniform(0,0.4).
predict_0 = zeros(num_patients)
predict_0[endpoint] = 0.9*rand( sum(endpoint))
predict_0[!endpoint] = 0.4*rand( sum(!endpoint))

# predict_1 and predict_2 are arrays of predictions for the patients for 1 day, and 2 or more days after the completion of surgery,
# respectively. Here, we will just us the same predictions as the ones from the day of surgery completion. However,
# we can put different values in these arrays in case we wish to allow for the predictions to be updated during the patient's stay.
predict_1 = copy(predict_0)
predict_2 = copy(predict_1)

HCUP_cost = 14000 * ones(length(endpoint))

day_list = collect(1:num_days)

dayofweek = repeat( collect(1:7), 30)

# Assume Monday (=1), Wednesday (=3) schedule with 8 patients per day.
visit_days = [1,3]
capacity_per_day = 8

# Assume that the provider is *not* restricted to only seeing
# patients on the day they are being discharged.
atdischarge = false

# Run the simulation
PS, RA, ERC  = fin_simulate(day_of_entry, rday_of_surgery_completion, rday_of_discharge, predict_0, predict_1, predict_2, endpoint, HCUP_cost, dayofweek, day_list, visit_days, capacity_per_day, atdischarge)

# Calculate the provider cost.
provider_wage_per_hour = 75.0
PC = ceil( num_days / 7) * capacity_per_day * length(visit_days) * provider_wage_per_hour
if ( length(visit_days) * capacity_per_day >= 20)
	PC *= 1.25
end

# Calculate expected readmission cost savings, and expected net cost savings.
effectiveness_constant = 0.10
ERCS = effectiveness_constant * ERC
ENCS = ERCS - PC

# Display all the metrics:
@show PS, RA, ERC, ERCS, PC, ENCS

# Example of bootstrapping
Random.seed!(200)

boot_samples = Array{Int64,1}[]
nB = 1000
for b in 1:nB
	push!(boot_samples, rand(1:length(endpoint), length(endpoint)))
end

PS_vec = zeros(nB)
RA_vec = zeros(nB)
ERC_vec = zeros(nB)

for b in 1:nB
	println("bootstrap example; b = ", b)
	bs = boot_samples[b]
	PS, RA, ERC = fin_simulate(day_of_entry[bs], rday_of_surgery_completion[bs], rday_of_discharge[bs], predict_0[bs], predict_1[bs], predict_2[bs], endpoint[bs], HCUP_cost[bs], dayofweek, day_list,
							visit_days, capacity_per_day, atdischarge)
	PS_vec[b] = PS
	RA_vec[b] = RA
	ERC_vec[b] = ERC
end

PS_025, PS_975 = quantile(PS_vec, [0.025, 0.975])
RA_025, RA_975 = quantile(RA_vec, [0.025, 0.975])
ERC_025, ERC_975 = quantile(ERC_vec, [0.025, 0.975])

max_stars_repo_name: vvmisic/finsim-code

max_issues_repo_name: vvmisic/finsim-code

max_forks_repo_name: vvmisic/finsim-code


text:
# phylo/dating.jl
# ==================
#
# Types and methods for computing coalescence times between two sequences.
#
# Part of the Bio.Phylo module.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Dating

using Distributions, Roots

export DatingEstimate, coaltime, SimpleEstimate, SpeedDating, SDResult, upper,
       middle, lower

"""
An abstract dating estimate type.

Types inheriting from DatingEstimate should have the following methods defined:

* lower: Returns the lower bound of the dating estimate.
* upper: Returns the upper bound of the dating estimate.

They may define their own show methods.
"""
abstract DatingEstimate
abstract DatingMethod

Base.print(io::IO, de::DatingEstimate) = println(io, "$(lower(de)) ... $(upper(de))")
Base.show(io::IO, de::DatingEstimate) = println(io, "Coalesence time estimate between two sequences:\nt lies between $(lower(de)) and $(upper(de))")
Base.in(val::Float64, de::DatingEstimate) = val <= upper(de) && val >= lower(de)

"""
A very simple expected divergence time estimate.
Assumes a strict molecult clock and that the divergence time is equal to

\$ t = d / (2\\mu) \$

Where \$d\$ is the evolutionary distance computed for two aligned sequences,
and \$\\mu\$ is the substitution rate.
"""
immutable SimpleEstimate <: DatingMethod end
restype(::Type{SimpleEstimate}) = Float64

"""
    coaltime(N::Int, K::Int, µ::Float64, ::Type{SimpleEstimate})

Compute the coalescence time between two sequences by the `SimpleEstimate`
method.

`N` is the length of the two aligned sequences, `K` is the number of mutations
and `µ` is the assumed mutation rate.

# Examples
```julia
coaltime(50, 17, 10e-9, SimpleEstimate)
```
"""
@inline function coaltime(N::Int, K::Int, µ::Float64, ::Type{SimpleEstimate})
    @assert N >= K >= 0 error("Condition: N >= K >= 0, not met.")
    return (K / N) / (2 * µ)
end

"""
    coaltime(N::Int, d::Float64, µ::Float64, ::Type{SimpleEstimate})

Compute the coalescence time between two sequences by the `SimpleEstimate`
method.

`N` is the length of the two aligned sequences, `d` is the evolutionary distance
and `µ` is the assumed mutation rate.

# Examples
```julia
coaltime(50, 0.12, 10e-9, SimpleEstimate)
```
"""
@inline function coaltime(N::Int, d::Float64, µ::Float64, ::Type{SimpleEstimate})
    @assert N >= d >= 0 error("Genetic distance `d` must be a value between 1 and 0.")
    return d / (2 * µ)
end

"""
`SpeedDate` is the name given to a method of estimating a divergence time between
two DNA sequence regions that was first implemented in the R package
HybridCheck in order to date regions of introgression in large sequence
contigs.

The coalescence time is estimated using the number of mutations that have
occurred between two aligned sequences. The calculation uses a strict molecular
clock which assumes a constant substitution rate, both through time and across
taxa.
Modelling the mutation accumulation process as a Bernoulli trial, the
probability of observing  \$k\$  or fewer mutations between two sequences of
length \$n\$ can be given as:

\$ Pr(X \\le k) = \\sum_{i=0}^{\\lfloor k \\rfloor} \\binom{n}{i} p^i (1 - p)^{n-i} \$

Where \$p\$ is the probability of observing a single mutation between the two
aligned sequences.
The value of \$p\$ depends on two key factors: the substitution rate and the
coalescence time.
If you assume a molecular clock, whereby two DNA sequences are both accumulating
mutations at a rate \$\\mu\$ for \$t\$ generations, then you may define
\$p = 2\\mu t\$.

Using these assumptions, the SpeedDate method finds the root of the following
formula for \$Pr(X \\le k) = 0.05\$, \$0.5\$, and \$0.95\$, and then divides the
three answers by twice the assumed substitution rate.

\$ f(n, k, 2\\mu t, Pr(X \\le k) = \\left( \\sum_{i=0}^{\\lfloor k \\rfloor} \\binom{n}{i} {2\\mu t}^i (1 - 2\\mu t)^{n-i}   \\right) - Pr(X \\le k) \$

This results in an upper, middle, and lower estimate of the coalescence time
\$t\$ of the two sequences (expressed as the number of generations).
"""
immutable SpeedDating <: DatingMethod end
restype(::Type{SpeedDating}) = SDResult

"""
SDResult is a simple datatype to store the results of the coaltime
function.

coaltime estimates date range that may be considered a 95% confidence range in
which the true coalsescence time between two sequences lies.

The type store the 5%, 50%, and 95% values for the range.
"""
immutable SDResult <: DatingEstimate
    lower::Float64
    middle::Float64
    upper::Float64
end

lower(x::Float64) = x
lower(x::SDResult) = x.lower
middle(x::Float64) = x
middle(x::SDResult) = x.middle
upper(c::Float64) = x
upper(x::SDResult) = x.upper

function Base.show(io::IO, de::SDResult)
    println(io, "Coalescence time estimate:\n5%: $(de.lower), 95%: $(de.upper)")
end

@inline function binomzero(p0::Float64, N::Int, B::Int)
    f(p::Float64) = cdf(Binomial(N, p), B) - p0
    return fzero(f, 0, 1)
end

"""
    coaltime(N::Int, K::Int, µ::Float64, ::Type{SpeedDating})

Compute the coalescence time between two sequences by modelling the process of
mutation accumulation between two sequences as a bernoulli process.

This method was first described in the paper:
Ward, B. J., & van Oosterhout, C. (2016). Hybridcheck: Software for the rapid detection, visualization and dating of recombinant regions in genome sequence data. Molecular Ecology Resources, 16(2), 534–539.

`N` is the length of the two aligned sequences, `K` is the number of mutations
and `µ` is the assumed mutation rate.

# Examples
```julia
coaltime(50, 17, 10e-9, SpeedDating)
```
"""
@inline function coaltime(N::Int, K::Int, µ::Float64, ::Type{SpeedDating})
    @assert N >= K >= 0 error("Condition: N >= K >= 0, not met.")
    div = 2 * µ
    ninetyfive = ceil(binomzero(0.05, N, K) / div)
    fifty = ceil(binomzero(0.5, N, K) / div)
    five = ceil(binomzero(0.95, N, K) / div)
    return SDResult(five, fifty, ninetyfive)
end

"""
    coaltime(N::Int, p::Float64, µ::Float64, ::Type{SpeedDating})

Compute the coalescence time between two sequences by modelling the process of
mutation accumulation between two sequences as a bernoulli process.

This method was first described in the paper:
Ward, B. J., & van Oosterhout, C. (2016). Hybridcheck: Software for the rapid
detection, visualization and dating of recombinant regions in genome sequence
data. Molecular Ecology Resources, 16(2), 534–539.

`N` is the length of the two aligned sequences, `d` is the evolutionary distance
and `µ` is the assumed mutation rate.

# Examples
```julia
coaltime(50, 17, 10e-9, SpeedDating)
```
"""
@inline function coaltime(N::Int, d::Float64, µ::Float64, ::Type{SpeedDating})
    @assert 1 >= d >= 0 error("Genetic distance `d` must be a value between 1 and 0.")
    K = convert(Int, ceil(d * N))
    return coaltime(N, K, µ, SpeedDating)
end

"""
    coaltime{M<:DatingMethod,R<:Real}(N::Int, arr::AbstractArray{R}, µ::Float64, ::Type{M})

Compute the coalescence time between two sequences by modelling the process of
mutation accumulation between two sequences as a bernoulli process.

This method was first described in the paper:
Ward, B. J., & van Oosterhout, C. (2016). Hybridcheck: Software for the rapid detection, visualization and dating of recombinant regions in genome sequence data. Molecular Ecology Resources, 16(2), 534–539.

In this specific method, `N` is the length of the aligned sequences,
`arr` is an array of either mutation counts (Integers) or genetic distances
(Floats), `µ` is the assumed mutation rate.

# Examples
```julia
coaltime(50, [17, 20, 10, 7], 10e-9, SpeedDating)
coaltime(50, [0.01, 0.09, 0.12, 0.20], 10e-9, SpeedDating)
```
"""
@inline function coaltime{M<:DatingMethod,R<:Real}(N::Int,
                                           arr::AbstractArray{R},
                                           µ::Float64,
                                           ::Type{M})
    Ts = similar(arr, restype(M))
    @inbounds for i in eachindex(arr)
        Ts[i] = coaltime(N, arr[i], µ, M)
    end
    return Ts
end

end


text:
using Test
import AoC2019

input_data = AoC2019.read_chars("./input16.txt")

function generate_multipliers(phase, length_of_input)
  base_pattern = Int8[0,1,0,-1]
  multipliers = zeros(Int8, length_of_input)
  stretched_pattern = Iterators.flatten(map(x -> Iterators.repeated(x, phase), base_pattern))
  whole_pattern = Iterators.cycle(stretched_pattern)
  #skip first
  for (i, v) in enumerate(whole_pattern)
    i == 1 && continue
    multipliers[i-1] = v
    i > length_of_input && break
  end
  return multipliers
end

@test generate_multipliers(1, 8) == [1,0,-1,0,1,0,-1,0]
@test generate_multipliers(2, 8) == [0,1,1,0,0,-1,-1,0]
@test generate_multipliers(8, 8) == [0,0,0,0,0,0,0,1]

function get_digit(input, num_entry::Int64)

  base_pattern = Int8[0,1,0,-1]

  digit::Int64 = 0
  #println(multipliers)
  for i in 1:num_entry
    digit += sum(input[i-1+num_entry:4*num_entry:end])
    digit -= sum(input[i-1+num_entry + 2*num_entry:4*num_entry:end])

  end
  return abs(digit % 10)
end

function fft_step(input)
  new_digits = Array{Int8, 1}(undef, length(input))
  for (i, digit) in enumerate(input)
    new_digits[i] = get_digit(input, i)
  end
  return new_digits
end

function my_fft(initial_input::Array{Int8, 1}, num_phases::Int64)
  # generate multipliers

  for i in 1:num_phases
    initial_input = fft_step(initial_input)
  end
  return initial_input
end

function my_fft_big(initial_input::Array{Int8, 1}, num_phases::Int64, replication_factor::Int64)
  my_big_input = collect(Iterators.flatten(Iterators.repeated(initial_input, replication_factor)))
  for i in 1:num_phases
    println("Iteration: $i")
    my_big_input = fft_step(my_big_input)
    println(my_big_input)
  end
  return my_big_input
end

test_input1 = Int8[1,2,3,4,5,6,7,8]
@test fft_step(test_input1) == Int8[4,8,2,2,6,1,5,8]
@test my_fft(test_input1, 4) == Int8[0,1,0,2,9,4,9,8]
test_input2 = Int8[8,0,8,7,1,2,2,4,5,8,5,9,1,4,5,4,6,6,1,9,0,8,3,2,1,8,6,4,5,5,9,5]
test_input3 = Int8[1,9,6,1,7,8,0,4,2,0,7,2,0,2,2,0,9,1,4,4,9,1,6,0,4,4,1,8,9,9,1,7]
@test my_fft(test_input2, 100)[1:8] == [2,4,1,7,6,1,7,6]
my_fft(input_data, 100)
#my_fft_big(input_data, 100, 1)
#@assert my_fft(input_data, 100) == my_fft_big(input_data, 100,1)

#test_input4 = Int8[0,3,0,3,6,7,3,2,5,7,7,2,1,2,9,4,4,0,6,3,4,9,1,5,6,5,4,7,4,6,6,4]
#test_input5 = Int8[1,2,3,4,5,6]
#@assert my_fft_big(input_data, 1, 2)[303673+1:313673+8] == [8,4,4,6,2,0,2,6]

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
using NCDatasets
using Formatting
using PyPlot

#
# Global Diagnose
#
# X v.s. latitude
#
# - Heat transport
# - Precipitation
#
# X v.s. time
#
# - Jet stream position
# - Precipitation
#
casename = "SMARTSLAB_SOM_lowres"
folder = joinpath("atm", "hist")
getFn0 = (y, m) -> joinpath(folder, format("$casename.cam.h0.00{:02d}-{:02d}.nc", y, m))

y_rngs = [(1, 20), (1,10), (11, 20)]

for (y_beg, y_end) in y_rngs

    local mht = nothing

    for y = y_beg:y_end, m = 1:12

        ds = Dataset(fn, "r")

        if mht == nothing
            mht = zeros(Float64, size(ds["VT"])...)
        end

        mht

        fn = getFn(y, m)

        close(ds)
    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
struct BeamSection
    id::String
    hid::Int
    A::Float64
    I₂::Float64
    I₃::Float64
    J::Float64
    As₂::Float64
    As₃::Float64
    W₂::Float64
    W₃::Float64
    sec_type::SectionType
    sizes::Vector
    BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃,sec_type=SectionType(0),sizes=[])=new(string(id),hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃,sec_type,sizes)
end

# function BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃)
#     id=string(id)
# end

function RectangleSection(id,hid,h,b)::BeamSection
    id=string(id)

    A=h*b
    I₃=b*h^3/12
    I₂=h*b^3/12

    bb,aa=sort([h,b])
    J=aa*bb^3*(1/3-0.21*bb/aa*(1-bb^4/12/aa^4))

    As₂=A#wrong
    As₃=A#wrong
    W₃=I₃/h*2
    W₂=I₂/b*2

    BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃,Enums.RECTANGLE,[h,b])
end

function HSection(id,hid,h,b,tw,tf)::BeamSection
    id=string(id)

    A=b*tf*2+tw*(h-2*tf)
    J=(b*tf^3*2+(h-tf)*tw^3)/3 #should be confirm!!!!!!!!!!
    I₃=b*h^3/12-(b-tw)*(h-2*tf)^3/12
    I₂=2*tf*b^3/12+(h-2*tf)*tw^3/12
    As₂=tw*(h-2tf)#wrong
    As₃=2*b*tf#wrong
    W₃=I₃/h*2
    W₂=I₂/b*2
    BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃,Enums.HSECTION,[h,b,tw,tf])
end

function ISection(id,hid,h,b1,b2,tw,tf1,tf2)::BeamSection
    id=string(id)

    hw=h-tf1-tf2
    A=b1*tf1+b2*tf2+tw*hw
    y0=(b1*tf1*(h-tf1/2)+b2*tf2*tf2/2+hw*tw*(hw/2+tf2))/A

    J=(b1*tf1^3+b2*tf2^3+(h-tf1/2-tf2/2)*tw^3)/3 #should be confirm!!!!!!!!!!

    I₃=tw*hw^3/12
    I₃+=b1*tf1^3/12+b1*tf1*(hw/2+tf1/2)^2
    I₃+=b2*tf2^3/12+b2*tf2*(hw/2+tf2/2)^2
    I₃-=A*(y0-h/2)^2

    I₂=b1^3*tf1/12+b2^3*tf2/12+tw^3*hw/12

    As₂=tw*hw#wrong
    As₃=b1*tf1+b2*tf2#wrong
    W₃=I₃/max(y0,h-y0)
    W₂=I₂/max(b1/2,b2/2)

    BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃,Enums.I_SECTION,[h,b1,b2,tw,tf1,tf2])
end

function BoxSection(id,hid,h,b,tw,tf)::BeamSection
    id=string(id)

    A=h*b-(h-2*tf)*(b-2*tw)
    J=(2*tw*(h-tf)/tw+2*tf*(b-tw)/tf)*2*A
    I₃=b*h^3/12-(b-2*tw)*(h-2*tf)^3/12
    I₂=h*b^3/12-(h-2*tf)*(b-2*tw)^3/12

    As₂=2*h*tw#wrong
    As₃=2*b*tf#wrong
    W₃=I₃/h*2
    W₂=I₂/b*2
    BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃,Enums.BOX,[h,b,tw,tf])
end

function PipeSection(id,hid,d,t)::BeamSection
    id=string(id)

    A=π*d^2/4-π*(d-2*t)^2/4
    J=π*(d-t)/t*2*A
    I₃=π*d^4/64*(1-((d-2*t)/d)^4)
    I₂=I₃

    As₂=A#wrong
    As₃=A#wrong
    W₃=I₃/d*2
    W₂=W₃
    r=d/2
    J=2/3*π*r*t^3
    BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃,Enums.PIPE,[d,t])
end

function CircleSection(id,hid,d)::BeamSection
    id=string(id)

    A=π*d^2/4
    J=π*d^4/32
    I₃=π*d^4/64
    I₂=I₃

    As₂=A#wrong
    As₃=A#wrong
    W₃=I₃/d*2
    W₂=W₃
    J=π*d^4/32
    BeamSection(id,hid,A,I₂,I₃,J,As₂,As₃,W₂,W₃,Enums.CIRCLE,[d])
end

# end


text:
#jacobian2by2.jl
# Code 8.1 of Random Eigenvalues by Alan Edelman

#Experiment:    Compute the Jacobian of a 2x2 matrix function
#Comment:       Symbolic tools are not perfect. The author
#               exercised care in choosing the variables

warn("Julia doesn't currently have a good symbolic math package")
warn("This may use features of PyCall and/or SymPy which are not publicly available yet")
using PyCall
pyinitialize("/opt/local/bin/python2")
using SymPy
jacobian(Y,X) = Sym(Y[:reshape](1,length(Y))[:jacobian](X[:reshape](1,length(X))))

p,q,r,s, a,b,c,d, t, e1,e2 = Sym(:p,:q,:r,:s, :a,:b,:c,:d, :t, :e1,:e2)
X = SymMatrix([p q; r s])
A = SymMatrix([a b; c d])

Y = X^2; J = jacobian(Y,X); JAC_square = Sym(J[:det]()[:factor]()); println("square:\n", J, "\n\n", JAC_square, "\n\n")
Y = X^3; J = jacobian(Y,X); JAC_cube = Sym(J[:det]()[:factor]()); println("cube:\n", J, "\n\n", JAC_cube, "\n\n")
Y = X[:inv](); J = jacobian(Y,X); JAC_inv = Sym(J[:det]()[:factor]()); println("inv:\n", J, "\n\n", JAC_inv, "\n\n")
Y = A*X; J = jacobian(Y,X); JAC_linear = Sym(J[:det]()[:factor]()); println("linear:\n", J, "\n\n", JAC_linear, "\n\n")
Y = SymMatrix([p q; r/p X[:det]()/p]); J = jacobian(Y,X); JAC_lu = Sym(J[:det]()[:factor]()); println("lu:\n", J, "\n\n", JAC_lu, "\n\n")

x = SymMatrix([p, s, r])
y = SymMatrix([sqrt(p), sqrt(s), r/(sqrt(p)*sqrt(s))])
J = jacobian(y,x)
JAC_DMD = Sym(J[:det]()[:factor]())
println("DMD:\n", J, "\n\n", JAC_DMD, "\n\n")

x = SymMatrix([p, s])
y = SymMatrix([sqrt(p^2 + s^2), atan(s/p)])
J = jacobian(y,x)
JAC_notrace = Sym(J[:det]()[:factor]())
println("notrace:\n", J, "\n\n", JAC_notrace, "\n\n")

Q = SymMatrix([cos(t)  -sin(t); sin(t)  cos(t)])
D = SymMatrix([e1 0; 0 e2])
Y = Q*D*Q[:transpose]()
y = SymMatrix([Y[1,1], Y[2,2], Y[1,2]])
x = SymMatrix([t, e1, e2])
J = jacobian(Y,X)
JAC_symeig = Sym(J[:det]()[:factor]())
println("symeig:\n", J, "\n\n", JAC_symeig, "\n\n")

X = SymMatrix([p s; s r])
Y = Sym(A[:transpose]())*X*A
y = SymMatrix([Y[1,1], Y[2,2], Y[1,2]])
x = SymMatrix([p, r, s])
J = jacobian(Y,X)
JAC_symcong = Sym(J[:det]()[:factor]())
println("symcong:\n", J, "\n\n", JAC_symcong, "\n\n")


text:
module DNABarcode

export
    HammingCode,
    genbarcodes,
    demultiplex

import Bio.Seq: DNAKmer, DNASequence, DNA_A, mismatches, isambiguous

include("abstractcode.jl")
include("hamming.jl")
include("sequence_levenshtein.jl")

end # module


text:
using Literate

"""
    notebook_filter(str)

Filter to correctly render inline math (two backquotes) in notebooks.
"""
function notebook_filter(str)
    re = r"(?<!`)``(?!`)"  # Two backquotes not preceded by nor followed by another
    replace(str, re => "\$")
end

"""
    notebook(input_file)

Create notebooks from `input_file`, applying needed filters and write output to
the root folder of our typical Jupyter session.
"""
function notebook(input_file)
    output_dir = "."
    Literate.notebook(input_file, output_dir; preprocess=notebook_filter)
end

nothing;


text:
# Get root logger.
# Nothing appears as we haven't set any config on any loggers.
root_logger = get_logger()
println(root_logger)

# This is also root logger.
println(get_logger(""))

# Set root config.
# It inserts a handler that outputs message to `STDERR`.
basic_config(MiniLogging.INFO, date_format="")

# It changes the root logger level.
get_logger("")
println(root_logger)

@warn(root_logger, "Hello", " world")

# Get a logger.
logger = get_logger("a.b")
println(logger)

# Since the level of `logger` is unset, it inherits its nearest ancestor's level.
# Its effective level is `INFO` (from `root_logger`) now.
@info(logger, "Hello", " world")

# Since `DEBUG < INFO`, no message is generated.
# Note the argument expressions are not evaluated in this case to increase performance.
@debug(logger, "Hello", " world", error("no error"))

# Explicitly set the level.
# The error is triggered.
logger.level = MiniLogging.DEBUG
try
    @debug(logger, "Hello", " world", error("has error"))
catch e
    println(e)
end

# Get a child logger.
logger2 = get_logger("a.b.c")
println(logger2)

# Its effective level now is `DEBUG` (from `logger`) now.
@debug(logger2, "Hello", " world")

# Add A New Level
MiniLogging.define_new_level(:trace, 25, :yellow)
@trace(logger, "Hello", " world")

max_stars_repo_name: colinfang/MiniLogging

max_issues_repo_name: colinfang/MiniLogging

max_forks_repo_name: colinfang/MiniLogging


text:
function collatz(n::Integer)
    n <= 0 && throw(DomainError())
    iseven(n) ? div(n, 2) : 3n + 1
end

function collatz_steps(n::Integer)
    i = 0
    while n != 1
        n = collatz(n)
        i += 1
    end

    i
end


text:
function unit_test()
    helpers_test_unit()
    windowing_test_unit()
    unassigned_test_unit()
    week_test_unit()

    Logging.info("Scheduler Unit Tests passing. Good job.")
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
using Documenter
using UnitfulEquivalences

DocMeta.setdocmeta!(UnitfulEquivalences, :DocTestSetup, :(using Unitful, UnitfulEquivalences))

makedocs(
    sitename = "UnitfulEquivalences.jl",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [UnitfulEquivalences],
    pages = [
             "Home" => "index.md",
             "Converting between equivalent quantities" => "usage.md",
             "Defining new equivalences" => "newequivalences.md"
            ]
)

deploydocs(
    repo = "github.com/sostock/UnitfulEquivalences.jl.git",
    devbranch = "main"
)


text:
# This file contains necessary code to compute weights according to a sum
# decomposition

export reachable_set, propagate_reachable_set, channel_capacities_by_target,
  channel_capacities_by_target_time, channel_capacities_mcts

# Compute reachable set on the grid for n steps as a set of states
function reachable_set(grid::Grid, state::State, steps::Integer)
  states = Set([state])

  for ii = 1:steps
    states = propagate_reachable_set(grid, states)
  end

  states
end
# Helper function for computation of reachable sets
function propagate_reachable_set(grid::Grid, states)
  union(map(x->Set(neighbors(grid, x)), collect(states))...)
end

# Computes reachable sets for each step in a horizon and returns the result in a
# named tuple to avoid ambiguity with the form
# (step=<step>, states=Set{State})
struct ReachableAtTime
  step::Int64
  states::Vector{State}
end
function reachable_sets(grid::Grid, state::State, horizon::Integer)
  sets = Vector{ReachableAtTime}(undef, horizon)

  for step in 1:horizon
    local states

    if step == 1
      states = propagate_reachable_set(grid, Set([state]))
    else
      states = propagate_reachable_set(grid, sets[step-1].states)
    end

    sets[step] = ReachableAtTime(step, collect(states))
  end

  sets
end

# We will define a partition matroid for an array of ReachableAtTime
#
# To do so, we provide access to an iterable object that represents the
# collection of set actions in a given block of the partition matroid
get_block(x::ReachableAtTime) = x.states
# get element will determine what goes into the objective
function get_element(partition_matroid::Vector{ReachableAtTime}, x)
  set = partition_matroid[x[1]]

  TimedObservation(set.step, set.states[x[2]])
end

#
# Channel capacity bounds
#
# The use of the sum decomposition of the objective relies on computation of
# upper-bounds of the individual summands. In the case of mutual information
# objectives, these bounds are channel capacities
#

#
# Compute an individual channel capacity bound between a robot and a target,
# optionally at a specific time-step
#

function channel_capacity_by_target_time(p::MultiRobotTargetTrackingProblem,
                                         filter::Filter;
                                         step, robot_index)
  configs = p.configs

  # Note: observations is a vector of TimedObservation
  function objective(observations)
    finite_horizon_information(configs.grid, filter, configs.sensor,
                               observations, step;
                               num_samples=
                                configs.objective_information_samples,
                               entropy_only_at_end=true
                              ).reward
  end

  # Due to the objective design, states after the given time-step do not matter
  partition_matroid = reachable_sets(configs.grid,
                                     p.partition_matroid[robot_index],
                                     step)

  # Construct and solve the optimization problem
  problem = ExplicitPartitionProblem(objective, partition_matroid)
  solution = solve_sequential(problem)

  # Assuming that the sampling of the objective is adequate
  # opt <= 2 * (solution.value + epsilon)
  # and epsilon ~= 0
  #
  # Hence we treat the following as an upper bound on the weight.
  #
  # Note: Other, tighter bounds are possible. These include:
  #   * Online bounds based on submodularity (Sum of greatest marginals in each
  #     block)
  #   * Possibly some channel capacity bounds, but that is a different story
  2 * solution.value
end

function channel_capacity_by_target(p::MultiRobotTargetTrackingProblem,
                                    filter::Filter; robot_index)
  configs = p.configs

  function objective(observations)
    finite_horizon_information(configs.grid, filter, configs.sensor,
                               observations, configs.horizon;
                               num_samples=
                                configs.objective_information_samples,
                               entropy_only_at_end=false
                              ).reward
  end

  partition_matroid = reachable_sets(configs.grid,
                                     p.partition_matroid[robot_index],
                                     configs.horizon)

  # Construct and solve the optimization problem
  problem = ExplicitPartitionProblem(objective, partition_matroid)
  solution = solve_sequential(problem)

  # See discussion in the target/time implementation
  2 * solution.value
end

# Partition matroids turned out to be super slow as I probably could have
# guessed. Using MCTS is a bit tenuous, but I should be able to get results in a
# small fraction of the time I would require otherwise
#
# On the plus side, this should be a relatively easy problem for MCTS to solve
function channel_capacity_mcts(p::MultiRobotTargetTrackingProblem,
                               filter::AbstractFilter; robot_index)
  configs=p.configs

  # Construct a solver for the individual filter
  problem = SingleRobotTargetTrackingProblem(configs.grid, configs.sensor,
                                           configs.horizon,
                                           [filter],
                                           num_information_samples=
                                             configs.solver_information_samples)
  function objective(observations)
    finite_horizon_information(configs.grid, filter, configs.sensor,
                               observations, configs.horizon;
                               num_samples=
                                configs.objective_information_samples,
                               entropy_only_at_end=false
                              ).reward
  end

  state = p.partition_matroid[robot_index]

  trajectory = solve_single_robot(problem, state,
                                  n_iterations=
                                    configs.solver_iterations).trajectory

  # Sample again to get an accurate estimate of the reward
  #
  # The tuple is a hack to put this in a format tha the multi-robot code can
  # recognize
  objective([trajectory])
end

#
# Compute all channel capacity bounds (by target and or time-step) for a given
# robot and store in an array
#

function channel_capacities_by_target_time(p::MultiRobotTargetTrackingProblem,
                                           robot_index::Integer,
                                           target_in_range =
                                           (robot_index, target_index) -> true
                                          )
  configs = p.configs

  map(product(enumerate(p.target_filters), 1:configs.horizon)) do
    ((target_index, filter), step)

    if target_in_range(robot_index, target_index)
      channel_capacity_by_target_time(p, filter, step=step,
                                      robot_index=robot_index)
    else
      0.0
    end
  end
end

function channel_capacities_by_target(p::MultiRobotTargetTrackingProblem,
                                      robot_index::Integer,
                                      target_in_range = (robot_index,
                                                         target_index) -> true
                                     )
  map(enumerate(p.target_filters)) do (target_index, filter)
    if target_in_range(robot_index, target_index)
      channel_capacity_by_target(p, filter,robot_index=robot_index)
    else
      0.0
    end
  end
end

function channel_capacities_mcts(p::MultiRobotTargetTrackingProblem,
                                 robot_index::Integer,
                                 target_in_range = (robot_index,
                                                    target_index) -> true
                                )
  map(enumerate(p.target_filters)) do (target_index, filter)
    if target_in_range(robot_index, target_index)
      channel_capacity_mcts(p, filter, robot_index=robot_index)
    else
      0.0
    end
  end
end

# Compute a weight based on channel capacities from the sum decomposition
#
# Note: here the weight corresponds to the maximum reduction of one robot's
# marginal gain given the other
#
# Here, we take advantage of the fact that the marginal gain with respect to a
# target must always be positive. As such, marginal can be reduced by at most
# the lesser channel capacity. Summing produces a bound over all components of
# the decomposition
compute_weight(a::Array{Float64}, b::Array{Float64}) = sum(min.(a, b))

function compute_weight_matrix(p::MultiRobotTargetTrackingProblem;
                               channel_capacity_method=channel_capacities_mcts,
                               robot_target_range_limit=Inf,
                               threaded = false
                              )
  local capacities

  n = get_num_agents(p)

  function range_test(robot_index, target_index)
    target_mean_distance(p, robot_index,
                         target_index) < robot_target_range_limit
  end

  # Compute channel capacities for each robot
  channel_capacity(x) = channel_capacity_method(p, x, range_test)

  if threaded
    capacities = thread_map(channel_capacity, 1:n, Vector{Float64})
  else
    capacities = map(channel_capacity, 1:n)
  end

  weights = zeros(n,n)

  for ii in 2:n, jj in 1:ii-1
    w = compute_weight(capacities[ii], capacities[jj])

    weights[ii, jj] = w
    weights[jj, ii] = w
  end

  weights
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
"""Utilities for working with Bravais lattices"""

"""
    lattice_params(lat_vecs::Mat3)

Compute the lattice parameters ``(a, b, c, α, β, γ)`` from a set of lattice vectors,
 which form the columns of `lat_vecs`.
"""
function lattice_params(lat_vecs::Mat3) :: NTuple{6, Float64}
    v1, v2, v3 = eachcol(lat_vecs)
    a, b, c = norm(v1), norm(v2), norm(v3)
    α = acosd((v2 ⋅ v3) / (b * c))
    β = acosd((v1 ⋅ v3) / (a * c))
    γ = acosd((v1 ⋅ v2) / (a * b))
    return (a, b, c, α, β, γ)
end

"""
    lattice_vectors(a, b, c, α, β, γ)

Return the lattice vectors, as columns of the ``3×3`` output matrix, that
correspond to the conventional unit cell defined by the lattice constants ``(a,
b, c)`` and the angles ``(α, β, γ)``.
"""
function lattice_vectors(a, b, c, α, β, γ) :: Mat3
    @assert all(0 < x < 180 for x in (α, β, γ))

    sγ, cγ = sind(γ), cosd(γ)
    cβ, cα = cosd(β), cosd(α)
    v1 = Vec3(a, 0, 0)
    v2 = Vec3(b * cγ, b * sγ, 0)
    v3x = c * cβ
    v3y = c / sγ * (cα - cβ * cγ)
    v3z = c / sγ * √(sγ^2 - cα^2 - cβ^2 + 2 * cα * cβ * cγ)
    v3 = Vec3(v3x, v3y, v3z)

    @assert norm(v1) ≈ a
    @assert norm(v2) ≈ b
    @assert norm(v3) ≈ c
    @assert acosd(v1⋅v2 / (a*b)) ≈ γ
    @assert acosd(v1⋅v3 / (a*c)) ≈ β
    @assert acosd(v2⋅v3 / (b*c)) ≈ α

    return [v1 v2 v3]
end

function is_standard_form(lat_vecs::Mat3)
    lat_params = lattice_params(lat_vecs)
    conventional_lat_vecs = lattice_vectors(lat_params...)
    return lat_vecs ≈ conventional_lat_vecs
end

# Return true if lattice vectors are compatible with a monoclinic space group
# "setting" for a Hall number.
function is_compatible_monoclinic_cell(lat_vecs, hall_number)
    @assert cell_type(hall_number) == monoclinic

    # Special handling of monoclinic space groups. There are three possible
    # conventions for the unit cell, depending on which of α, β, or γ is
    # special.
    _, _, _, α, β, γ = lattice_params(lat_vecs)
    choice = Spglib.get_spacegroup_type(hall_number).choice
    x = first(replace(choice, "-" => ""))
    if x == 'a'
        return β≈90 && γ≈90
    elseif x == 'b'
        return α≈90 && γ≈90
    elseif x == 'c'
        return α≈90 && β≈90
    else
        error()
    end
end

"""
    CellType

An enumeration over the different types of 3D Bravais unit cells.
"""
@enum CellType begin
    triclinic
    monoclinic
    orthorhombic
    tetragonal
    # Rhombohedral is a special case. It is a lattice type (a=b=c, α=β=γ) but
    # not a spacegroup type. Trigonal space groups are conventionally described
    # using either hexagonal or rhombohedral lattices.
    rhombohedral
    hexagonal
    cubic
end

"""
    cell_type(lat_vecs::Mat3)

Infer the `CellType` of a unit cell from its lattice vectors, i.e. the columns
of `lat_vecs`. Report an error if lattice vectors are not in conventional form.
"""
function cell_type(lat_vecs::Mat3)
    a, b, c, α, β, γ = lattice_params(lat_vecs)

    if !(lat_vecs ≈ lattice_vectors(a, b, c, α, β, γ))
        error("Lattice vectors are not in conventional form. Consider using `lattice_vectors(a, b, c, α, β, γ)`.")
    end

    if a ≈ b ≈ c
        if α ≈ β ≈ γ ≈ 90
            return cubic
        elseif α ≈ β ≈ γ
            return rhombohedral
        end
    end

    if α ≈ β ≈ γ ≈ 90
        if a ≈ b
            return tetragonal
        elseif b ≈ c || c ≈ a
            error("Found a nonconventional tetragonal unit cell. Use `lattice_vectors(a, a, c, 90, 90, 90)` instead.")
        else
            return orthorhombic
        end
    end

    if (a ≈ b && α ≈ β ≈ 90 && (γ ≈ 60 || γ ≈ 120)) ||
       (b ≈ c && β ≈ γ ≈ 90 && (α ≈ 60 || α ≈ 120)) ||
       (c ≈ a && γ ≈ α ≈ 90 && (β ≈ 60 || β ≈ 120))
        if γ ≈ 120
            return hexagonal
        else
            error("Found a nonconventional hexagonal unit cell. Use `lattice_vectors(a, a, c, 90, 90, 120)` instead.")
        end
    end

    # Accept any of three possible permutations for monoclinic unit cell
    if α ≈ β ≈ 90 || β ≈ γ ≈ 90 || α ≈ γ ≈ 90
        return monoclinic
    end

    return triclinic
end

"Return the standard cell convention for a given Hall number"
# Using the convention of spglib, listed at
# http://pmsl.planet.sci.kobe-u.ac.jp/~seto/?page_id=37
function cell_type(hall_number::Int)
    if 1 <= hall_number <= 2
        triclinic
    elseif 3 <= hall_number <= 107
        monoclinic
    elseif 108 <= hall_number <= 348
        orthorhombic
    elseif 349 <= hall_number <= 429
        tetragonal
    elseif 430 <= hall_number <= 461
        # The trigonal space groups require either rhombohedral or hexagonal
        # cells. The Hall numbers below have "setting" R.
        hall_number in [434, 437, 445, 451, 453, 459, 461] ? rhombohedral : hexagonal
    elseif 462 <= hall_number <= 488
        hexagonal
    elseif 489 <= hall_number <= 530
        cubic
    else
        error("Invalid Hall number $hall_number. Allowed range is 1..530")
    end
end

function all_compatible_cells(cell::CellType)
    if cell == triclinic
        [triclinic, monoclinic, orthorhombic, tetragonal, rhombohedral, hexagonal, cubic]
    elseif cell == monoclinic
        [monoclinic, orthorhombic, tetragonal, hexagonal, cubic]
    elseif cell == orthorhombic
        [orthorhombic, tetragonal, cubic]
    elseif cell == tetragonal
        [tetragonal, cubic]
    elseif cell == rhombohedral
        [rhombohedral]
    elseif cell == hexagonal
        [hexagonal]
    elseif cell == cubic
        [cubic]
    else
        error()
    end
end

function is_trigonal_symmetry(hall_number::Int)
    return 430 <= hall_number <= 461
end


text:
# This file was generated, do not modify it. # hide
model_info=DataFrame(ModelName=model_names,Accuracy=loss_acc,CrossEntropyLoss=loss_ce,F1Score=loss_f1);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE

using NodeNumbering

adjacency = Dict(
    1 => [3, 8, 9],
    2 => [3, 8, 7, 6],
    3 => [1, 2],
    4 => [8, 9],
    5 => [7, 8],
    6 => [2, 7],
    7 => [5, 2, 6],
    8 => [1, 2, 4, 5],
    9 => [1, 4])

@testset "node degrees" begin
    E = Dict(
        1 => 3,
        2 => 4,
        3 => 2,
        4 => 2,
        5 => 2,
        6 => 2,
        7 => 3,
        8 => 4,
        9 => 2)
    result = NodeNumbering.node_degrees(adjacency)
    @test result == E
end


text:
module structs
    using DataStructures

    include("io/internal.jl")
    include("io/output.jl")

    mutable struct MOEUM
        """
        An Tabular Struct for Massive Data Processing

        **Common methods**
        """

        ###parameters
        #arguments
        dat_sori::Array
        hol_sori::Array
        name::String
        #hidden parameters
        _sok_dict::DataStructures.OrderedDict
        _sok_arr::Array

        ###functions
        #internal functions
        select::Function
        describe::Function

        #modify
        insert::Function
        append::Function
        #equation::Function

        #query
        #where::Function

        #output
        to_string::Function
        to_dataframe::Function
        to_dict::Function
        to_csv::Function
        to_db::Function

        #init
        function MOEUM(;dat_sori = [], hol_sori = [], name = "moeum.MOEUM")
            #init instance
            instance = new(dat_sori, hol_sori, name)

            #=key_list = instance.dat_sori; sok_dict = Dict(dat => [] for dat in instance.dat_sori)
            for hol in instance.hol_sori
                for i in eachindex(instance.dat_sori)
                    append!(sok_dict[instance.dat_sori[i]], hol[i])
                end
            end=#
            instance._sok_dict, instance._sok_arr = internal.set_sok(instance.dat_sori, instance.hol_sori, instance.name)

            #instance._sok_dict = DataStructures.OrderedDict(key_list[i]=> ColumnSet(sok_dict[key_list[i]], i, key_list[i], instance.name) for i in eachindex(key_list))
            #instance._sok_arr = [RowSet(DataStructures.OrderedDict(instance.dat_sori[j] => instance.hol_sori[i][j] for j in eachindex(instance.dat_sori)), i, instance.name) for i in eachindex(instance.hol_sori)]

            #modify
            instance.insert = function(;dat_pos::Int = length(instance.dat_sori) + 1, dat_name::Any = "", hol_value::Array = [])
                return internal.insert(instance, dat_pos, dat_name, hol_value)
            end

            instance.append = function(;dat_name::Any = "", hol_value::Array = [])
                return internal.append(instance, dat_name, hol_value)
            end

            #internal functions
            instance.select = function(key)
                return internal.select(instance, key)
            end

            instance.describe = function()
                return internal.describe(instance)
            end

            #output functions
            instance.to_string = function()
                return output.to_string(instance)
            end

            instance.to_dataframe = function()
                return output.to_dataframe(instance)
            end

            instance.to_dict = function(;orientation::String = "Dict")
                return output.to_dict(instance, orientation)
            end

            instance.to_csv = function(;csv_path::String = joinpath(pwd(), string(instance.name, ".csv")))
                return output.to_csv(instance, csv_path)
            end

            return instance
        end
    end
    function Base.show(io::Base.IO, instance::MOEUM)
        println(io, instance.to_string())
    end

    mutable struct ColumnSet
        ###parameters
        #arguments
        data::Array
        #hidden parameters
        _data_type::String
        _parent_name::String
        _col_num::Int
        _col_name::String

        ###functions
        #descriptive statistics/common
        describe::Function
        count::Function
        #descriptive statistics/numeric
        sum::Function
        avg::Function
        std::Function
        min::Function
        max::Function
        #descriptive statistics/string
        #top::Function
        unique::Function
        #descriptive statistics/invalid
        _invalid::Function

        #eval equation
        equation::Function

        #init
        function ColumnSet(data = [], col_num = 0, col_name = "dat", parent_name = "moeum.MOEUM")
            instance = new(data) ; instance._parent_name, instance._col_num, instance._col_name = col_name = parent_name, col_num, col_name
            type_count = Dict("Numeric"=>0, "String"=>0)
            for i in eachindex(instance.data)
                item = string(instance.data[i])
                if !isnull(tryparse(Float64, item))
                    if contains(item, ".")
                        parse_type = Float64
                    else
                        parse_type = Int
                    end

                    instance.data[i] = parse(parse_type, item)
                    type_count["Numeric"] += 1
                else
                    type_count["String"] += 1
                end
            end

            #descriptive statistics/invalid
            instance._invalid = function()
                println(string("invalid for data-type is ", instance._data_type))
                return false
            end

            #set default invalid to functions
            instance.equation, instance.sum, instance.avg, instance.std, instance.min, instance.max, instance.unique, instance.describe = [instance._invalid for i in 1:8]
            if type_count["Numeric"] >= type_count["String"]
                instance._data_type = "Numeric"
                #instance.data = collect(float(item) for item in instance.data)

                #descriptive statistics/numeric
                parsed_data = [parse(Float64, string(item)) for item in instance.data]
                instance.sum = function()
                    return Base.sum(instance.data)
                end

                instance.avg = function()
                    return Base.mean(parsed_data)
                end

                instance.std = function()
                    return Base.std(parsed_data)
                end

                instance.max = function()
                    Base.maximum(instance.data)
                end

                instance.min = function()
                    Base.minimum(instance.data)
                end

                #eval equation
                instance.equation = function(equation_string)
                    for item in instance.data
                        eval(parse(equation_string))
                    end
                end
            else
                instance._data_type = "String"
                #instance.data = collect(string(item) for item in instance.data)

                #descriptive statistics/string
                #instance.top = function()
                #    Base.top(instance.data)
                #end

                instance.unique = function()
                    return Base.unique(instance.data)
                end
            end

            #descriptive statistics/common
            instance.describe = function(;return_value = true, header_name = true)
                describe_dict = Dict()
                if instance._data_type == "Numeric"
                    describe_dict["count"] = instance.count()
                    describe_dict["min"] = instance.min()
                    describe_dict["max"] = instance.max()
                    describe_dict["sum"] = instance.sum()
                    describe_dict["std"] = instance.std()
                    describe_dict["avg"] = instance.avg()
                else
                    describe_dict["unique"] = instance.unique()
                end

                describe_string = join(["$key : $value" for (key, value) in describe_dict], "\n")
                if header_name
                    describe_string = string(instance._col_name, "@", instance._parent_name, "\n", describe_string)
                else
                    describe_string = string("ColumnSet-", instance._col_num, "@", instance._parent_name, "\n", describe_string)
                end
                if return_value
                    println(describe_string)
                    return describe_dict
                else
                    return describe_string
                end
            end

            instance.count = function()
                return length(instance.data)
            end

            return instance
        end
    end
    function Base.show(io::Base.IO, instance::ColumnSet)
        println(io, string(instance.describe(return_value = false, header_name = false), "\n\n"))
    end

    mutable struct RowSet
        ###parameters
        #arguments
        data::DataStructures.OrderedDict
        #hidden parameters
        _parent_name::String
        _row_num::Int

        ###functions
        null::Function
        keys::Function

        #init
        function RowSet(data = DataStructures.OrderedDict(), row_num = 0, parent_name = "moeum.MOEUM")
            instance = new(data) ; instance._parent_name, instance._row_num = parent_name, row_num
            instance.null = function()
                result = []
                for (key, value) in instance.data
                    if value == nothing
                        Base.append!(result, key)
                    end
                end

                return result
            end

            instance.keys = function()
                result = Base.keys(instance.data)
                return result
            end

            return instance
        end
    end
    function Base.show(io::Base.IO, instance::RowSet)
        print_string = string("RowSet-", instance._row_num, "@", instance._parent_name, "\n")
        print_string = string(
            print_string,
            string("keys : ", instance.keys(), "\n")
        )

        println(io, print_string[1:end - 1])
    end
end


text:
### 4.1 - Packages

# Q1
#=
using Pkg
Pkg.add("Primes")
=#
# Alternatively, use the package manager in the REPL

# Q2
using Primes
#primes_lst = primes(1000000) # This takes a while.

### 4.2 - Plots

using Plots
#=
# Define data
globaltemperatures = [14.4, 14.5, 14.8, 15.2, 15.5, 15.8];
numpirates = [45000, 20000, 15000, 5000, 400, 17];

# Activate a plotting backend
gr()

# Plot data
p = plot(numpirates, globaltemperatures, label="line", size = (1200, 800));

scatter!(p, numpirates, globaltemperatures, label="points") ;
=#

# Q1

gr()

x = -10:0.2:10;
y = x.^2;

p = plot(x,y, legend = false);
scatter!(p, x, y)

# Q2

p1 = plot(x, x);
p2 = plot(x, x.^2);
p3 = plot(x, x.^3);
p4 = plot(x, x.^4);
plot(p1, p2, p3, p4, layout = (4, 1), legend = false)

max_stars_repo_path: Exercise Solutions/4 - Packages and Plotting.jl