text:
module System_parameters
    using Random
    export Params

    # Physical setting
    printlist_physical = ["L","β","NC","Nthermalization","Nsteps","initial","initialtrj","update_method","useOR","numOR","Nwing"]
    # Physical setting(fermions)
    printlist_fermions  = ["quench","Dirac_operator","Clover_coefficient","r","hop","Nf","mass",
                            "Domainwall_M","Domainwall_m","Domainwall_N5","Domainwall_ωs","Domainwall_b","Domainwall_c","Domainwall_r",
                            "BoundaryCondition",
                                "smearing_for_fermion","stout_numlayers","stout_ρ","stout_loops"]
    # System Control
    printlist_systemcontrol = ["log_dir","logfile","loadU_format","loadU_dir","loadU_fromfile","loadU_filename","saveU_dir","saveU_format",
                                "saveU_every","verboselevel","randomseed","measurement_basedir","measurement_dir","julian_random_number"]
    # HMC related
    printlist_HMCrelated = ["Δτ","SextonWeingargten","N_SextonWeingargten","MDsteps","eps","MaxCGstep"]
    # Action parameter for SLMC
    printlist_parameters_SLMC = ["βeff","firstlearn","use_autogeneratedstaples","couplingcoeff","couplinglist","coupling_loops"]
    # Measurement set
    printlist_measurement = ["measurement_methods"]

    const printlists = [printlist_physical,printlist_fermions,printlist_systemcontrol,printlist_HMCrelated,printlist_parameters_SLMC,printlist_measurement]
    const printlists_header = ["# Physical setting ","# Physical setting(fermions)","# System Control","# HMC related","# Action parameter for SLMC","# Measurement set"]

    defaultmeasures = Array{Dict,1}(undef,2)
    for i=1:length(defaultmeasures)
        defaultmeasures[i] = Dict()
    end
    defaultmeasures[1]["methodname"] = "Plaquette"
    defaultmeasures[1]["measure_every"] = 1
    defaultmeasures[1]["fermiontype"] = nothing
    defaultmeasures[2]["methodname"] = "Polyakov_loop"
    defaultmeasures[2]["measure_every"] = 1
    defaultmeasures[2]["fermiontype"] = nothing

    system = Dict()
    actions = Dict()
    md = Dict()
    cg = Dict()
    wilson = Dict()
    staggered = Dict()
    measurement = Dict()

    system["Dirac_operator"]  = "Wilson"
    system["quench"] = false
    system["initial"] = "cold"
    system["BoundaryCondition"]=[1,1,1,-1]
    system["Nwing"]= 1
    system["randomseed"] = 111
    #system["IntegratedHMC"] = false
    #system["Heatbath"] = false
    system["verboselevel"] = 2

    system["saveU_format"] = "JLD"
    system["saveU_every"] = 1
    system["saveU_dir"] = "./confs"

    system["initialtrj"] = 1

    #system["loadU_format"] = "JLD"
    #system["loadU_dir"] = "./confs"

    system["update_method"] = "HMC"

    actions["use_autogeneratedstaples"] = false
    #actions["couplinglist"] = []
    actions["couplingcoeff"] = []
    actions["coupling_loops"] = nothing

    md["Δτ"] = 0.05
    md["MDsteps"] = 20
    md["SextonWeingargten"] = false
    system["Nsteps"] = 100
    system["Nthermalization"] = 0

    #md["N_SextonWeingargten"]

    cg["eps"] = 1e-19
    cg["MaxCGstep"] = 3000

    #wilson["hop"]
    wilson["r"] = 1
    wilson["Clover_coefficient"] =  1.5612

    #staggered["mass"]
    staggered["Nf"] = 4

    measurement["measurement_methods"] = defaultmeasures

    system["julian_random_number"] = false

    mutable struct Params_set
        system::Dict
        actions::Dict
        md::Dict
        cg::Dict
        wilson::Dict
        staggered::Dict
        measurement::Dict

        function Params_set()
            return new(system,actions,md,cg,wilson,staggered,measurement)
        end

        function Params_set(system,actions,md,cg,wilson,staggered,measurement)
            return new(system,actions,md,cg,wilson,staggered,measurement)
        end

        function Params_set(system,actions,md,cg,measurement;wilson=Dict(),staggered=Dict())
            return new(system,actions,md,cg,wilson,staggered,measurement)
        end
    end

    function make_parameters(system,actions,md,cg,measurement;wilson=Dict(),staggered=Dict())
        return Params_set(system,actions,md,cg,measurement,wilson=wilson,staggered=staggered)
    end

    struct Params
        L::Tuple  # = (2,2,2,2) # Mandatory
        β::Float64 # = 6 # Mandatory
        NC::Int64 # = 3
        Dirac_operator::Union{Nothing,String} # = "Wilson"
        quench::Bool # = false
        initial::String # = "cold"
        BoundaryCondition::Array{Int8,1} # =[1,1,1,-1]
        Nwing::Int64  # = 1
        randomseed::Int64 # = 111

        #For actions
        use_autogeneratedstaples::Bool # = false
        couplinglist::Array{String,1} # = []
        couplingcoeff::Array{Float64,1} # = []
        coupling_loops::Union{Nothing,Array{Array{Array{Tuple{Int,Int},1},1},1}}

        #For MD
        Δτ::Float64 # = 0.05
        MDsteps::Int64 # = 20
        SextonWeingargten::Bool # = false
        Nsteps::Int64 # = 100
        Nthermalization::Int64 # = 10

        #For Sexton-Weingargten
        N_SextonWeingargten::Int64 # = 10 # Mandatory

        #For IntegratedHMC
        #IntegratedHMC::Bool # = false

        #For Heatbath
        #Heatbath::Bool # = false

        #For CG
        eps::Float64 # = 1e-19
        MaxCGstep::Int64 # = 3000

        #For Wilson

        hop::Float64 # = 0.141139#Hopping parameter # Mandatory
        r::Float64 # = 1 #Wilson term
        #For WilsonClover
        Clover_coefficient::Float64  # = 1.5612

        #For Staggered
        mass::Float64 # = 0.5 # Mandatory
        Nf::Int8 # = 4

        #verbose
        verboselevel::Int8 # = 2

        #For measurement
        measurement_methods::Array{Dict,1} # = defaultmeasures

        saveU_format::Union{String,Nothing}
        saveU_every::Int64
        saveU_dir::String

        loadU_format::Union{String,Nothing}
        loadU_dir::String

        update_method::String
        βeff::Union{Float64,Array{Float64,1}}
        firstlearn::Int64
        log_dir::String
        logfile::String
        load_fp::IOStream
        measuredir::String
        integratedFermionAction::Bool
        training_data_name::String
        useOR::Bool
        numOR::Int64
        initialtrj::Int64
        loadU_fromfile::Bool
        loadU_filename::String

        smearing_for_fermion::String
        stout_numlayers::Union{Nothing,Int64}
        stout_ρ::Union{Nothing,Array{Float64,1}}
        stout_loops::Union{Nothing,Array{String,1}}

        Domainwall_M::Union{Nothing,Float64}
        Domainwall_m::Union{Nothing,Float64}
        Domainwall_N5::Union{Nothing,Int64}
        Domainwall_b::Union{Nothing,Float64}
        Domainwall_c::Union{Nothing,Float64}
        Domainwall_ωs::Union{Nothing,Array{Float64,1}}
        Domainwall_r::Union{Nothing,Float64}

        julian_random_number::Bool

        function Params(system,actions,md,cg,wilson,staggered,measurement)
            #pnames = fieldnames(Params)
            #display(pnames)

            if haskey(system,"L")
                L = system["L"] #::Tuple  = (2,2,2,2) # Mandatory
            else
                error("system[\"L\"] should be set. eg. L = (4,4,4,4)")
            end

            if haskey(system,"β")
                β = system["β"] #::Float64 # = 6 # Mandatory
            else
                error("system[\"β\"] should be set. eg. β = 2.1")
            end

            NC = system["NC"] #::Int64 = 3
            Dirac_operator = system["Dirac_operator"] #::Union{Nothing,String} = "Wilson"
            quench = system["quench"] #::Bool = false
            initial = system["initial"] #::String = "cold"
            BoundaryCondition = system["BoundaryCondition"] #::Array{Int8,1}=[1,1,1,-1]
            Nwing = system["Nwing"] #::Int64 = 1
            randomseed = system["randomseed"] #::Int64 = 111

            use_autogeneratedstaples = actions["use_autogeneratedstaples"] #::Bool = false
            if use_autogeneratedstaples

                couplingcoeff = actions["couplingcoeff"] #::Array{Float64,1} = []
                #println(couplingcoeff)

                if haskey(actions,"coupling_loops")
                    coupling_loops = actions["coupling_loops"]
                else
                    error("actions[\"coupuling_loops\"] should be set if actions[\"use_autogeneratedstaples\"] is true")
                end
                if haskey(actions,"couplinglist")
                    @assert length(actions["couplinglist"]) == length(couplingcoeff) "the length of actions[\"couplinglist\"] should be equivalent to the length of actions[\"couplingcoeff\"]"
                    #println("1" )
                    couplinglist = actions["couplinglist"] #::Array{String,1} = []
                else
                    #println("2" )
                    couplinglist = String["$n" for n=1:length(couplingcoeff)]
                end
                #println(couplinglist )
            else
                couplinglist = []
                couplingcoeff = []
                coupling_loops = nothing
            end

            Δτ = md["Δτ"] #::Float64 = 0.05
            MDsteps = md["MDsteps"]#::Int64 = 20
            SextonWeingargten = md["SextonWeingargten"] #::Bool = false
            Nsteps = system["Nsteps"] #::Int64 = 100
            Nthermalization = system["Nthermalization"] #::Int64 = 10
            if SextonWeingargten
                if haskey(md,"N_SextonWeingargten")
                    N_SextonWeingargten = md["N_SextonWeingargten"]
                else
                    error("md[\"N_SextonWeingargten\"] should be set.")
                end
            else
                N_SextonWeingargten = 0
            end

            #For Heatbath
            #Heatbath = system["Heatbath"] #::Bool = false

            #For CG
            if quench == false
                eps = cg["eps"] #::Float64= 1e-19
                MaxCGstep = cg["MaxCGstep"] #::Int64 = 3000
            else
                eps = 1e-19
                MaxCGstep = 3000
            end

            #For Wilson
            if Dirac_operator == "Wilson" || Dirac_operator == "WilsonClover"
                println("$Dirac_operator fermion is used")
                if haskey(wilson,"hop")
                    hop = wilson["hop"] #::Float64 # = 0.141139#Hopping parameter # Mandatory
                else
                    error("wilson[\"hop\"] should be set.")
                end
                r = wilson["r"] #::Float64 = 1 #Wilson term

            else
                wilson["hop"]  = 0
                hop = 0
                r = 0
            end

            if Dirac_operator == "Domainwall"
                println("$Dirac_operator fermion is used")
                Domainwall_M = wilson["Domainwall_M"]
                Domainwall_m = wilson["Domainwall_m"]
                Domainwall_N5 = wilson["Domainwall_N5"]
                Domainwall_b = wilson["Domainwall_b"]
                Domainwall_c = wilson["Domainwall_c"]
                Domainwall_ωs =wilson["Domainwall_ωs"]
                Domainwall_r = wilson["Domainwall_r"]
                #r = wilson["r"]

            else
                Domainwall_N5 = nothing
                Domainwall_b = nothing
                Domainwall_c = nothing
                Domainwall_ωs = nothing
                Domainwall_r = nothing
                Domainwall_M = nothing
                Domainwall_m = nothing
            end

            #For WilsonClover
            if Dirac_operator == "WilsonClover"
                Clover_coefficient = wilson["Clover_coefficient"] #::Float64  = 1.5612
            else
                wilson["Clover_coefficient"] = 0
                Clover_coefficient = 0
            end

            #For Staggered
            if Dirac_operator == "Staggered"
                println("Staggered fermion is used")
                if haskey(staggered,"mass")
                    mass = staggered["mass"] #::Float64 # = 0.5 # Mandatory
                else
                    error("staggered[\"mass\"] should be set.")
                end
                Nf = staggered["Nf"] #::Int8 = 4
            else
                staggered["mass"] = 0
                mass = 0
                staggered["Nf"] = 0
                Nf = staggered["Nf"] #::Int8 = 4
            end

            verboselevel = system["verboselevel"]

            measurement_methods = measurement["measurement_methods"]

            saveU_format = system["saveU_format"]
            seveU_every = system["saveU_every"]
            saveU_dir = system["saveU_dir"]

            if saveU_format ≠ nothing
                if isdir(saveU_dir) == false
                    mkdir(saveU_dir)
                end
                println("$saveU_dir is used for saving configurations")
            end

            βeff = 0.0
            firstlearn = 0

            update_method = system["update_method"]
            useOR = false
            numOR = 0

            if update_method == "HMC"
                println("HMC will be used")
            elseif update_method == "Heatbath"
                println("Heatbath will be used")

                if system["quench"] == false
                    error("system[\"quench\"] = false. The Heatbath method needs the quench update. Put the other system[\"update_method\"] != \"Heatbath\" or system[\"quench\"] = true")
                end

                useOR = set_params(system,"useOR",false)
                numOR = set_params(system,"numOR",0)

            elseif update_method == "IntegratedHMC"
                println("IntegratedHMC will be used")

                #IntegratedHMC = system["IntegratedHMC"]
                if system["quench"] == false
                    println("system[\"quench\"] = true is set")
                    system["quench"] = true
                    #error("system[\"quench\"] = false. The IntegratedHMC needs the quench update. Put the other system[\"update_method\"] != \"IntegratedHMC\" or system[\"quench\"] = true")
                end

            elseif update_method == "Fileloading"
                println("No update will be used (read-measure mode)")
                system["quench"] = true
            elseif update_method == "SLHMC"
                println("SLHMC will be used")
                if system["quench"] == false
                    println("system[\"quench\"] = true is set")
                    system["quench"] = true
                    #error("system[\"quench\"] = false. The SLHMC needs the quench update. Put the other system[\"update_method\"] != \"SLHMC\" or system[\"quench\"] = true")
                end

                if haskey(system,"βeff")
                    βeff = system["βeff"]
                else
                    error("system[\"βeff\"] should be set when you want to do SLHMC.")
                end

                if haskey(system,"firstlearn")
                    firstlearn = system["firstlearn"]
                else
                    error("system[\"firstlearn\"] should be set when you want to do SLHMC.")
                end
            elseif update_method == "SLMC"
                println("SLMC will be used")
                if system["quench"] == false
                    println("system[\"quench\"] = true is set")
                    system["quench"] = true
                    #error("system[\"quench\"] = false. The SLMC needs the quench update. Put the other system[\"update_method\"] != \"SLMC\" or system[\"quench\"] = true")
                end

                if haskey(system,"βeff")
                    βeff = system["βeff"]
                else
                    error("system[\"βeff\"] should be set when you want to do SLMC.")
                end

                if haskey(system,"firstlearn")
                    firstlearn = system["firstlearn"]
                else
                    error("system[\"firstlearn\"] should be set when you want to do SLMC.")
                end
            elseif update_method == "IntegratedHB"
                println("IntegratedHB will be used")

                #IntegratedHMC = system["IntegratedHMC"]
                if system["quench"] == false
                    println("system[\"quench\"] = true is set")
                    system["quench"] = true
                    #error("system[\"quench\"] = false. The IntegratedHB needs the quench update. Put the other system[\"update_method\"] != \"IntegratedHB\" or system[\"quench\"] = true")
                end
            else
                error("""
                system[\"update_method\"] = $update_method is not supported.
                Supported methods are
                HMC
                Heatbath
                IntegratedHMC
                Fileloading
                """)
            end

            if update_method == "Fileloading"
                if haskey(system,"loadU_format")
                    loadU_format = system["loadU_format"]
                else
                    error("""system[\"loadU_format\"] should be set.
                    Supported formats are
                    JLD
                    """)
                end

                if haskey(system,"loadU_dir")
                    loadU_dir = system["loadU_dir"]
                else
                    error("system[\"loadU_dir\"] should be set")
                end

            else
                if haskey(system,"loadU_format")
                    loadU_format = system["loadU_format"]
                else
                    loadU_format = nothing
                end

                loadU_dir = "."
            end

            if haskey(system,"log_dir")
                log_dir = system["log_dir"]
                if isdir(log_dir) == false
                    mkdir(log_dir)
                end
            else
                error("system[\"log_dir\"] should be set")
            end

            if haskey(system,"logfile")
                logfile = pwd()*"/"*log_dir*"/"*system["logfile"]
            else
                error("system[\"logfile\"] should be set")
            end
            load_fp = open(logfile,"a")

            if haskey(measurement,"measurement_basedir")
                measurement_basedir = measurement["measurement_basedir"]
                if isdir(measurement_basedir) == false
                    mkdir(measurement_basedir)
                end
            else
                measurement_basedir = "nothing"
            end

            if haskey(measurement,"measurement_dir")
                measurement_dir = measurement["measurement_dir"]
                if isdir(pwd()*"/"*measurement_basedir*"/"*measurement_dir) == false
                    mkdir(pwd()*"/"*measurement_basedir*"/"*measurement_dir)
                end
            else
                measurement_dir = "nothing"
            end
            measuredir = pwd()*"/"*measurement_basedir*"/"*measurement_dir

            if haskey(actions,"IntegratedFermionAction")
                integratedFermionAction = actions["IntegratedFermionAction"]
            else
                integratedFermionAction = false
            end

            if haskey(actions,"training_data_name")
                training_data_name = pwd()*"/"*actions["training_data_name"]
            else
                training_data_name = pwd()*"/"*"trainingdata.txt"
            end

            if haskey(system,"initialtrj")
                initialtrj = system["initialtrj"]
            else
                initialtrj  = 1
            end

            if haskey(system,"loadU_fromfile")
                loadU_fromfile = system["loadU_fromfile"]
            else
                loadU_fromfile = false
            end

            #=
            if haskey(system,"loadU_filename")
                loadU_filename = system["loadU_filename"]
            else
                loadU_filename = ""
            end
            =#

            loadU_filename =set_params(system,"loadU_filename","")
            smearing_for_fermion = set_params(system,"smearing_for_fermion","nothing")
            stout_numlayers = set_params(system,"stout_numlayers",nothing)
            stout_ρ = set_params(system,"stout_ρ",nothing)
            stout_loops = set_params(system,"stout_loops",nothing)

            julian_random_number = set_params(system,"julian_random_number",true)

            #error("r = $r desu")
            return new(
                L,# = system["L"] #::Tuple  = (2,2,2,2) # Mandatory
                β,# = system["β"] #::Float64 # = 6 # Mandatory
                NC,# = system["NC"] #::Int64 = 3
                Dirac_operator,# = system["Dirac_operator"] #::Union{Nothing,String} = "Wilson"
                quench,# = system["quench"] #::Bool = false
                initial,# = system["initial"] #::String = "cold"
                BoundaryCondition,# = system["BoundaryCondition"] #::Array{Int8,1}=[1,1,1,-1]
                Nwing,# = system["Nwing"] #::Int64 = 1
                randomseed,# = system["randomseed"] #::Int64 = 111
                use_autogeneratedstaples,# = actions["use_autogeneratedstaples"] #::Bool = false
                couplinglist,# = actions["couplinglist"] #::Array{String,1} = []
                couplingcoeff,# = actions["couplingcoeff"] #::Array{Float64,1} = []
                coupling_loops,
                Δτ,# = md["Δτ"] #::Float64 = 0.05
                MDsteps,# = md["MDsteps"]#::Int64 = 20
                SextonWeingargten,# = md["SextonWeingargten"] #::Bool = false
                Nsteps,# = system["Nsteps"] #::Int64 = 100
                Nthermalization,# = system["Nthermalization"] #::Int64 = 10
                N_SextonWeingargten,# = md["N_SextonWeingargten"]
                #IntegratedHMC,# = system["IntegratedHMC"]
                #Heatbath,# = system["Heatbath"] #::Bool = false
                eps,# = cg["cg"] #::Float64= 1e-19
                MaxCGstep,# = cg["MaxCGstep"] #::Int64 = 3000
                hop,# = wilson["hop"] #::Float64 # = 0.141139#Hopping parameter # Mandatory
                r,# = wilson["r"] #::Float64 = 1 #Wilson term
                Clover_coefficient,# = wilson["Clover_coefficient"] #::Float64  = 1.5612
                mass,# = staggered["mass"] #::Float64 # = 0.5 # Mandatory
                Nf,# = staggered["Nf"] #::Int8 = 4
                verboselevel,
                measurement_methods,
                saveU_format,
                seveU_every,
                saveU_dir,
                loadU_format,
                loadU_dir,
                update_method,
                βeff,
                firstlearn,
                log_dir,
                logfile,
                load_fp,
                measuredir,
                integratedFermionAction,
                training_data_name,
                useOR,
                numOR,
                initialtrj,
                loadU_fromfile,
                loadU_filename,
                smearing_for_fermion,
                stout_numlayers,
                stout_ρ,
                stout_loops,
                Domainwall_M,#::Float64
                Domainwall_m,#::Float64
                Domainwall_N5,#::Int64
                Domainwall_b,#::Float64
                Domainwall_c,#::Float64
                Domainwall_ωs,#::Union{Nothing,Array{Float64,1}}
                Domainwall_r,
                julian_random_number
            )

        end

        function Params(params_set::Params_set)
            return Params(params_set.system,params_set.actions,params_set.md,params_set.cg,params_set.wilson,params_set.staggered,params_set.measurement)
        end

    end

    function set_params(dict,string,default)
        if haskey(dict,string)
            return dict[string]
        else
            return default
        end
    end

    function make_parametersdict(p::T) where T <: Any
        pnames = fieldnames(T)
        pdict = Dict()
        for i=1:length(pnames)
            pdict[String(pnames[i])] = getfield(p,pnames[i])
        end
        return pdict,pnames
    end

    function print_parameters_file(p)
        filename = p.logfile*"_parameters.jl"# "parameters_used.jl"
        fp = open(filename,"w")
        println(fp,"# - - parameters - - - - - - - - - - - ")
        pdict,pnames = make_parametersdict(p)
        for param in pdict
            if typeof(param[2]) == String
                println(fp,"$(param[1]) = \"$(param[2])\"")
            else
                println(fp,"$(param[1]) = $(param[2])")
            end
        end
        println(fp,"# - - - - - - - - - - - - - - - - - - -")
        close(fp)
        return
    end

    #=
    function gather_allparameters(params_set::Params_set)
        allparams = Dict()
        for (name,key) in params_set.system
            allparams[name] = key
        end

        for (name,key) in params_set.actions
            allparams[name] = key
        end

        for (name,key) in params_set.md
            allparams[name] = key
        end

        for (name,key) in params_set.cg
            allparams[name] = key
        end

        for (name,key) in params_set.wilson
            allparams[name] = key
        end

        for (name,key) in params_set.staggered
            allparams[name] = key
        end

        for (name,key) in params_set.measurement
            allparams[name] = key
        end
        return allparams
    end
    =#

    function setprint(fp,fp2,string)
        println(fp,string)
        if fp2 != nothing
            println(fp2,string)
        end
        println(string)
    end

    function get_stringfromkey(key)
        if typeof(key) == String
            string = "\"$key\""
        else
            string = "$key"
        end
    end

    function get_header(params_set::Params_set,inputname)
        if haskey(params_set.system,inputname)
            return "system",params_set.system[inputname]
        elseif haskey(params_set.actions,inputname)
            return "actions",params_set.actions[inputname]
        elseif haskey(params_set.md,inputname)
            return "md",params_set.md[inputname]
        elseif haskey(params_set.cg,inputname)
            return "cg",params_set.cg[inputname]
        elseif haskey(params_set.wilson,inputname)
            return "wilson",params_set.wilson[inputname]
        elseif haskey(params_set.staggered,inputname)
            return "staggered",params_set.staggered[inputname]
        elseif haskey(params_set.measurement,inputname)
            return "measurement",params_set.measurement[inputname]
        else
            return nothing,nothing
        end
    end

    function print_measurementinfo(fp,fp2,key)
        string = "measurement[\"measurement_methods\"] = Dict[ "
        setprint(fp,fp2,string)
        for (i,data) in enumerate(key)
            string = "  Dict{Any,Any}(\"methodname\" => \"$(data["methodname"])\","
            setprint(fp,fp2,string)
            count = 0
            #println(data)
            for (name,key_i) in data
                #println(name)
                #println(key_i)
                if name != "methodname"
                    count += 1
                    paramstring = get_stringfromkey(key_i)
                    string = "    \"$(name)\" => "*paramstring
                    if count != length(data) -1
                        string *= ","
                    end
                    if name == "fermiontype" && key_i == nothing
                    else
                        setprint(fp,fp2,string)
                    end

                end
            end
            string = "  )"
            if i != length(key)
                string *= ","
            end
            setprint(fp,fp2,string)

        end
        string = "]"
        setprint(fp,fp2,string)
    end

    function print_parameters_list(params_set::Params_set,p=nothing;filename=nothing)
        if filename == nothing
            @assert p != nothing "wrong input!"

            filename = p.logfile*"_parameters.jl"
            fp2 = p.load_fp
        else
            fp2 = nothing
        end
        fp = open(filename,"w")

        setprint(fp,fp2,"# - - parameters - - - - - - - - - - - ")
        for (i,printlist_i) in enumerate(printlists)
            setprint(fp,fp2,printlists_header[i] )
            for name in printlist_i
                headstring,key = get_header(params_set,name)
                if headstring != nothing
                    if name == "measurement_methods"
                        print_measurementinfo(fp,fp2,key)
                    else
                        string = get_stringfromkey(key)
                        paramstring = headstring*"[\"$name\"] = "*string
                        setprint(fp,fp2,paramstring)
                    end
                end
            end
            setprint(fp,fp2,"\t" )
        end

        setprint(fp,fp2,"# - - - - - - - - - - - - - - - - - - -")

        close(fp)

        println("""
        # Your parameters were written in $filename
        # If you want to do the simulation with same parameters,
        # Just do
        # julia run.jl $filename
        """)
        flush(stdout)

    end

    function print_parameters(params_set::Params_set,p)
        #allparams = gather_allparameters(params_set)
        print_parameters_list(params_set,p)
        return
    end

    function print_parameters(filename,params_set::Params_set)
        #allparams = gather_allparameters(params_set)
        print_parameters_list(params_set,filename=filename)
    end

    function print_parameters(p)
        println("# - - parameters - - - - - - - - - - - ")

        pdict,pnames = make_parametersdict(p)
        for param in pdict
            if typeof(param[2]) == String
                println("$(param[1]) = \"$(param[2])\"")
            else
                println("$(param[1]) = $(param[2])")
            end
        end
        println("# - - - - - - - - - - - - - - - - - - -")

        print_parameters_file(p)
        println("""
        # Your parameters were written in parameters_used.jl
        # If you want to do the simulation with same parameters,
        # Just do
        # julia run.jl parameters_used.jl
        """)
        flush(stdout)
        return
    end

    #=
    system  = System_params()
    md = MDparams()
    fermimon = Wilson_parames()
    measurement = Measurement_parames()
    =#

    #=
    if length(ARGS) > 0
        include(pwd()*"/"*ARGS[1])
    end
    =#

    function defaultparameters()
        return Params_set(system,actions,md,cg,wilson,staggered,measurement)
    end

    #display(system)
    #exit()

    function parameterloading(system,actions,md,cg,wilson,staggered,measurement)
        param_set = Params_set(system,actions,md,cg,wilson,staggered,measurement)
        if param_set.system["Dirac_operator"] == nothing
            param_set.system["quench"] = true
        end
        p = Params(param_set)
        Random.seed!(p.randomseed)

        print_parameters(param_set,p)
        return p
    end

    function parameterloading(param_set::Params_set)
        if param_set.system["Dirac_operator"] == nothing
            param_set.system["quench"] = true
        end
        p = Params(param_set)
        Random.seed!(p.randomseed)

        print_parameters(param_set,p)
        return p
    end

    function parameterloading()
        param_set = defaultparameters()

        if param_set.system["Dirac_operator"] == nothing
            param_set.system["quench"] = true
        end
        p = Params(param_set)
        Random.seed!(p.randomseed)

        print_parameters(param_set)
        return p

    end

end

#using .System_parameters
#System_parameters.parameterloading(ARGS[1])


text:
using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkitStandardLibrary.Blocks: Dₜ, t

#=
Testing strategy:
The general strategy is to test systems using simple intputs where the solution
is known on closed form. For algebraic systems (without differential variables),
an integrator with a constant input is often used together with the system under test.
=#

@testset "Constant" begin
    @info "Testing Constant"
    @named c = Constant(42)
    sys = structural_simplify(c)
    prob = ODEProblem(sys, Pair[], (0.0, 1.0))
    @test_broken begin # StackOverflowError:
        sol = solve(prob, Rosenbrock23())
        @test all(==(42), sol[c.y])
    end
    @named int = Integrator(; k=1)
    @named iosys = ODESystem([int.u~c.y], t, systems=[c, int])
    sys = structural_simplify(iosys)
    prob = ODEProblem(sys, Pair[], (0.0, 1.0))
    sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1)
    @test sol[int.y] ≈ 42 .* (0:0.1:1)
end

@testset "Integrator" begin
    @info "Testing Integrator"
    for k = [0.1, 1, 10]
        @named int = Integrator(; k)
        @test count(ModelingToolkit.isinput, states(int)) == 1
        @test count(ModelingToolkit.isoutput, states(int)) == 1
        @named iosys = ODESystem([int.u~1], t, systems=[int])
        sys = structural_simplify(iosys)
        stateind(sym) = findfirst(isequal(sym),states(sys))
        prob = ODEProblem(sys, Pair[int.u=>1.], (0.0, 1.0))
        sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1)
        @test sol[int.y] ≈ k .* (0:0.1:1)
    end
end

@testset "Derivative" begin
    @info "Testing Derivative"

    #= Derivative
    The test output below is generated by
    using ControlSystems
    sys = ss(-1/T, 1/T, -k/T, k/T)
    tv = 0:0.5:10
    u = (x,t)->[sin(t)]
    y = vec(lsim(sys, u, tv, alg=Rosenbrock23())[1])
    =#
    k = 1
    T = 0.1

    y01 = [0.0, 0.9096604391560481, 0.6179369162956885, 0.16723968919320775, -0.3239425882305049, -0.7344654437585882, -0.9662915429884467, -0.9619031643363591, -0.7219189996926385, -0.3046471954239838, 0.18896274787342904, 0.6325612488150467, 0.923147635361496, 0.9882186461533009, 0.8113758856575801, 0.4355269842556595, -0.05054266121798534, -0.5180957852231662, -0.8615644854197235, -0.994752654263345, -0.8845724777509947]

    y1 = [0.0, 0.37523930001382705, 0.5069379343173124, 0.422447016206449, 0.17842742193310424, -0.14287580928455357, -0.44972307981519677, -0.6589741190943343, -0.7145299845867902, -0.5997749247850142, -0.34070236779586216, -5.95731929625698e-5, 0.33950710748637825, 0.595360048429, 0.7051403889991136, 0.6421181090255983, 0.4214753349401378, 0.09771852881756515, -0.24995564964733635, -0.5364893060362096, -0.6917461951831227]

    y10 = [0.0, 0.04673868865158038, 0.07970450452536708, 0.09093906605247397, 0.07779607227750623, 0.04360203242101193, -0.0031749143460660587, -0.050989771426848074, -0.08804727520541561, -0.10519046453331109, -0.09814083278784949, -0.06855209962041757, -0.023592611490189652, 0.025798926487949535, 0.0675952553752348, 0.0916256775597053, 0.09206230764744555, 0.06885879535935949, 0.027748930190142837, -0.021151336671582116, -0.06582115823326284]

    for k = [0.1, 1, 10], (T,y) = zip([0.1, 1, 10], [y01, y1, y10])
        @named der = Derivative(; k, T)
        @test count(ModelingToolkit.isinput, states(der)) == 1
        @test count(ModelingToolkit.isoutput, states(der)) == 1
        @named iosys = ODESystem([der.u~sin(t)], t, systems=[der])
        sys = structural_simplify(iosys)
        stateind(sym) = findfirst(isequal(sym),states(sys))
        prob = ODEProblem(sys, Pair[der.u=>0., der.x=>0], (0.0, 10.0))
        sol = solve(prob, Rosenbrock23(), saveat=0:0.5:10)
        # plot([sol[der.y] k.*y]) |> display
        @test sol[der.y] ≈ k .* y rtol=1e-2
    end
end

@testset "FirstOrder" begin
    @info "Testing FirstOrder"

    for k = [0.1, 1, 10], T = [0.1, 1, 10]
        @named fo = FirstOrder(; k, T)
        @test count(ModelingToolkit.isinput, states(fo)) == 1
        @test count(ModelingToolkit.isoutput, states(fo)) == 1
        @named iosys = ODESystem([fo.u~1], t, systems=[fo])
        sys = structural_simplify(iosys)
        prob = ODEProblem(sys, Pair[fo.u=>1., fo.x=>0], (0.0, 10.0))
        sol = solve(prob, Rosenbrock23(), saveat=0:0.1:10)
        y = k .* (1 .- exp.(.-sol.t ./ T)) # Known solution to first-order system
        # plot([sol[fo.y] y]) |> display
        @test sol[fo.y] ≈ y rtol=1e-3
    end

end

@testset "SecondOrder" begin
    @info "Testing SecondOrder"

    # The impulse response of a second-order system with damping d follows the equations below
    function so(t,w,d)
        val = if d == 0
            1/w * sin(w*t)
        elseif d < 1
            1/(w*sqrt(1-d^2)) * exp(-d*w*t) * sin(w*sqrt(1-d^2)*t)
        elseif d == 1
            t*exp(-w*t)
        else
            1/(w*sqrt(d^2-1)) * exp(-d*w*t) * sinh(w*sqrt(d^2-1)*t)
        end
        val
    end

    w = 1
    d = 0.5
    for k = [0.1, 1, 10], w = [0.1, 1, 10], d = [0, 0.01, 0.1, 1, 1.1]
        @named sos = SecondOrder(; k, w, d)
        @test count(ModelingToolkit.isinput, states(sos)) == 1
        @test count(ModelingToolkit.isoutput, states(sos)) == 1
        @named iosys = ODESystem([sos.u~0], t, systems=[sos])
        sys = structural_simplify(iosys)
        prob = ODEProblem(sys, Pair[sos.u=>0.0, sos.xd=>1.0], (0.0, 10.0)) # set initial derivative state to 1 to simulate an impulse response
        sol = solve(prob, Rosenbrock23(), saveat=0:0.1:10, reltol=1e-6)
        y =  so.(sol.t,w,d)# Known solution to second-order system
        # plot([sol[sos.y] y]) |> display
        @test sum(abs2, sol[sos.y] - y) < 1e-4
    end
end

@testset "PID" begin
    @info "Testing PID"

    k = 2
    Ti = 0.5
    Td = 0.7
    wp = 1
    wd = 1
    Ni = √(Td / Ti)
    Nd = 12
    y_max = Inf
    y_min = -Inf
    u_r = sin(t)
    u_y = 0
    function solve_with_input(; u_r, u_y,
        controller = PID(; k, Ti, Td, wp, wd, Ni, Nd, y_max, y_min, name=:controller)
    )
        @test count(ModelingToolkit.isinput, states(controller)) == 5 # 2 in PID, 1 sat, 1 I, 1 D
        @test count(ModelingToolkit.isoutput, states(controller)) == 4
        # TODO: check number of unbound inputs when available, should be 2
        @named iosys = ODESystem([controller.u_r~u_r, controller.u_y~u_y], t, systems=[controller])
        sys = structural_simplify(iosys)
        prob = ODEProblem(sys, Pair[], (0.0, 10.0))
        sol = solve(prob, Rosenbrock23(), saveat=0:0.1:10)
        controller, sys, sol
    end

    # linearity in u_r
    controller, sys, sol1 = solve_with_input(u_r=sin(t), u_y=0)
    controller, sys, sol2 = solve_with_input(u_r=2sin(t), u_y=0)
    @test sum(abs, sol1[controller.ea]) < eps() # This is the acutator model error due to saturation
    @test 2sol1[controller.y] ≈ sol2[controller.y] rtol=1e-3 # linearity in u_r

    # linearity in u_y
    controller, sys, sol1 = solve_with_input(u_y=sin(t), u_r=0)
    controller, sys, sol2 = solve_with_input(u_y=2sin(t), u_r=0)
    @test sum(abs, sol1[controller.ea]) < eps() # This is the acutator model error due to saturation
    @test 2sol1[controller.y] ≈ sol2[controller.y] rtol=1e-3 # linearity in u_y

    # zero error
    controller, sys, sol1 = solve_with_input(u_y=sin(t), u_r=sin(t))
    @test sum(abs, sol1[controller.y]) ≈ 0 atol=sqrt(eps())

    # test saturation
    controller, sys, sol1 = solve_with_input(; u_r=10sin(t), u_y=0,
        controller = PID(; k, Ti, Td, wp, wd=0, Ni, Nd, y_max=10, y_min=-10, name=:controller)
    )
    @test extrema(sol1[controller.y]) == (-10, 10)

    # test P set-point weighting
    controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0,
        controller = PID(; k, Ti, Td, wp=0, wd, Ni, Nd, y_max, y_min, name=:controller)
    )
    @test sum(abs, sol1[controller.ep]) ≈ 0 atol=sqrt(eps())

    # test D set-point weighting
    controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0,
        controller = PID(; k, Ti, Td, wp, wd=0, Ni, Nd, y_max, y_min, name=:controller)
    )
    @test sum(abs, sol1[controller.ed]) ≈ 0 atol=sqrt(eps())

    # zero integral gain
    controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0,
        controller = PID(; k, Ti=false, Td, wp, wd, Ni, Nd, y_max, y_min, name=:controller)
    )
    @test isapprox(sum(abs, sol1[controller.I.y]), 0, atol=sqrt(eps()))

    # zero derivative gain
    @test_skip begin # During the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number: [5]
        controller, sys, sol1 = solve_with_input(; u_r=sin(t), u_y=0,
            controller = PID(; k, Ti, Td=false, wp, wd, Ni, Nd, y_max, y_min, name=:controller)
        )
        @test isapprox(sum(abs, sol1[controller.D.y]), 0, atol=sqrt(eps()))
    end

    # Tests below can be activated when the concept of unbound_inputs exists in MTK
    # @test isequal(Set(unbound_inputs(controller)), @nonamespace(Set([controller.u_r, controller.u_y])))
    # @test isempty(unbound_inputs(sys))
    # @test isequal(bound_inputs(sys), inputs(sys))
    # @test isequal(
    #     Set(bound_inputs(sys)),
    #     Set([controller.u_r, controller.u_y, controller.I.u, controller.D.u, controller.sat.u])
    #     )
end

## Additional test of PID controller using ControlSystems
# using ControlSystems
# kd = 1
# Nd = 12
# Td = 1
# T = Td/Nd
# Cd = ss(-1/T, 1/T, -kd/T, kd/T) |> tf

# C = ControlSystems.pid(kp=10, ki=1, kd=0, series=true, time=true) + 10*Cd
# P = tf(1,[1, 0])^2
# L = ss(P*C)

# @named controller = PID(k=10, Ti=1, Td=1)
# @named plant = Blocks.StateSpace(ssdata(ss(P))...)
# @named iosys = ODESystem([
#     controller.u_r~1,
#     controller.u_y~plant.y[1],
#     controller.y~plant.u[1]
# ], t, systems=[controller, plant])
# sys = structural_simplify(iosys)
# prob = ODEProblem(sys, Pair[], (0.0, 6))
# sol = solve(prob, Rosenbrock23())

# res = step(feedback(L), sol.t)
# y = res.y[:]
# plot(res)
# plot!(sol, vars=[plant.y[1]])
# @test sol[plant.y[1]] ≈ y rtol = 1e-3
##

@testset "StateSpace" begin
    @info "Testing StateSpace"

    A = [0 1; 0 0]
    B = [0, 1]
    C = [1 0]
    D = 0
    @named sys = Blocks.StateSpace(A,B,C,D)
    @test count(ModelingToolkit.isinput, states(sys)) == 1
    @test count(ModelingToolkit.isoutput, states(sys)) == 1
    @named iosys = ODESystem([sys.u[1] ~ 1], t, systems=[sys])
    iosys = structural_simplify(iosys)
    prob = ODEProblem(iosys, Pair[], (0.0, 1.0))
    sol = solve(prob, Rosenbrock23(), saveat=0:0.1:1)
    @test sol[sys.x[2]] ≈ (0:0.1:1)
    @test sol[sys.x[1]] ≈ sol[sys.y[1]]

    D = randn(2, 2) # If there's only a `D` matrix, the result is a matrix gain
    @named sys = Blocks.StateSpace([],[],[],D)
    gain = Blocks.Gain(D, name=:sys)
    @test sys == gain
end

ext: jl


text:
# Note: Test environment is setup by https://github.com/JuliaCI/PkgEvalHadoopEnv
include("hdfstests.jl")
include("yarntests.jl")

test_ugi()
test_hdfs()
test_credential_tokens()
test_container_id()
test_yarn_client()
test_unmanaged_yarn_clustermanager()
test_managed_yarn_clustermanager()
@info("runtests.jl done")


text:
## Module associated with the MeshArrays package

module MeshArrays

using Printf

include("gcmfaces_type.jl");
include("gcmfaces_grids.jl");
include("gcmfaces_calc.jl");
include("gcmfaces_exch.jl");
include("gcmfaces_convert.jl");
include("gcmfaces_IO.jl");
include("gcmfaces_demo.jl");

export gcmfaces, exchange, gradient, convergence, smooth, mask, fsize
export GCMGridSpec, GCMGridLoad, GCMGridOnes
export demo1, demo2
#The following functions rely on grid specs; currently via global vars.
export read_bin, convert2array, convert2gcmfaces
#The following exch_UV differs from normal exchange; incl. exch_UV_N.
export exch_UV
#The following codes add dependencies to Plots & NetCDF.
#include("gcmfaces_plot.jl");
#include("gcmfaces_nctiles.jl");

end # module


text:
##### Moisture tendencies

export CreateClouds

#####
##### First order fluxes
#####

function flux(::TotalMoisture, ::Advect, atmos, args)
    @unpack state = args
    u = state.ρu / state.ρ
    return u * state.moisture.ρq_tot
end

function flux(::LiquidMoisture, ::Advect, atmos, args)
    @unpack state = args
    u = state.ρu / state.ρ
    return u * state.moisture.ρq_liq
end

function flux(::IceMoisture, ::Advect, atmos, args)
    @unpack state = args
    u = state.ρu / state.ρ
    return u * state.moisture.ρq_ice
end

#####
##### Second order fluxes
#####

function flux(::TotalMoisture, ::MoistureDiffusion, atmos, args)
    @unpack state, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    return d_q_tot * state.ρ
end

function flux(::LiquidMoisture, ::MoistureDiffusion, atmos, args)
    @unpack state, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    d_q_liq = (-D_t) .* diffusive.moisture.∇q_liq
    return d_q_liq * state.ρ
end

function flux(::IceMoisture, ::MoistureDiffusion, atmos, args)
    @unpack state, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    d_q_ice = (-D_t) .* diffusive.moisture.∇q_ice
    return d_q_ice * state.ρ
end

function flux(::TotalMoisture, ::HyperdiffViscousFlux, atmos, args)
    @unpack state, hyperdiffusive = args
    return hyperdiffusive.hyperdiffusion.ν∇³q_tot * state.ρ
end

#####
##### Sources
#####

function source(::TotalMoisture, s::Subsidence, m, args)
    @unpack state, aux, diffusive = args
    z = altitude(m, aux)
    w_sub = subsidence_velocity(s, z)
    k̂ = vertical_unit_vector(m, aux)
    return -state.ρ * w_sub * dot(k̂, diffusive.moisture.∇q_tot)
end

"""
    CreateClouds <: TendencyDef{Source}

A source/sink to `q_liq` and `q_ice` implemented as a relaxation towards
equilibrium in the Microphysics module.
The default relaxation timescales are defined in CLIMAParameters.jl.
"""
struct CreateClouds <: TendencyDef{Source} end

prognostic_vars(::CreateClouds) = (LiquidMoisture(), IceMoisture())

function source(::LiquidMoisture, s::CreateClouds, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    # get current temperature and phase partition
    FT = eltype(state)
    q = PhasePartition(ts)
    T = air_temperature(ts)
    param_set = parameter_set(m)

    # phase partition corresponding to the current T and q.tot
    # (this is not the same as phase partition from saturation adjustment)
    ts_eq = PhaseEquil_ρTq(param_set, state.ρ, T, q.tot)
    q_eq = PhasePartition(ts_eq)

    # cloud condensate as relaxation source terms
    S_q_liq = conv_q_vap_to_q_liq_ice(param_set.microphys.liq, q_eq, q)

    return state.ρ * S_q_liq
end

function source(::IceMoisture, s::CreateClouds, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    # get current temperature and phase partition
    FT = eltype(state)
    q = PhasePartition(ts)
    T = air_temperature(ts)
    param_set = parameter_set(m)

    # phase partition corresponding to the current T and q.tot
    # (this is not the same as phase partition from saturation adjustment)
    ts_eq = PhaseEquil_ρTq(param_set, state.ρ, T, q.tot)
    q_eq = PhasePartition(ts_eq)

    # cloud condensate as relaxation source terms
    S_q_ice = conv_q_vap_to_q_liq_ice(param_set.microphys.ice, q_eq, q)

    return state.ρ * S_q_ice
end

function source(::TotalMoisture, s::RemovePrecipitation, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    if has_condensate(ts)
        nt = remove_precipitation_sources(s, m, args)
        return nt.S_ρ_qt
    else
        FT = eltype(state)
        return FT(0)
    end
end

function source(::TotalMoisture, s::WarmRain_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_qt
end

function source(::LiquidMoisture, s::WarmRain_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_ql
end

function source(::TotalMoisture, s::RainSnow_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_qt
end

function source(::LiquidMoisture, s::RainSnow_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_ql
end

function source(::IceMoisture, s::RainSnow_1M, m, args)
    @unpack cache = args.precomputed.precipitation
    return cache.S_ρ_qi
end


text:
using FdeSolver
using Plots
using SpecialFunctions

## inputs
tSpan = [0, 1]     # [intial time, final time]
y0 = 0             # intial value
β = 0.5            # order of the derivative
h=2^-5

# Equation
par = β
F(t, y, par) = (40320 ./ gamma(9 - par) .* t .^ (8 - par) .- 3 .* gamma(5 + par / 2)
           ./ gamma(5 - par / 2) .* t .^ (4 - par / 2) .+ 9/4 * gamma(par + 1) .+
           (3 / 2 .* t .^ (par / 2) .- t .^ 4) .^ 3 .- y .^ (3 / 2))
# Jacobian
JacobF(t, y, par) = -(3 / 2) .* y .^ (1 / 2)

## Numerical solution
t, Yapp = FDEsolver(F, tSpan, y0, β, par, h = h, StopIt = "Convergence", tol = 10e-8, itmax = 30)

t1, Yapp1 = FDEsolver(F, tSpan, y0, β, par, JF = JacobF, tol = 10e-8, itmax = 30)

#plot
plot(t, Yapp, linewidth = 5, title = "Solution of a 1D fractional IVP",
     xaxis = "Time (t)", yaxis = "y(t)", label = "Approximation")
plot!(t1, Yapp1, linewidth = 5, ls = :dot, label = "Approximation with Jacobian")
plot!(t, t -> (t.^8 - 3 * t .^ (4 + β / 2) + 9 / 4 * t.^β),
      lw = 3, ls = :dash, label = "Exact solution")


text:
using FusionSystems
using Documenter

makedocs(;
    modules=[FusionSystems],
    authors="Dan Segal <dansegal2@gmail.com> and contributors",
    repo="https://github.com/djsegal/FusionSystems.jl/blob/{commit}{path}#L{line}",
    sitename="FusionSystems.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://djsegal.github.io/FusionSystems.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/djsegal/FusionSystems.jl",
)


text:
device_gc() = af_device_gc()

function device_mem_info()
  a_bytes = Csize_t[0]
  a_buffers = Csize_t[0]
  l_bytes = Csize_t[0]
  l_buffers = Csize_t[0]
  af_device_mem_info(a_bytes, a_buffers, l_bytes, l_buffers)
  Int(a_bytes[1]), Int(a_buffers[1]), Int(l_bytes[1]), Int(l_buffers[1])
end

function device_ptr(a::AFArray)
  p = Ptr{Void}[0]
  af_get_device_ptr(p, a)
  #finalizer(a, unlock_device_ptr)
  p[1]
end

function mem_stepsize()
  p = Csize_t[0]
  af_get_mem_step_size(p)
  Int(p[1])
end

lock_device_ptr(a::AFArray) = af_lock_device_ptr(a)

unlock_device_ptr(a::AFArray) = af_unlock_device_ptr(a)

set_mem_stepsize(bytes::Int) = af_set_mem_step_size(bytes)

free_device(a::AFArray) = af_free_device(a)

function device_count()
  p = Cint[0]
  af_get_device_count(p)
  Int(p[1])
end

function device_info()
  name = UInt8[0]
  platform = UInt8[0]
  toolkit = UInt8[0]
  compute = UInt8[0]
  af_device_info(name, platform, toolkit, compute)
  bytestring(name), bytestring(platform)
  #error("Not implemented yet.")
end

function getdevice(device::Int)
  p = Cint[0]
  af_get_device(p)
  Int(p[1])
end

setdevice(device::Int) = af_set_device(device)

sync(device::Int=-1) = af_sync(device)

function available_backends()
  p = Cint[0]
  af_get_available_backends(p)
  v = Int(p[1])
  cpu = (v & AF_BACKEND_CPU) == 0 ? "off" : "on"
  cuda = (v & AF_BACKEND_CUDA) == 0 ? "off" : "on"
  opencl = (v & AF_BACKEND_OPENCL) == 0 ? "off" : "on"
  "cpu:$(cpu), cuda:$(cuda), opencl:$(opencl)"
end

function set_backend(str::ASCIIString)
  if str == "cpu"
    b = AF_BACKEND_CPU
  elseif str == "cuda"
    b = AF_BACKEND_CUDA
  elseif str == "opencl"
    b = AF_BACKEND_OPENCL
  else
    throw("No backend: $(str)")
  end
  af_set_backend(b)
end


text:
using DiffEqBase: DiscreteCallback, terminate!
using NQCBase: apply_cell_boundaries!, check_atoms_in_cell

outside_cell(u,t,integrator) = !check_atoms_in_cell(integrator.p.cell, get_positions(u))
function enforce_periodicity!(integrator)
    apply_cell_boundaries!(integrator.p.cell, get_positions(integrator.u))
end

"""
    CellBoundaryCallback()

Whenever atoms leave the simulation cell, enforce the periodicity by wrapping the positions
at the cell boundary.
"""
CellBoundaryCallback() = DiscreteCallback(outside_cell, enforce_periodicity!)

"""
    TerminateCellCallback()

If the atoms leave the simulation cell, terminate the simulation.
"""
TerminateCellCallback() = DiscreteCallback(outside_cell, terminate!)

"""
    TerminatingCallback(func)

Provide a function that returns true when the simulation should terminate.
"""
TerminatingCallback(func) = DiscreteCallback(func, terminate!)


text:
import Base: copy
import Random: AbstractRNG, rand, SamplerType

export Family,
    fitness,
    genome,
    Individual,
    ongenome,
    Population

mutable struct Individual{G}
    genome::G
    fitnessvalue::Union{FitnessValue, Nothing}
end

Individual(genome::G) where {G} = Individual{G}(genome, nothing)
Individual(genome::G, fitnessvalue) where G =
    Individual{G}(genome, convert(FitnessValue, fitnessvalue))

genome(i::Individual) = i.genome
fitness(i::Individual) = i.fitnessvalue

copy(individual::Individual) = Individual(copy(individual.genome), individual.fitnessvalue)

rand(rng::AbstractRNG, ::SamplerType{Individual{T}}) where {T} = Individual{T}(rand(T), nothing)

"""
    ongenome(f, i[, is])

Apply `f` on the `genome` fields of all individuals: `f(i.genome, is[1].genome, ..., is[N].genome)`.
"""
function ongenome(f, i::Individual)
    f(i.genome)
end

function ongenome(f, i::Individual, is::Vararg{Individual, N}) where {N}
    if @generated
        :(f(i.genome, $((:(is[$n].genome) for n in 1:N)...)))
    else
        f(getfield.((i, is...), :genome)...)
    end
end

const Population{T} = Vector{Individual{T}}
(::Type{Population})(genomes) = Individual.(collect(genomes))

const Family{T, N} = NTuple{N, Individual{T}}
(::Type{Family{T, N}})(args::Vararg{T, N}) where {T,N} = args


text:
module TestEuclid
import ProjectiveGeometricAlgebra3d as PGA
using .PGA.Euclid
using .PGA.Euclid: pga, norm, join, random
using .PGA: ∨
using Test

E0 = Point(0,0,0)
E1 = Point(1,0,0)
E2 = Point(0,1,0)
E3 = Point(0,0,1)

@testset "distance" begin
    @test distance(E0, E1) ≈ 1
    @test distance(E0, E2) ≈ 1
    @test distance(E0, E3) ≈ 1

    @test distance(E1, E0) ≈ 1
    @test distance(E1, E2) ≈ sqrt(2)
    @test distance(E1, E3) ≈ sqrt(2)

    @test distance(E2, E0) ≈ 1
    @test distance(E2, E1) ≈ sqrt(2)
    @test distance(E2, E3) ≈ sqrt(2)

    @test distance(E3, E0) ≈ 1
    @test distance(E3, E1) ≈ sqrt(2)
    @test distance(E3, E2) ≈ sqrt(2)
end

@testset "norm" begin
    @test norm(pga(E0) ∨ pga(E1)) ≈ 1
    @test norm(pga(E1) ∨ pga(E2)) ≈ sqrt(2)
    c1 = randn(3)
    c2 = randn(3)
    @test norm(c1 - c2) ≈ norm(pga(Point(c1)) ∨ pga(Point(c2)))
end

@testset "angle between planes" begin
    p_yz = Plane(point=randn(3), directions=([0,0,rand()], (0,rand(),0)))
    p_2  = Plane(point=randn(3), directions=([0,0,rand()], rand()*[1,1,0]))
    @test Euclid.cos_angle(p_yz, p_2) ≈ cos(pi/4)

end

@testset "join and meet" begin
    @test join(E0,E1) ≈ Line(point=[0,0,0], direction=[1,0,0])
    @test join(E0,E1,E2) ≈ Plane(point=[0,0,0], directions=([1,0,0], [0,1,0]))

    pxy = join(E0,E1,E2)
    pyz = join(E0,E2,E3)
    @test meet(pxy, pyz) ≈ join(E0, E2)
end

@testset "project" begin
    x,y,z = randn(3)
    pt = Point(x,y,z)
    qt = project(join(E0,E1,E2), pt)
    @test qt ≈ Point([x,y,0])
    l1 = join(Point(0,0,0), Point(1,2,3))
    l2 = project(join(E0,E1,E2), l1)
    l3 = join(Point(0,0,0), Point(1,2,0))
    @test l2 ≈ l3
end

@testset "Motor" begin
    m = rotator(axis=join(E0, E1), angle=pi/2)
    @test m(Point(0,0,1)) ≈ Point(0,-1,0)

    for _ in 1:100
        m1 = random(Motor{Float64})
        pt1, pt2 = random(Point{Float64}, 2)
        m1(join(pt1, pt2)) ≈ join(m1(pt1), m1(pt2))
    end
    for Geo in [Point{Float32}, Line{Float32}, Plane{Float32}, Motor{Float32}]
        m1, m2 = random(Motor{Float32}, 2)
        geo = random(Geo)
        @test m1(m2(geo)) ≈ (m1∘m2)(geo)
    end
end

end#module


text:
using TopOpt
import Makie
using TopOpt.TopOptProblems.Visualization: visualize

using FromFile
@from "new_problems.jl" import NewPointLoadCantilever

using TimerOutputs

# function run_topopt()
    println("Start running.")
    # https://github.com/KristofferC/TimerOutputs.jl
    to = TimerOutput()
    reset_timer!(to)
    Nonconvex.show_residuals[] = true

    # Define the problem
    E = 1.0 # Young’s modulus
    v = 0.3 # Poisson’s ratio
    f = 1.0; # downward force

    # Parameter settings
    V = 0.3 # volume fraction
    # xmin = 0.001 # minimum density
    xmin = 1e-6 # minimum density
    rmin = 2.0; # density filter radius

    nels = (20, 10, 2)
    sizes = (1.0, 1.0, 1.0)
    @timeit to "problem def" problem = NewPointLoadCantilever(Val{:Linear}, nels, sizes, E, v, f);

    # Define a finite element solver
    @timeit to "penalty def" penalty = TopOpt.PowerPenalty(3.0)
    @timeit to "solver def" solver = FEASolver(Displacement, Direct, problem, xmin = xmin,
        penalty = penalty);

    # Define compliance objective
    @timeit to "objective def" begin
        # Define compliance objective
        comp = Compliance(problem, solver)
        filter = DensityFilter(solver, rmin = rmin)
        obj = x -> comp(filter(x))
    end

    # Define volume constraint
    @timeit to "constraint def" begin
        volfrac = TopOpt.Volume(problem, solver)
        constr = x -> volfrac(filter(x)) - V
    end

    # Define subproblem optimizer
    # TODO MMA02 converge to weird results under this criteria
    mma_options = options = Nonconvex.MMAOptions(
        maxiter = 400,
        tol = Nonconvex.Tolerance(x = 1e-3, f = 1e-6),
        )
    ipopt_options = Nonconvex.IpoptOptions(
        # maxiter = 100,
        # tol = Nonconvex.Tolerance(x = 1e-3, f = 1e-6),
    )

    convcriteria = Nonconvex.GenericCriteria()
    # convcriteria = Nonconvex.KKTCriteria()

    x0 = fill(V, length(solver.vars))
    @timeit to "optimizer def" optimizer = Optimizer(obj, constr, x0,
        Nonconvex.MMA87(),
        options = mma_options,
        # Nonconvex.IpoptAlg(),
        # options = ipopt_options,
        convcriteria = convcriteria
        );
    # Nonconvex.MMA02(),

    # Define SIMP optimizer
    @timeit to "simp def" simp = SIMP(optimizer, solver, penalty.p);

    # Solve
    # initial solution, critical to set it to volfrac! (blame non-convexity :)
    @timeit to "simp run" result = simp(x0);

    # Print the timings in the default way
    show(to)
    println()

    @show result.convstate
    @show result.objval
    try
        @show optimizer.workspace.iter
    catch
        # IpoptWorkspace has no field iter
    end

    # # Visualize the result using Makie.jl
    fig = visualize(problem; topology=result.topology,
        default_exagg_scale=0.07, scale_range=10.0, vector_linewidth=3, vector_arrowsize=0.5)
    Makie.display(fig)

    # return problem, result
# end

# problem, result = run_topopt();


text:
using EquivalenceClassesConstructor
using Test

using Printf
using DataStructures
using JLD
using Lazy
include("../src/Types.jl")
include("../src/GraphTools.jl")
include("../src/helpers.jl")
include("../src/IO.jl")

pred1 = Predicate((x,y) -> x == y)
pred2 = Predicate((x,y) -> x != y)
pred3 = Predicate((x,y) -> all(x .== -y))
pred4 = Predicate((x,y) -> x == y+3)

m1 = Mapping(x -> [x])
m2 = Mapping(x -> [-x])
vl1 = collect(1:5)
vl2 = collect(-3:5)
cl1 = find_classes(m1, vl1)
cl2 = find_classes(m2, vl2)
cl2_res_index  = Dict(3 => UInt32(1), -3 => UInt32(1), 2 => UInt32(2), -2=> UInt32(2), 1 => UInt32(3), -1 => UInt32(3), 0 => UInt32(4), 4 => UInt32(8), 5 => UInt32(9))
cl2_res_direct = Dict(3 => -3, -3 => -3, 2 => -2, -2=>-2, 1 => -1, -1 => -1, 0 => 0, 4 => 4, 5 => 5)

@testset "EquivalenceClassesConstructor.jl" begin
    @testset "helpers" begin
        include("./helpers.jl")
    end
    @testset "Types" begin
        include("./Types.jl")
    end
    @testset "graph" begin
        include("./GraphTools.jl")
    end
    @testset "IO" begin
        include("./IO.jl")
    end
end


text:
module ClausenFunctions

include("range_reduction.jl")
include("Cl1.jl")
include("Cl2.jl")
include("Cl3.jl")
include("Cl4.jl")
include("Cl5.jl")
include("Cl6.jl")
include("Cl.jl")
include("Sl.jl")

export cl1
export cl2
export cl3
export cl4
export cl5
export cl6
export cl
export sl

end


text:
using Base.Test
import COBRAUtils

@testset "E. Coli" begin
    model = COBRAUtils.loadEColiTestModel()

    @test COBRAUtils.metabolite_index(model, "6pgc[c]") == 4
    @test COBRAUtils.reaction_index(model, "ACONTa") == 4
    @test sort(model.mets[COBRAUtils.reaction_metabolites(model, "ACALD")]) == sort(["acald[c]", "accoa[c]", "coa[c]", "h[c]", "nad[c]", "nadh[c]"])
    @test sort(model.mets[COBRAUtils.reaction_substrates(model, "ACALD")]) == sort(["acald[c]", "coa[c]", "nad[c]"])
    @test sort(model.mets[COBRAUtils.reaction_products(model, "ACALD")]) == sort(["accoa[c]", "h[c]", "nadh[c]"])
    @test sort(model.rxns[COBRAUtils.metabolite_reactions(model, "ac[c]")]) == sort(["ACt2r", "ACKr"])

    @test COBRAUtils.reaction_substrates(model, "EX_ac(e)") == [7]
    @test isempty(COBRAUtils.reaction_products(model, "EX_ac(e)"))

    @test COBRAUtils.exchange_reactions(model, "ac[e]") == [20]

    status, μ, v = COBRAUtils.optimize(model)
    @test μ ≈ 0.8739215069684305
end

@testset "reduce_model E. Coli" begin
    ecoli = COBRAUtils.loadEColiTestModel();
    ecolired = COBRAUtils.reduce_model(ecoli);
    @test size(ecolired.S) == (68, 87)
    status0, μ0, v0 = COBRAUtils.optimize(ecoli)
    status1, μ1, v1 = COBRAUtils.optimize(ecolired)
    @test μ0 ≈ μ1
end


text:
#=
# testing
# Julia version:
Created by admin-bellei at 2018-02-18
=#
println("test")

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
@testset "Sphere" begin
    # Simple constructor calls
    @test Sphere(zeros(4)) ≈ 0.0 # This is the minimum
    @test Sphere(ones(4)) ≈ 4.0
    @test Sphere(0.0) ≈ 0.0
    @test Sphere(5.0) ≈ 25.0
    # Using multiple dispatch
    @test isa(evaluate(Sphere(), 0.0), Number)
end

@testset "Easom" begin
    @test_throws AssertionError Easom(zeros(4))
    @test Easom(fill(3.0, 2)) ≈ -0.94156415
    minimum_val = @SVector [π, π]
    @test Easom(minimum_val) ≈ -1.0 # Check for minimum
    @test isa(evaluate(Easom(), [2.0, 4.0]), Number)
end

@testset "Ackley" begin
    @test_throws AssertionError Ackley([])
    @test isapprox(Ackley(zeros(5)), 0.0, atol = 1e-15) # Check for minimum
    @test isa(evaluate(Ackley(), [2.0, 4.0]), Number)
end

@testset "Rosenbrock" begin
    @test_throws AssertionError Rosenbrock(zeros(1))
    @test isapprox(Rosenbrock(fill(1.0, 4)), 0.0, atol = 1e-15) # Check for minimum
end

@testset "GoldsteinPrice" begin
    ground_truth = 3.0
    true_vector = @SVector [0.0, -1.0]
    @test_throws AssertionError GoldsteinPrice(zeros(3))
    # Check for minimum
    @test isapprox(GoldsteinPrice(true_vector), ground_truth, atol = 1e-15)
end

@testset "Beale" begin
    ground_truth = 0
    true_vector = @SVector [3, 0.5]
    @test_throws AssertionError Beale(zeros(3))
    # Check for minimum
    @test isapprox(Beale(true_vector), ground_truth, atol = 1e-15)
end

@testset "Levy" begin
    ground_truth = 0
    dimension = 10
    true_vector = fill(1.0, dimension)
    # Check for minimum
    @test isapprox(Levy(true_vector), ground_truth, atol = 1e-15)
end


text:
module GXBeam

using ArnoldiMethod
using FLOWMath
using LinearAlgebra
using LinearMaps
using NLsolve
using SparseArrays
using StaticArrays
using WriteVTK

export curve_length, discretize_beam
export System, Assembly, PrescribedConditions, DistributedLoads
export reset_state!
export static_analysis, steady_state_analysis, eigenvalue_analysis, time_domain_analysis
export static_analysis!, steady_state_analysis!, eigenvalue_analysis!, time_domain_analysis!
export AssemblyState, left_eigenvectors, wiener_milenkovic, write_vtk, correlate_eigenmodes

# Constant used for scaling forces/moments
# this is needed because sparse arrays don't seem to be able to handle
# mismatched units very well
const FORCE_SCALING = nextpow(2, 1e3)

const GAUSS_NODES = SVector(-0.8611363115940526, -0.3399810435848563, 0.3399810435848563, 0.8611363115940526)
const GAUSS_WEIGHTS = SVector(0.34785484513745385, 0.6521451548625462, 0.6521451548625462, 0.34785484513745385)

include("math.jl")
include("element.jl")
include("point.jl")
include("loads.jl")
include("assembly.jl")
include("system.jl")
include("analyses.jl")
include("postprocess.jl")

end


text:
using Parquet
using Decimals
using Test

const decimal_encoding_testdata = [
    (
        precision=Int32(19),
        scale=Int32(0),
        datatype=Int128,
        byte_data=[
            UInt8[0,0,0,0,1],
            UInt8[0,0,0,1,1],
            UInt8[0,0,1,1,1],
            UInt8[0,1,1,1,1],
            UInt8[1,1,1,1,1],
        ],
        converted_data=Int128[1, 257, 65793, 16843009, 4311810305]
    ),
    (
        precision=Int32(10),
        scale=Int32(0),
        datatype=Int64,
        byte_data=[
            UInt8[0,0,0,0,1],
            UInt8[0,0,0,1,1],
            UInt8[0,0,1,1,1],
            UInt8[0,1,1,1,1],
            UInt8[1,1,1,1,1],
        ],
        converted_data=Int64[1, 257, 65793, 16843009, 4311810305]
    ),
    (
        precision=Int32(5),
        scale=Int32(0),
        datatype=Int32,
        byte_data=[
            UInt8[0,0,0,1],
            UInt8[0,0,1,1],
            UInt8[0,1,1,1],
        ],
        converted_data=Int32[1, 257, 65793]
    ),
    (
        precision=Int32(4),
        scale=Int32(0),
        datatype=Int16,
        byte_data=[
            UInt8[0,1],
            UInt8[1,1],
        ],
        converted_data=Int16[1, 257]
    ),
    (
        precision=Int32(4),
        scale=Int32(2),
        datatype=Decimal,
        byte_data=[
            UInt8[0,1],
            UInt8[1,1],
        ],
        converted_data=Decimal[Decimal(0, 1, -2), Decimal(0, 257, -2)]
    ),
    (
        precision=Int32(19),
        scale=Int32(5),
        datatype=Decimal,
        byte_data=[
            UInt8[0,0,0,0,1],
            UInt8[1,1,1,1,1],
        ],
        converted_data=Decimal[Decimal(0, 1, -5), Decimal(0, 4311810305, -5)]
    ),
    (
        precision=Int32(10),
        scale=Int32(-2),
        datatype=Decimal,
        byte_data=[
            UInt8[0,0,0,0,1],
            UInt8[0,0,0,1,1],
        ],
        converted_data=[Decimal(0, 1, 2), Decimal(0, 257, 2)]
    ),
    (
        precision=Int32(4),
        scale=Int32(2),
        datatype=Decimal,
        byte_data=Int64[
            200,
            1234,
        ],
        converted_data=Decimal[Decimal(0, 200, -2), Decimal(0, 1234, -2)]
    ),
    (
        precision=Int32(4),
        scale=Int32(2),
        datatype=Decimal,
        byte_data=Int32[
            200,
            1234,
        ],
        converted_data=Decimal[Decimal(0, 200, -2), Decimal(0, 1234, -2)]
    ),
]

function test_codec()
    @debug("testing reading bitpacked run (old scheme)")
    let data = UInt8[0x05, 0x39, 0x77]
        byte_width = Parquet.@bit2bytewidth(UInt8(3))
        typ = Parquet.@byt2itype(byte_width)
        #arr = Array{typ}(undef, 8)
        inp = Parquet.InputState(data, 0)
        out = Parquet.OutputState(typ, 8)
        Parquet.read_bitpacked_run_old(inp, out, Int32(8), UInt8(3))
        @test out.data == Int32[0:7;]
    end

    @debug("testing reading bitpacked run")
    let data = UInt8[0x88, 0xc6, 0xfa]
        byte_width = Parquet.@bit2bytewidth(UInt8(3))
        typ = Parquet.@byt2itype(byte_width)
        #arr = Array{typ}(undef, 8)
        inp = Parquet.InputState(data, 0)
        out = Parquet.OutputState(typ, 8)
        Parquet.read_bitpacked_run(inp, out, 8, UInt8(3), byte_width)
        @test out.data == Int32[0:7;]
    end

    @debug("testing decimal decoding")
    for data in decimal_encoding_testdata
        (d, f1) = Parquet.map_logical_decimal(Int32(data.precision), Int32(data.scale))
        f2 = (bytes)->Parquet.logical_decimal(bytes, data.precision, data.scale)
        f3 = (bytes)->Parquet.logical_decimal(bytes, data.precision, data.scale; use_float=true)
        @test d === data.datatype
        if isbitstype(d)
            @test all(map(f1, data.byte_data) .=== data.converted_data)
            @test all(map(f2, data.byte_data) .=== data.converted_data)
            @test all(map(f3, data.byte_data) .=== data.converted_data)
        else
            @test all(map(f1, data.byte_data) .== data.converted_data)
            @test all(map(f2, data.byte_data) .== data.converted_data)
            @test all(map(f3, data.byte_data) .== convert(Vector{Float64}, data.converted_data))
        end
    end

    iob = IOBuffer()
    write(iob, Int64(1000))
    write(iob, Int64(1234))
    is = Parquet.InputState(take!(iob), 0)
    os = Parquet.OutputState(Decimals.Decimal, 2)
    Parquet.read_plain_values(is, os, Int32(2), (v)->Parquet.logical_decimal(v, 4, 2), Parquet.PAR2._Type.INT64)
    @test os.data == Decimal[Decimal(0, 1000, -2), Decimal(0, 1234, -2)]

    iob = IOBuffer()
    write(iob, Int32(1000))
    write(iob, Int32(1234))
    is = Parquet.InputState(take!(iob), 0)
    os = Parquet.OutputState(Decimals.Decimal, 2)
    Parquet.read_plain_values(is, os, Int32(2), (v)->Parquet.logical_decimal(v, 4, 2), Parquet.PAR2._Type.INT32)
    @test os.data == Decimal[Decimal(0, 1000, -2), Decimal(0, 1234, -2)]

    iob = IOBuffer(UInt8[0,0,0,0,1,0,0,0,1,1])
    is = Parquet.InputState(take!(iob), 0)
    os = Parquet.OutputState(Decimals.Decimal, 2)
    Parquet.read_plain_values(is, os, Int32(2), (v)->Parquet.logical_decimal(v, 4, 2), Parquet.PAR2._Type.FIXED_LEN_BYTE_ARRAY)
    @test os.data == Decimal[Decimal(0, 1, -2), Decimal(0, 257, -2)]
end

@testset "codec" begin
    test_codec()
end


text:
"""
Higams' ldlt_symm translated from matlab. Performs  a so called
 BKK  Bounded Bunch Kaufman factorization of A0, that means
 ||L|| is bounded and bounded away from zero.
"""
function  ldlt_symm(A0 :: Array{Float64,2}, piv :: Char='r')

    #LDLT_SYMM  Block LDL^T factorization for a symmetric indefinite matrix.
    #     Given a Hermitian matrix A,
    #     [L, D, P, RHO, NCOMP] = LDLT_SYMM(A, PIV) computes a permutation P,
    #     a unit lower triangular L, and a real block diagonal D
    #     with 1x1 and 2x2 diagonal blocks, such that  P*A*P' = L*D*L'.
    #     PIV controls the pivoting strategy:
    #       PIV = 'p': partial pivoting (Bunch and Kaufman),
    #       PIV = 'r': rook pivoting (Ashcraft, Grimes and Lewis).
    #     The default is partial pivoting.
    #     RHO is the growth factor.
    #     For PIV = 'r' only, NCOMP is the total number of comparisons.

    #     References:
    #     J. R. Bunch and L. Kaufman, Some stable methods for calculating
    #        inertia and solving symmetric linear systems, Math. Comp.,
    #        31(137):163-179, 1977.
    #     C. Ashcraft, R. G. Grimes and J. G. Lewis, Accurate symmetric
    #        indefinite linear equation solvers. SIAM J. Matrix Anal. Appl.,
    #        20(2):513-561, 1998.
    #     N. J. Higham, Accuracy and Stability of Numerical Algorithms,
    #        Second edition, Society for Industrial and Applied Mathematics,
    #        Philadelphia, PA, 2002; chap. 11.

    #    This routine does not exploit symmetry and is not designed to be
    #     efficient.
    #
    #    Adapted from N. Higham Matlab code.
    #    JPD april 23 2015

    # Since array contents are mutable and modified,
    # copy the input matrix to keep it unchanged outside the
    # scope of the function

    A = copy(A0)

    # minimal checks for conforming inputs
    isequal(triu(A)',tril(A)) || error("Must supply Hermitian matrix.")
    piv in ['p','r'] || error("Pivoting must be 'p' or 'r'.")

    n, = size(A)

    k = 1
    D = eye(n,n);
    L = eye(n,n);
    if n == 1   D = A; end

    pp = collect(1:n)

    maxA = norm(A, Inf)
    ρ = maxA;

    ncomp = 0;
    s=1

    α = (1 + sqrt(17))/8
    while k < n
        (λ, vr) = findmax( abs.(A[k+1:n,k]) );
        r = vr[1] + k;
        if λ > 0
            swap = false;
            if abs(A[k,k]) >= α*λ
                s = 1;
            else
                if piv == 'p'
                    temp = A[k:n,r]; temp[r-k+1] = 0
                    σ = norm(temp, Inf)
                    if α*λ^2 <= abs(A[k,k])*σ
                        s = 1;
                    elseif abs(A[r,r]) >= α*σ
                        swap = true;
                        m1 = k; m2 = r;
                        s = 1;
                    else
                        swap = true;
                        m1 = k+1; m2 = r;
                        s = 2;
                    end
                    if swap
                        A[[m1,m2],:] = A[[m2,m1],:]
                        L[[m1,m2],:] = L[[m2,m1],:]
                        A[:,[m1,m2]] = A[:,[m2,m1]]
                        L[:,[m1,m2]] = L[:,[m2,m1]]

                        pp[[m1,m2]] = pp[[m2,m1]]
                    end
                elseif piv == 'r'
                    j = k;
                    pivot = false;
                    λ_j = λ;
                    while ~pivot
                        (temp1,vr) = findmax( abs.(A[k:n,j]) );
                        ncomp = ncomp + n-k;
                        r = vr[1] + k - 1;
                        temp = A[k:n,r]; temp[r-k+1] = 0.0;
                        λᵣ = max(maximum(temp), -minimum(temp))
                        ncomp = ncomp + n-k;
                        if α*λᵣ <= abs(A[r,r])
                            pivot = true;
                            s = 1;
                            A[k,:], A[r,:] = A[r,:], A[k,:]
                            L[k,:], L[r,:] = L[r,:], L[k,:]
                            A[:,k], A[:,r] = A[:,r], A[:,k]
                            L[:,k], L[:,r] = L[:,r], L[:,k]
                            pp[k], pp[r]  = pp[r], pp[k]
                        elseif λ_j == λᵣ
                            pivot = true;
                            s = 2;
                            A[k,:], A[j,:] = A[j,:], A[k,:]
                            L[k,:], L[j,:] = L[j,:], L[k,:]
                            A[:,k], A[:,j] = A[:,j], A[:,k]
                            L[:,k], L[:,j] = L[:,j], L[:,k]
                            pp[k], pp[j]  = pp[j], pp[k]
                            k1 = k+1;
                            A[k1,:], A[r,:] = A[r,:], A[k1,:]
                            L[k1,:], L[r,:] = L[r,:], L[k1,:]
                            A[:,k1], A[:,r] = A[:,r], A[:,k1]
                            L[:,k1], L[:,r] = L[:,r], L[:,k1]
                            pp[k1], pp[r]  = pp[r], pp[k1]
                        else
                            j = r;
                            s = 1
                            λ_j = λᵣ;
                        end
                    end
                end
            end

            if s == 1

                D[k,k] = A[k,k];
                A[k+1:n,k] = A[k+1:n,k]/A[k,k];
                L[k+1:n,k] = A[k+1:n,k];
                i = k+1:n;
                A[i,i] = A[i,i] - A[i,k:k] * A[k:k,i];
                A[i,i] = 0.5 * (A[i,i] + A[i,i]');

            elseif s == 2

                E = A[k:k+1,k:k+1];
                D[k:k+1,k:k+1] = E;
                i= k+2:n;
                C = A[i,k:k+1];
                temp = C/E;
                L[i,k:k+1] = temp;
                A[i,k+2:n] = A[i,k+2:n] - temp*C';
                A[i,i] = 0.5 * (A[i,i] + A[i,i]');
            end

            # Ensure diagonal real (see LINPACK User's Guide, p. 5.17).
            for i=k+s:n
                A[i,i] = real(A[i,i]);
            end

            if  k+s <= n
                val, = findmax(abs.(A[k+s:n,k+s:n]))
                ρ = max(ρ, val );
            end

        else  # Nothing to do.
            s = 1;
            D[k,k] = A[k,k];
        end

        k = k + s;

        if k == n
            D[n,n] = A[n,n];
            break
        end

    end
    #P=eye(n)
    #P = P[pp,:]
    return L, D, pp, ρ, ncomp
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
module SSLTestDataRead
using SimpleSDMLayers
using Test

l = SimpleSDMPredictor(WorldClim, BioClim, 1)

f1 = tempname()
geotiff(f1, l)
mp1 = geotiff(SimpleSDMResponse, f1)

@test typeof(mp1) <: SimpleSDMResponse
@test size(mp1) == size(l)
@test mp1 == l

f2 = tempname()
geotiff(f2, l; nodata=-3.4f38)
mp2 = geotiff(SimpleSDMPredictor, f2)

@test typeof(mp2) <: SimpleSDMPredictor
@test size(mp2) == size(l)
@test mp2 == l

end


text:
using HealthBase
using Test

struct Foo end

@testset "HealthBase.jl" begin
    include("smart_authorization.jl")
end


text:
module LMDB

    @isdefined(Docile) && eval(:(@document))

    import Base: open, close, getindex, setindex!, put!, reset,
                 isopen, count, delete!, keys, get, show, convert, show
    import Base.Iterators: drop

    export Environment, create, open, close, sync, set!, unset!, getindex, setindex!, path, info, show,
           Transaction, start, abort, commit, reset, renew, environment,
           DBI, drop, delete!, keys, get, put!,
           Cursor, count, transaction, database,
           isflagset, isopen,
           LMDBError, CursorOps, LMDBDict

    """LMDB exception type"""
    struct LMDBError <: Exception
        code::Cint
        msg::AbstractString
        LMDBError(code::Cint) = new(code, errormsg(code))
    end
    show(io::IO, err::LMDBError) = print(io, "Code[$(err.code)]: $(err.msg)")
    function checked_call(f,args...)
        ret = f(args...)
        ret === zero(Cint) || throw(LMDBError(ret))
        ret
    end

    include("LibLMDB.jl")
    using .LibLMDB
    include("common.jl")
    include("env.jl")
    include("txn.jl")
    include("dbi.jl")
    include("cur.jl")
    include("dicts.jl")
end


text:
export HRay, RandomHRay, get_vertex, point_on_ray

"""
`HRay(P::HPoint, t::Real)` returns a ray with vertex `P`
pointing to `exp(im*t)`.
"""
struct HRay <: HObject
    pt::HPoint
    t::Float64
    attr::Dict{Symbol,Any}
    function HRay(P::HPoint, theta::Real)
        th = mod(theta,2pi)
        PP = HPoint(getz(P))
        R = new(PP,th,Dict{Symbol,Any}())
        set_color(R)
        set_thickness(R)
        set_line_style(R)
        return R
    end
end
HRay(R::HRay) = HRay(R.pt, R.t)  # copy constructor
HRay() = HRay(HPoint(),0.0)

"""
`HRay(A,B)` where `A` and `B` are points creates the ray
with vertex `A` passing through `B`.
"""
function HRay(A::HPoint, B::HPoint)
    @assert A!=B "Points must be distinct to define a ray"
    f = move2xplus(A,B)
    g = inv(f)
    z = g(1)
    t = angle(z)
    return HRay(A,t)
end

"""
`RandomHRay()` creates a random ray. `RandomRay(P::HPoint)` creates a
random ray with vertex `P`.
"""
RandomHRay(P::HPoint) = HRay(P, 2*pi*rand())
RandomHRay() = RandomHRay( RandomHPoint() )

"""
`get_vertex(R::HRay)` returns the vertex (end point) of the ray.
"""
get_vertex(R::HRay) = R.pt

function (==)(R::HRay, RR::HRay)
    R.pt == RR.pt && abs(exp(im*R.t) - exp(im*RR.t)) <= THRESHOLD * eps(1.0)
end

function show(io::IO, R::HRay)
    P = R.pt
    t = R.t
    print(io,"HRay($P,$t)")
end

function HLine(R::HRay)
    f = move2xplus(R)
    g = inv(f)
    a = g(-1)
    b = g(1)
    s = angle(a)
    t = angle(b)
    return HLine(s,t)
end

"""
`point_on_ray(R::HRay)` returns a point in the interior of the ray.
"""
function point_on_ray(R::HRay)
    f = move2xplus(R)
    g = inv(f)
    P = g(HPoint(1,0))
    return P
end

"""
`in(P::HPoint, R::HRay)` determine if `P` lies on the ray `R`.
"""
function in(P::HPoint, R::HRay)
    f = move2xplus(R)
    z = f(getz(P))
    return abs(imag(z)) <= THRESHOLD*eps(1.0) && real(z) >= -THRESHOLD*eps(1.0)
end

## TO DO LIST

# meet of rays with each other, lines, segments
# reflect_across


text:
"""
	sphdistance(cmd0::String="", arg1=nothing, kwargs...)

Create Voronoi distance, node, or natural nearest-neighbor grid on a sphere

Full option list at [`sphdistance`]($(GMTdoc)sphdistance.html)

Parameters
----------

- **C** | **save_mem** :: [Type => Bool]

    For large data sets you can save some memory (at the expense of more processing).
    ($(GMTdoc)sphdistance.html#a)
- **D** | **duplicates** :: [Type => Bool]

    Delete any duplicate points [Default assumes there are no duplicates].
    ($(GMTdoc)sphdistance.html#d)
- **E** | **what_quantity** :: [Type => Str]   ``Arg = d|n|z[dist]``

    Specify the quantity that should be assigned to the grid nodes.
    ($(GMTdoc)sphdistance.html#e)
- **G** | **grid** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = sphdistance(....) form.
    ($(GMTdoc)sphdistance.html#g)
- **I** | **inc** :: [Type => Str or Number]

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    ($(GMTdoc)sphdistance.html#i)
- **L** | **dist_unit** :: [Type => Str]      ``Arg = d|e|f|k|M|n|u``

    Specify the unit used for distance calculations.
    ($(GMTdoc)sphdistance.html#l)
- **N** | **nodes** :: [Type => Str]      ``Arg = nodes``

    Read the information pertaining to each Voronoi polygon (the unique node lon, lat and polygon area)
    from a separate file.
    ($(GMTdoc)sphdistance.html#n)
- **Q** | **voronoi** :: [Type => Str]     ``Arg = voronoifile``

    Append the name of a file with pre-calculated Voronoi polygons.
    ($(GMTdoc)sphdistance.html#q)
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function sphdistance(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("sphdistance ", cmd0, arg1, arg2)

	d = KW(kwargs)
	help_show_options(d)		# Check if user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:I :R :V_params :b :d :e :h :i :r :yx])
	cmd  = parse_these_opts(cmd, d, [[:C :save_mem], [:D :duplicates], [:E :what_quantity], [:G :grid :outgrid], [:L :dist_unit]])
	cmd, arg1, arg2 = parse_QN_sphdst(d, [[:Q :voronoi], [:N :nodes]], cmd, arg1, arg2)

	common_grd(d, "sphdistance " * cmd, arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
function parse_QN_sphdst(d::Dict, symbs::Array{Array{Symbol,2},1}, cmd::String, arg1, arg2)
	(show_kwargs[1]) && print_kwarg_opts(symbs[1], "Array{GMTdataset} | GMTdataset | Array | Number | String")
	(show_kwargs[1]) && return print_kwarg_opts(symbs[3], "Array{GMTdataset} | GMTdataset | Array | Number | String")
	N_used = (arg1 === nothing) ?  0 : 1
	flags ="QN"			# Process option -Q & -N
	for k = 1:2
		if ((val = find_in_dict(d, symbs[k])[1]) !== nothing)
			cmd *= " -" * flags[k]
			if (isa(val, GMTdataset) || isa(val, Array{<:GMTdataset}) || (isa(val, Array{<:Number}) && k == 2) )
				(N_used == 0) ? arg1 = val : arg2 = val
			else
				cmd *= arg2str(val)
			end
		end
	end
	return cmd, arg1, arg2
end

# ---------------------------------------------------------------------------------------------------
sphdistance(arg1; kw...) = sphdistance("", arg1; kw...)