text:
# using MPI
using QXContexts
using ArgParse
using TimerOutputs

"""
    parse_commandline(ARGS)

Parse command line arguments and return argument dictionary
"""
function parse_commandline(ARGS)
    s = ArgParseSettings("QXContexts")

    @add_arg_table! s begin
        "--dsl", "-d"
            help = "DSL file path"
            required = true
            arg_type = String
        "--parameter-file", "-p"
            help = "Parameter file path, default is to use dsl filename with .yml suffix"
            default = nothing
            arg_type = String
        "--input-file", "-i"
            help = "Input data file path, default is to use dsl filename with .jld2 suffix"
            default = nothing
            arg_type = String
        "--output-file", "-o"
            help = "Output data file path"
            required = true
            arg_type = String
        "--number-amplitudes", "-a"
            help = "Number of amplitudes to calculate out of number in parameter file"
            default = nothing
            arg_type = Int
        "--number-slices", "-n"
            help = "The number of slices to use out of number given in parameter file"
            default = nothing
            arg_type = Int
        "--sub-comm-size", "-s"
            help = "The number of ranks to assign to each sub-communicator for partitions"
            default = 1
            arg_type = Int
        "--mpi", "-m"
            help = "Use MPI"
            action = :store_true
        "--gpu", "-g"
            help = "Use GPU if available"
            action = :store_true
        "--timings", "-t"
            help = "Enable output of timings with warmup run"
            action = :store_true
    end
    return parse_args(ARGS, s)
end

"""
    main(ARGS)

QXContexts entry point
"""
function main(ARGS)
    args = parse_commandline(ARGS)

    dsl_file       = args["dsl"]
    input_file     = args["input-file"]
    param_file = args["parameter-file"]
    output_file    = args["output-file"]
    number_amplitudes = args["number-amplitudes"]
    number_slices  = args["number-slices"]
    sub_comm_size  = args["sub-comm-size"]
    use_mpi        = args["mpi"]
    use_gpu        = args["gpu"]
    timings        = args["timings"]

    results = execute(dsl_file, input_file, param_file, output_file;
                      use_mpi=use_mpi, sub_comm_size=sub_comm_size,
                      use_gpu=use_gpu, max_amplitudes=number_amplitudes,
                      max_slices=number_slices,
                      timings=timings)

    if timings
        reset_timer!()
        results = execute(dsl_file, input_file, param_file, output_file;
                          use_mpi=use_mpi, sub_comm_size=sub_comm_size,
                          use_gpu=use_gpu, max_amplitudes=number_amplitudes,
                          max_slices=number_slices,
                          timings=timings)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end


text:
function _add_parameters_to_url_query_string(base::String, options::Dict{String,Any})::String

    # init -
    url_string = base

    parameters = ""
    for (key, value) in options
        parameters *= "$(key)=$(value)&"
    end

    # cut off trailing &
    query_parameters = parameters[1:end-1]

    # return -
    return url_string * query_parameters
end

function model(userModelType::Type{PQPolygonSDKUserModel},
    options::Dict{String,Any})::PQPolygonSDKUserModel

    # initialize -
    user_model = PQPolygonSDKUserModel()

    # for the result of the fields, let's lookup in the dictionary.
    # error state: if the dictionary is missing a value -
    for field_name_symbol ∈ fieldnames(userModelType)

        # convert the field_name_symbol to a string -
        field_name_string = string(field_name_symbol)

        # check the for the key -
        if (haskey(options,field_name_string) == false)
            throw(ArgumentError("dictionary is missing: $(field_name_string)"))
        end

        # get the value -
        value = options[field_name_string]

        # set -
        setproperty!(user_model,field_name_symbol,value)
    end

    # return -
    return user_model
end

function model(apiModelType::Type{T},
    userModel::PQPolygonSDKUserModel,
    options::Dict{String,Any})::AbstractPolygonEndpointModel where T<:AbstractPolygonEndpointModel

    # initialize -
    model = eval(Meta.parse("$(apiModelType)()")) # empty endpoint model -

    # grab the apikey from the user object -
    apikey = userModel.apikey
    model.apikey = apikey

    # for the result of the fields, let's lookup in the dictionary.
    # error state: if the dictionary is missing a value -
    for field_name_symbol ∈ fieldnames(apiModelType)

        # skip apikey - we already set that
        if (field_name_symbol != :apikey)

            # convert the field_name_symbol to a string -
            field_name_string = string(field_name_symbol)

            # check the for the key -
            if (haskey(options,field_name_string) == false)
                throw(ArgumentError("dictionary is missing: $(field_name_string)"))
            end

            # get the value -
            value = options[field_name_string]

            # set -
            setproperty!(model,field_name_symbol,value)
        end
    end

    # return -
    return model
end

# -- URL FACTORY METHODS BELOW HERE --------------------------------------------------- #
function url(base::String, model::PolygonTickerNewsEndpointModel;
    apiversion::Int = 2)::String

    # get data from the API call data -
    apikey = model.apikey
    ticker = model.ticker
    limit = model.limit
    order = model.order

    # check: do we have a published_utc parameter?
    published_utc = nothing
    if (isnothing(model.published_utc) == false)
        published_utc = model.published_utc
    end

    # check: do we have a sort parameter?
    sort = nothing
    if (isnothing(model.sort) == false)
        sort = model.sort
    end

    # build up the base string -
    base_url = "$(base)/v$(apiversion)/reference/news?"

    # what keys are passed as parameters?
    options_dictionary = Dict{String,Any}()
	options_dictionary["apiKey"] = apikey
    options_dictionary["limit"] = limit
    options_dictionary["order"] = order
    options_dictionary["ticker"] = ticker

    # do we have sort data?
    if (isnothing(sort) == false)
        options_dictionary["order"] = sort
    end

    # do we have published_utc data?
    if (isnothing(published_utc) == false)
        options_dictionary["published_utc"] = published_utc
    end

    # return -
    return _add_parameters_to_url_query_string(base_url, options_dictionary)
end

function url(base::String, model::PolygonDailyOpenCloseEndpointModel;
    apiversion::Int = 1)::String

    # get data from the API call data -
    adjusted = model.adjusted
    apikey = model.apikey
    date = model.date
    ticker = model.ticker

    # build up the base string -
    base_url = "$(base)/v$(apiversion)/open-close/$(ticker)/$(date)?"

    # what keys are passed as parameters?
    options_dictionary = Dict{String,Any}()
	options_dictionary["adjusted"] = adjusted
	options_dictionary["apiKey"] = apikey

    # return -
    return _add_parameters_to_url_query_string(base_url, options_dictionary)
end

function url(base::String, model::PolygonGroupedDailyEndpointModel;
    apiversion::Int = 2)::String

    # get data from the API call data -
    adjusted = model.adjusted
    apikey = model.apikey
    date = model.date

    # build up the base string -
    base_url = "$(base)/v$(apiversion)/aggs/grouped/locale/us/market/stocks/$(date)?"

    # what keys are passed as parameters?
    options_dictionary = Dict{String,Any}()
	options_dictionary["adjusted"] = adjusted
	options_dictionary["apiKey"] = apikey

    # return -
    return _add_parameters_to_url_query_string(base_url, options_dictionary)
end

function url(base::String, model::PolygonPreviousCloseEndpointModel;
    apiversion::Int = 2)::String

    # get data from the API call data -
    adjusted = model.adjusted
    apikey = model.apikey
    ticker = model.ticker

    # build up the base string -
    base_url = "$(base)/v$(apiversion)/aggs/ticker/$(ticker)/prev?"

    # what keys are passed as parameters?
    options_dictionary = Dict{String,Any}()
	options_dictionary["adjusted"] = adjusted
	options_dictionary["apiKey"] = apikey

    # return -
    return _add_parameters_to_url_query_string(base_url, options_dictionary)
end

function url(base::String, model::PolygonAggregatesEndpointModel;
    apiversion::Int = 2)::String

    # get data from the API call data -
    adjusted = model.adjusted
    limit = model.limit
    sortdirection = model.sortdirection
    apikey = model.apikey
    ticker = model.ticker
    to = model.to
    from = model.from
    multiplier = model.multiplier
    timespan = model.timespan

    # build up the base string -
    base_url = "$(base)/v$(apiversion)/aggs/ticker/$(ticker)/range/$(multiplier)/$(timespan)/$(from)/$(to)?"

    # what keys are passed as parameters?
    options_dictionary = Dict{String,Any}()
	options_dictionary["adjusted"] = adjusted
	options_dictionary["sort"] = sortdirection
	options_dictionary["limit"] = limit
	options_dictionary["apiKey"] = apikey

    # return -
    return _add_parameters_to_url_query_string(base_url, options_dictionary)
end

function url(base::String, model::PolygonTickerDetailsEndpointModel;
    apiversion::Int = 1)::String

    # get data from the API call data -
    apikey = model.apikey
    ticker = model.ticker

    # build up the base string -
    base_url = "$(base)/v$(apiversion)/meta/symbols/$(ticker)/company?"

    # what keys are passed as parameters?
    options_dictionary = Dict{String,Any}()
    options_dictionary["apiKey"] = apikey

    # return -
    return _add_parameters_to_url_query_string(base_url, options_dictionary)
end
# -- URL FACTORY METHODS ABOVE HERE --------------------------------------------------- #

ext: jl


text:
using Agents
using Agents.Models: Wolf, Sheep
using BenchmarkTools
using Test

a = @benchmark step!(model, agent_step!, model_step!, 500) setup = (
    (model, agent_step!, model_step!) = Models.predator_prey(
        n_wolves = 40,
        n_sheep = 60,
        dims = (25, 25),
        Δenergy_sheep = 5,
        Δenergy_wolf = 13,
        sheep_reproduce = 0.2,
        wolf_reproduce = 0.1,
        regrowth_time = 20,
    )
) samples = 100 teardown = (@test count(i -> i isa Sheep, allagents(model)) > 0 &&
       count(i -> i isa Wolf, allagents(model)) > 0)
println("Agents.jl WolfSheep (ms): ", minimum(a.times) * 1e-6)

a = @benchmark step!(model, agent_step!, model_step!, 100) setup = (
    (model, agent_step!, model_step!) = Models.flocking(
        n_birds = 300,
        separation = 1,
        cohere_factor = 0.03,
        separate_factor = 0.015,
        match_factor = 0.05,
    )
)
println("Agents.jl Flocking (ms): ", minimum(a.times) * 1e-6)

a = @benchmark step!(model, agent_step!, model_step!, 10) setup = (
    (model, agent_step!, model_step!) =
        Models.schelling(griddims = (50, 50), numagents = 2000)
) samples = 100
println("Agents.jl Schelling (ms): ", minimum(a.times) * 1e-6)

a = @benchmark step!(model, agent_step!, model_step!, 100) setup =
    ((model, agent_step!, model_step!) = Models.forest_fire())
println("Agents.jl ForestFire (ms): ", minimum(a.times) * 1e-6)


text:
@testset "exports" begin
    # Ensure that every exported name is actually defined
    for name in names(RigidBodyDynamics)
        @test isdefined(RigidBodyDynamics, name)
    end
end


text:
include("ukbb.jl")
include("estimators.jl")
include("utils.jl")
include("models.jl")


text:
# src/transactions.jl
# Transactions for Wallets

export AbstractTransaction,
       AbstractSendTransaction,
       AbstractDestroyTransaction,
       statetype,
       stateof,
       associatedfee,
       history,
       transactions_between,
       withdrawls,
       deposits,
       build_send_transaction,
       build_destroy_transaction,
       sign!,
       execute!,
       cancel!,
       rollback!,
       cancelback!,
       is_pending,
       has_executed,
       was_completed,
       was_cancelled,
       failed,
       send!,
       destroy!

"""```
AbstractTransaction{S}
```
Abstraction for Transactions with state type `S`.
"""
abstract type AbstractTransaction{S} end

"""```
AbstractSendTransaction{S} <: AbstractTransaction{S}
```
Abstraction for Transactions which transfer assets.
"""
abstract type AbstractSendTransaction{S} <: AbstractTransaction{S} end

"""```
AbstractDestroyTransaction{S} <: AbstractTransaction{S}
```
Abstraction for Transactions which destroy assets.
"""
abstract type AbstractDestroyTransaction{S} <: AbstractTransaction{S} end

"""```
statetype(::Type{AbstractTransaction{S}}) === S
```
Return the state type `S` of the Transaction abstraction.
"""
statetype(::Type{AbstractTransaction{S}}) where {S} = S

"""```
statetype(::AbstractTransaction{S}) === S
```
Return the state type `S` of the Transaction abstraction.
"""
statetype(transaction::AbstractTransaction) = statetype(typeof(transaction))

"""```
stateof(transaction::AbstractTransaction{S})::S
```
Return the current state (of type `S`) of the transaction.
"""
function stateof(transaction::AbstractTransaction{S})::S where {S}
    missing_api("stateof", transaction)
end

"""```
associatedfee(transaction::AbstractTransaction)
```
Return the fee associated to the given `transaction`.
"""
function associatedfee(transaction::AbstractTransaction)::AbstractDestroyTransaction
    missing_api("associatedfee", transaction)
end

"""```
history(wallet::AbstractWallet)
```
Return the transaction history of the given `wallet`.
"""
function history(wallet::AbstractWallet)
    missing_api("history", wallet)
end

"""```
transactions_between(source::AbstractWallet, target::AbstractWalletOrAddress)
```
Return all transactions between `source` and `target`.
"""
function transactions_between(source::AbstractWallet, target::AbstractWalletOrAddress)
    missing_api("transactions_between", source, target)
end

"""```
transactions_between(source::AbstractWalletOrAddress, target::AbstractWallet)
```
Return all transactions between `source` and `target`.
"""
function transactions_between(source::AbstractWalletOrAddress, target::AbstractWallet)
    missing_api("transactions_between", source, target)
end

"""```
withdrawls(wallet::AbstractWalletOrAddress, target::AbstractWallet)
```
Return all send transactions to `target` from `wallet`.
"""
function withdrawls(wallet::AbstractWallet, target::AbstractWalletOrAddress)
    missing_api("withdrawls", wallet, target)
end

"""```
withdrawls(wallet::AbstractWallet)
```
Return all withdrawls out of `wallet`.
"""
function withdrawls(wallet::AbstractWallet)
    missing_api("withdrawls", wallet)
end

"""```
deposits(wallet::AbstractWalletOrAddress, source::AbstractWallet)
```
Return all send transactions from `source` to `wallet`.
"""
function deposits(wallet::AbstractWallet, source::AbstractWalletOrAddress)
    missing_api("deposits", wallet, source)
end

"""```
deposits(wallet::AbstractWallet)
```
Return all deposits into `wallet`.
"""
function deposits(wallet::AbstractWallet)
    missing_api("deposits", wallet)
end

"""```
build_send_transaction(source::AbstractWallet, target, asset, quantity;
                       fulfill_maximum=false)
```
Prepare transaction which will send a certain amount of the given `asset` to
the `target` if it is contained in this `quantity` in the original `wallet`.
If `fulfill_maximum` is `true`, then fulfill as much of the transaction as
possible.
"""
function build_send_transaction(
    source::AbstractWallet{K,N},
    target::AbstractWalletOrAddress,
    asset::K,
    quantity::N;
    fulfill_maximum = false,
)::AbstractSendTransaction where {K,N}
    missing_api(
        "build_send_transaction",
        source,
        target,
        asset,
        quantity;
        fulfill_maximum = fulfill_maximum,
    )
end

"""```
build_destroy_transaction(wallet::AbstractWallet, asset, quantity;
                          fulfill_maximum=false)
```
Prepare transaction which will destroy a certain amount of the given `asset` if
it is contained in this `quantity` in the original `wallet`. If
`fulfill_maximum` is `true`, then fulfill as much of the transaction as
possible.
"""
function build_destroy_transaction(
    wallet::AbstractWallet{K,N},
    asset::K,
    quantity::N;
    fulfill_maximum = false,
)::AbstractDestroyTransaction where {K,N}
    missing_api(
        "build_destroy_transaction",
        wallet,
        asset,
        quantity,
        fulfill_maximum = fulfill_maximum,
    )
end

"""```
sign!(transaction::AbstractTransaction)
```
Sign the given `transaction`.
"""
function sign!(transaction::AbstractTransaction, parties)
    missing_api("sign!", transaction, parties)
end

"""```
execute!(transaction::AbstractTransaction)
```
Execute the given `transaction` and wait for a return signal.
"""
function execute!(transaction::AbstractTransaction)
    missing_api("execute!", transaction)
end

"""```
cancel!(transaction::AbstractTransaction)
```
Cancel a pending transaction.
"""
function cancel!(transaction::AbstractTransaction)
    missing_api("cancel!", transaction)
end

"""```
rollback!(transaction::AbstractTransaction)
```
Rollback a completed transaction.
"""
function rollback!(transaction::AbstractTransaction)
    missing_api("rollback!", transaction)
end

"""```
cancelback!(transaction::AbstractTransaction; onsuccess, onfailure)
```
Rollback a fulfilled transaction or cancel a pending one.
"""
function cancelback!(
    transaction::AbstractTransaction;
    onsuccess = identity,
    onfailure = identity,
)
    if was_completed(transaction)
        return rollback!(transaction)
    elseif is_pending(transaction)
        return cancel!(transaction)
    elseif was_cancelled(transaction)
        return onsuccess(transaction)
    else
        return onfailure(transaction)
    end
end

"""```
is_pending(transaction::AbstractTransaction)::Bool
```
Check if `transaction` has been executed but is pending.
"""
function is_pending(transaction::AbstractTransaction)::Bool
    missing_api("is_pending", transaction)
end

"""```
has_executed(transaction::AbstractTransaction)::Bool
```
Check if `transaction` has been executed.
"""
function has_executed(transaction::AbstractTransaction)::Bool
    missing_api("has_executed", transaction)
end

"""```
was_completed(transaction::AbstractTransaction)::Bool
```
Check if `transaction` was fulfilled.
"""
function was_completed(transaction::AbstractTransaction)::Bool
    missing_api("was_completed", transaction)
end

"""```
was_cancelled(transaction::AbstractTransaction)::Bool
```
Check if `transaction` was forcibly cancelled.
"""
function was_cancelled(transaction::AbstractTransaction)::Bool
    missing_api("was_cancelled", transaction)
end

"""```
failed(transaction::AbstractTransaction)::Bool
```
Check if `transaction` failed.
"""
function failed(transaction::AbstractTransaction)::Bool
    missing_api("failed", transaction)
end

"""```
send!(source::AbstractWallet, target, asset, quantity;
      fulfill_maximum=false, check, onfailure)
```
Send a certain amount of the given asset to a target if it is contained in this
quantity in the original wallet. If `fulfill_maximum` is `true`, then fulfill
as much of the transaction as possible.
"""
function send!(
    source::AbstractWallet{K,N},
    target::AbstractWalletOrAddress,
    asset::K,
    quantity::N;
    fulfill_maximum = false,
    check = x -> true,
    onfailure = identity,
) where {K,N}
    transaction = build_send_transaction(
        source,
        target,
        asset,
        quantity;
        fulfill_maximum = fulfill_maximum,
    )
    if !check(transaction)
        return onfailure(transaction)
    end
    return execute!(transaction)
end

"""```
destroy!(wallet::AbstractWallet, asset, quantity;
         fulfill_maximum=false, check, onfailure)
```
Destroy a certain amount of the given asset if it is contained in the wallet.
If `fulfill_maximum` is `true`, then fulfill as much of the transaction as
possible.
"""
function destroy!(
    wallet::AbstractWallet{K,N},
    asset::K,
    quantity::N;
    fulfill_maximum = false,
    check = x -> true,
    onfailure = identity,
) where {K,N}
    transaction = build_destroy_transaction(
        wallet,
        asset,
        quantity;
        fulfill_maximum = fulfill_maximum,
    )
    if !check(transaction)
        return onfailure(transaction)
    end
    return execute!(transaction)
end


text:
function initstate(domain::CompiledDomain, problem::GenericProblem)
    # Construct state using interpreter, then convert to compiled
    generic_state = initstate(get_source(domain), problem)
    compiled_state = statetype(domain)(generic_state)
    return compiled_state
end

function initstate(domain::CompiledDomain, objtypes::AbstractDict, fluents)
    # Construct state using interpreter, then convert to compiled
    generic_state = initstate(get_source(domain), objtypes, fluents)
    compiled_state = statetype(domain)(generic_state)
    return compiled_state
end


text:
# Test for solving div-grad system both complex and real and single and multiple rhs
using MUMPSjInv
using Test

include("getDivGrad.jl")
A = getDivGrad(32, 32, 16)
n = size(A, 1)

# REAL: test for single rhs
println("Test for real SPD matrix: one rhs")
rhs = randn(n)
x = solveMUMPS(A, rhs, 1)
err  =  norm(A * x - rhs) / norm(rhs)
@test err < 1e-14

# REAL: test for multiple rhs
println("Test for real SPD matrix: multiple rhs")
nrhs = 10
rhs = randn(n, nrhs)
x = solveMUMPS(A, rhs, 1)
err = zeros(nrhs)
for i= 1:nrhs
    err[i] =  norm(A * x[:, i] - rhs[:, i]) / norm(rhs[:, i])
end
@test eltype(x) == Float64
@test maximum(err) < 1e-14

# COMPLEX: test for multiple rhs
println("Test for complex symmetric matrix: one rhs")
r = rand(n)
A = A + im * spdiagm(0 => r)
rhs = randn(n) + im * randn(n)
x = solveMUMPS(A, rhs,1)
err=  norm(A * x - rhs) / norm(rhs)
@test eltype(x) == ComplexF64
@test err < 1e-14

# COMPLEX : test for multiple rhs
println("Test for complex symmetric matrix: multiple rhs")
nrhs = 10
rhs = randn(n, nrhs) + im * randn(n, nrhs)
x = solveMUMPS(A, rhs, 2)
err = zeros(nrhs)
for i = 1:nrhs
    err[i] = norm(A * x[:, i] - rhs[:, i]) / norm(rhs[:, i])
end
@test eltype(x) == ComplexF64
@test maximum(err) < 1e-14


text:
"""
    lyapunov_flow(vectorfield, r0, parameters)
Integrates the equation of motion of a given vector field with the integrator passed through the parameters.
"""
function lyapunov_flow(field::Function, r0::Vector{Float64}, p::Parameters)
    method = p.integrator.f
    t = 0.0:p.dt:p.dtsampling

    function extendedfield(time::Float64, r::Vector{Float64})
        field(r, p.potential, 1/p.T, p.thermo)
    end

    try
        t,pos = method(extendedfield, r0, t, points=:specified)
        return t, pos
    catch
        t,pos = method(extendedfield, r0, t)
        return t,pos
    end

end

"""
Return the lyapunovspectrum of the passed vector field computed under the conditions passed in the parameters
"""
function lyapunov_spectrum(field::Function, r::Vector{Float64}, p::Parameters, burninsteps = 100)
    n = 3 ##Dimension of the dynamical system considered
    w = eye(n)
    norms = zeros(p.nsteps,n)

    for i in 1:burninsteps
        r = lyapunov_flow(field, r, p)[2][end]
    end

    for i in 1:p.nsteps
        r = lyapunov_flow(field, r, p)[2][end]
        u = reshape(r[4:end],n,n)
        w = gramschmidt(u')
        for j in 1:n
            norms[i,j] = norm(w[:,j])
        end
        for j in 1:n
            w[:,j] = w[:,j]/norm(w[:,j])
        end
        r[4:end] = copy(w'[:])
    end
    exps = zeros(n)
    for k in 1:n
        exps[k] = sum(log(norms[:,k]))/(p.nsteps*p.dtsampling)
    end

    return exps
end

"""
Return the results of the simulation. Basically, the results are of two types: either the lyapunov spectra for a certain number of random initial conditions or the points belonging to the trajectory defined by a random initial condition.
"""
function lyapunov_simulation(p::Parameters)
    beta = 1./p.T
    r = zeros(12)
    init = initialcondition(beta,p.Q)
    r[1:3] = copy(init)
    r[4] = r[8] = r[12] = 1.0  #Entries of the identity matrix in the initial condition
    exps = lyapunov_spectrum(variationalDDfield,r,p)
    return vcat(exps, init)
end

max_stars_repo_name: dapias/STODs

max_issues_repo_name: dapias/STODs

max_forks_repo_name: dapias/STODs


text:
function andersonmix!( vin::Array{Float64,1}, vout::Array{Float64,1},
                       beta::Float64, df, dv, iter::Int64, mixdim::Int64)
    # output [vnew,df,dv]
    Npoints = size(vin)[1]

    # function evaluation overwrites vout,
    # use dvout instead
    dvout = vout - vin

    iterused = min(iter-1,mixdim)
    ipos = iter - 1 - floor(Int64, (iter-2)/mixdim)*mixdim

    #@printf("ipos = %d\n", ipos)

    if iter > 1
        # compute the changes in function evaluations and the step (changes in potentials)
        df[:,ipos] = df[:,ipos] - dvout
        dv[:,ipos] = dv[:,ipos] - vin
    end
    #
    vinsave  = copy(vin)
    dvoutsave = copy(dvout)

    if (iter > 1)
        #B = (df(:,1:iterused))'*(df(:,1:iterused));
        #B = (B+B')/2;
        #
        #work = df(:,1:iterused)'*reshape(vout,n123,1);
        #gammas = B\(work);
        #
        # the use of pseudo-inverse is not the most efficient
        # implementation.
        #
        gammas = pinv(df[:,1:iterused])*dvout
        #println("gammas = ", gammas[:])
        #
        # the following loop can be replaced by two simple gemvs
        # and some reshape operations.
        #
        for i = 1:iterused
            vin  = vin  - gammas[i] * dv[:,i]
            dvout = dvout - gammas[i] * df[:,i]
        end
    end

    inext = iter - floor( Int64, (iter - 1) / mixdim) * mixdim
    #@printf("inext = %d\n", inext)

    df[:,inext] = dvoutsave
    dv[:,inext] = vinsave

    return vin[:] + beta*dvout[:]
end


text:
using ForwardDiff:jacobian

const MC = 100000           # maximum of counts
const IVAL = 1e-2           # first variation
const RATE = 1e-3           # variation rate
const NEPS = 1e-12          # eps of Newton's method

function create_diffeq(model_path::String)
    lines::Vector{String} = []
    open(joinpath(model_path, "set_model.jl"), "r") do f
        append!(lines, readlines(f))
    end
    for (i, line) in enumerate(lines)
        if occursin("function diffeq!", line)
            lines[i] = "function diffeq(u::Vector)\n    du = similar(u)\n\n"
        elseif line == "end"
            lines[i] = "    return du\nend"
            lines = lines[1:i]
            break
        end
    end
    open(joinpath(model_path, "forwarddiff.jl"), "w") do f
        for line in lines
            write(f, line * "\n")
        end
    end
end

# matrix transformation (large diagonal elements move to upper)
function pivoting!(s::Matrix{Float64}, pivot::Int, dim_newton::Int)
    v0::Vector{Float64} = zeros(dim_newton + 1)
    v1::Vector{Float64} = zeros(dim_newton + 1)
    possess::Int = 0
    max_element::Float64 = 0.0

    for i in pivot:size(s, 1)
        current_element = abs(s[i, pivot])
        if max_element <= current_element
            max_element = current_element
            possess = i
        end
    end

    for j in 1:size(s, 2)
        v0[j] = s[possess, j]
        v1[j] = s[pivot, j]
    end

    for j in 1:size(s, 2)
        s[possess, j] = v1[j]
        s[pivot, j] = v0[j]
    end
end

# Gaussian elimination (row reduction)
function gaussian_elimination!(s::Matrix{Float64}, e::Vector{Float64}, dim_newton::Int)
    for i in 1:dim_newton
        pivoting!(s, i, dim_newton)
    end
    # forward
    for k in 1:size(s, 1)
        w = (s[k, k] != 0.0) ? 1.0 / s[k, k] : 1.0
        for j in k:size(s, 2)
            s[k, j] *= w
            for i in k:size(s, 1)
                s[i, j] -= s[i, k] * s[k, j]
            end
        end
    end
    # backward
    for i in size(s, 1):-1:1
        sum = 0.0
        for j in i:length(e)
            sum += s[i, j] * e[j]
        end
        e[i] = s[i, end] - sum
    end
end

# Newton's method
function newtons_method!(
        diffeq::Function,
        get_derivatives::Function,
        x::Vector{Float64},
        real_part::Vector{Float64},
        imaginary_part::Vector{Float64},
        fix_num::Int,
        p::Vector{Float64},
        successful::Bool,
        bifparam::Int,
        n_state::Int,
        dim_newton::Int,
        n_variable::Int)
    u::Vector{Float64} = zeros(n_state)
    vx::Vector{Float64} = zeros(dim_newton)
    s::Matrix{Float64} = zeros(dim_newton, dim_newton + 1)

    for i in eachindex(x)
        if fix_num == i
            for j in eachindex(vx)
                idx = i + j
                if idx > length(x)
                    idx -= length(x)
                end
                vx[j] = x[idx]
            end
            break
        else
            continue
        end
    end

    # initial error
    e::Vector{Float64} = zeros(dim_newton)
    error::Float64 = 1.0

    while error > NEPS
        for i in 1:n_variable
            if fix_num == i
                idx_param = n_variable - i
                p[bifparam] = (idx_param == 0) ? x[fix_num] : vx[idx_param]
                for j in eachindex(u)
                    idx = j - i
                    if idx == 0
                        u[j] = x[fix_num]
                    elseif idx < 0
                        u[j] = vx[n_variable + idx]
                    else
                        u[j] = vx[idx]
                    end
                end
                break
        else
                continue
        end
        end

        # initialization
        dFdx::Matrix{Float64} = jacobian(diffeq, u)
        dFdp::Vector{Float64} = get_derivatives(u, p)

        F::Vector{Float64} = diffeq(u)

        eigenvalues::Array{Complex{Float64},1} = eigvals(dFdx)
        for (i, eigenvalue) in enumerate(eigenvalues)
            real_part[i] = real(eigenvalue)
            imaginary_part[i] = imag(eigenvalue)
        end

        # s = [dF-F]
        for i in 1:n_variable
            if fix_num == i
                for k in 1:n_state
                    for j in 1:n_state
                        idx = i + j
                        if idx == n_variable
                            s[k, j] = dFdp[k]
                        elseif idx > n_variable
                            s[k, j] = dFdx[k, idx - n_variable]
                        else
                            s[k, j] = dFdx[k, idx]
                        end
                    end
                    s[k, n_variable] = -F[k]
                end
                break
            else
                continue
            end
        end

        gaussian_elimination!(s, e, dim_newton)

        # update error
        error = 0.0
        @inbounds for i in eachindex(e)
            vx[i] += e[i]
            error += e[i] * e[i]
        end
        error = sqrt(error)
        if isnan(error) || isinf(error)
            successful = false
            break
        end
    end

    for i in eachindex(x)
        if fix_num == i
            for j in eachindex(vx)
                idx = i + j
                if idx > length(x)
                    idx -= length(x)
                end
                x[idx] = vx[j]
            end
            break
        else
            continue
        end
    end
end

function new_curve!(
        model_path::Union{String,SubString{String}},
        p::Vector{Float64},
        diffeq::Function,
        get_derivatives::Function,
        get_steady_state::Function;
        direction::Bool=false,
        bifparam::Int,
        n_state::Int,
        n_param::Int=1,
        n_variable::Int=n_state + 1,
        dim_newton::Int=n_state)
    # bifparam : name(index) of bifurcation parameter
    # n_state : num of state variables
    # n_param : num of parameters
    # n_variable : num of variables
    # dim_newton : dim of Newton's method
    count::Int = 1
    x::Vector{Float64} = zeros(n_variable)
    dx::Vector{Float64} = zeros(n_variable)

    real_part::Vector{Float64} = zeros(n_state)
    imaginary_part::Vector{Float64} = zeros(n_state)

    # file
    if !isdir(
        joinpath(
            model_path,
            "data",
        )
    )
        mkdir(
            joinpath(
                model_path,
                "data",
            )
        )
    else
        files::Vector{String} = readdir(
            joinpath(
                model_path,
                "data",
            )
        )
        for file in files
            rm(
                joinpath(
                    model_path,
                    "data",
                    "$file",
                )
            )
        end
    end

    FOUT1 = open(joinpath(model_path, "data", "fp.dat"), "w") # file for fixed point
    FOUT2 = open(joinpath(model_path, "data", "ev.dat"), "w") # file for eigenvalues

    # initial condition
    x[1:n_state] = get_steady_state(p)
    x[end] = p[bifparam]  # x-axis

    # initial fixed
    fix_val::Float64 = x[end]
    fix_num::Int = n_variable
    x[fix_num] = fix_val

    # first Newton's method
    successful::Bool = true
    newtons_method!(
        diffeq, get_derivatives, x, real_part, imaginary_part, fix_num, p,
        successful, bifparam, n_state, dim_newton, n_variable
    )

    write(FOUT1, @sprintf("%d\t", count))
    for i in eachindex(x)
        write(FOUT1, @sprintf("%10.8e\t", x[i]))
    end
    write(FOUT1, @sprintf("%d\n", fix_num))
    write(FOUT2, @sprintf("%d\t", count))
    for i in 1:n_state
        write(
            FOUT2, @sprintf(
                "%10.8e\t%10.8e\t", real_part[i], imaginary_part[i]
            )
        )
    end
    write(FOUT2, @sprintf("%10.8e\t%d\n", p[bifparam], fix_num))
    count += 1

    # keep optimums
    px::Vector{Float64} = copy(x)

    # variation
    fix_val += ifelse(direction, +IVAL, -IVAL)

    # same fixed variable
    x[fix_num] = fix_val

    while count <= MC && successful
        newtons_method!(
            diffeq, get_derivatives, x, real_part, imaginary_part, fix_num, p,
            successful, bifparam, n_state, dim_newton, n_variable
        )

        # maximum variation
        for (i, prev) in enumerate(px)
        @inbounds dx[i] = x[i] - prev
        end
        sum::Float64 = 0.0
        for i in eachindex(dx)
            @inbounds sum += dx[i] * dx[i]
        end
        ave::Float64 = sqrt(sum)
        for i in eachindex(dx)
            @inbounds dx[i] /= ave
        end
        px = copy(x)
        for (i, diff) in enumerate(dx)
            @inbounds x[i] += abs(RATE) * diff
        end

        # fix variable with maximum variation
        fix_num = 1
        for i in 2:length(dx)
            if abs(dx[fix_num]) < abs(dx[i])
                fix_num = i
        end
        end

            # Stop calc.
        if x[end] <= 0.0
            successful = false
        end

        write(FOUT1, @sprintf("%d\t", count))
        for i in eachindex(x)
            write(FOUT1, @sprintf("%10.8e\t", x[i]))
        end
        write(FOUT1, @sprintf("%d\n", fix_num))
        write(FOUT2, @sprintf("%d\t", count))
        for i in 1:n_state
            write(
                FOUT2, @sprintf(
                    "%10.8e\t%10.8e\t", real_part[i], imaginary_part[i]
                )
            )
        end
        write(FOUT2, @sprintf("%10.8e\t%d\n", p[bifparam], fix_num))
        count += 1
    end

    close(FOUT1)
    close(FOUT2)
end

function get_bistable_regime(ev::Matrix{Float64}, n_state::Int)
    br::Vector{Int} = []
    for i in 1:size(ev, 1)
        if maximum(ev[i, [2j for j in 1:n_state]]) > 0.0
            push!(br, i)
        end
    end
    return br
end


text:
abstract type AbstractMarkovDecisionProcess end
## --------------------------------------------- ##
## GENERIC PART
## --------------------------------------------- ##

"""
    MarkovDecisionProcess is a MDP implementation of AbstractMarkovDecisionProcess.

A Markov decision process is a sequential decision problem with fully observable and stochastic environment with a transition model and rewards function.
The discount factor (gamma variable) describes the preference for current rewards over future rewards.
"""
struct MarkovDecisionProcess{T} <: AbstractMarkovDecisionProcess
    initial::T
    states::Set{T}
    actions::Set{T}
    terminal_states::Set{T}
    transitions::Dict
    gamma::Float64
    reward::Dict

    function MarkovDecisionProcess{T}(initial::T, actions_list::Set{T}, terminal_states::Set{T}, transitions::Dict, states::Union{Nothing, Set{T}}, gamma::Float64) where T
        !(0 < gamma <= 1) ? nothing : error("MarkovDecisionProcess():\nThe gamma variable of an MDP must be between 0 and 1, the constructor was given ", gamma, "!")
        new_states = typeof(states) <: Set ? states : Set{typeof(initial)}()

        return new(initial, new_states, actions_list, terminal_states, transitions, gamma, Dict())
    end
end

function MarkovDecisionProcess(initial, actions_list::Set, terminal_states::Set, transitions::Dict;
                                states::Union{Nothing, Set}=nothing, gamma::Float64=0.9)
    return MarkovDecisionProcess{typeof(initial)}(initial, actions_list, terminal_states, transitions, states, gamma)
end

"""
    reward(mdp::T, state) where {T <: AbstractMarkovDecisionProcess}

Return a reward based on the given 'state'.
"""
function reward(mdp::T, state) where {T <: AbstractMarkovDecisionProcess}
    return mdp.reward[state]
end

"""
    transition_model(mdp::T, state, action) where {T <: AbstractMarkovDecisionProcess}

Return a list of (P(s'|s, a), s') pairs given the state 's' and action 'a'.
"""
function transition_model(mdp::T, state, action) where {T <: AbstractMarkovDecisionProcess}
    length(mdp.transitions) == 0 ? nothing : error("transition_model(): The transition model for the given 'mdp' could not be found!")

    return mdp.transitions[state][action]
end

"""
    actions(mdp::T, state) where {T <: AbstractMarkovDecisionProcess}

Return a set of actions that are possible in the given state.
"""
function actions(mdp::T, state) where {T <: AbstractMarkovDecisionProcess}
    if state in mdp.terminal_states
        return Set{Nothing}([nothing])
    else
        return mdp.actions;
    end
end

## --------------------------------------------- ##
## SPECIFIC PART 4x3 GAME
## --------------------------------------------- ##

"""
    GridMarkovDecisionProcess is a MDP implementation of the grid from chapter 17.

Obstacles in the environment are representedby a nothing.

Arguments:

"""
struct GridMarkovDecisionProcess <: AbstractMarkovDecisionProcess
    initial::Tuple{Int64, Int64}
    states::Set{Tuple{Int64, Int64}}
    actions::Set{Tuple{Int64, Int64}}
    terminal_states::Set{Tuple{Int64, Int64}}
    grid::Array{Union{Nothing, Float64}, 2}
    gamma::Float64
    reward::Dict
    function GridMarkovDecisionProcess(initial::Tuple{Int64, Int64}, terminal_states::Set{Tuple{Int64, Int64}},
                                       grid::Array{Union{Nothing, Float64}, 2};
                                       states::Union{Nothing, Set{Tuple{Int64, Int64}}}=nothing, gamma::Float64=0.9)
        (0 < gamma <= 1) ? nothing : error("MarkovDecisionProcess():\nThe gamma variable of an MDP must be between 0 and 1, the constructor was given ", gamma, "!")
        new_states = typeof(states) <: Set ? states : Set{Tuple{Int64, Int64}}()
        orientations::Set = Set{Tuple{Int64, Int64}}([(1, 0), (0, 1), (-1, 0), (0, -1)])
        reward::Dict = Dict()
        for i in 1:size(grid, 1)
            for j in 1:size(grid, 2)
                reward[(i, j)] = grid[i, j]
                if !(grid[i, j] === nothing)
                    push!(new_states, (i, j))
                end
            end
        end

        return new(initial, new_states, orientations, terminal_states, grid, gamma, reward);
    end
end

"""
    go_to(gmdp::GridMarkovDecisionProcess, state::Tuple{Int64, Int64}, direction::Tuple{Int64, Int64})

Return the next state given the current state and direction. If the next state is not known
"""
function go_to(gmdp::GridMarkovDecisionProcess, state::Tuple{Int64, Int64}, direction::Tuple{Int64, Int64})
    next_state::Tuple{Int64, Int64} = state .+ direction

    return next_state in gmdp.states ? next_state : state
end

"""
    transition_model

Return an array of (P(s'|s, a), s') pairs given the state 's' and action 'a'.
"""
function transition_model(gmdp::GridMarkovDecisionProcess, state::Tuple{Int64, Int64}, action::Nothing)
    return [(0.0, state)];
end

"""
    transition_model

Return an array of (P(s'|s, a), s') pairs given the state 's' and action 'a'.
"""
function transition_model(gmdp::GridMarkovDecisionProcess, state::Tuple{Int64, Int64}, action::Tuple{Int64, Int64})
    return [
            (0.8, go_to(gmdp, state, action)),
            (0.1, go_to(gmdp, state, turn_heading(action, -1))),
            (0.1, go_to(gmdp, state, turn_heading(action, 1)))
            ];
end

# (0, 1) will move the agent rightward.
# (-1, 0) will move the agent upward.
# (0, -1) will move the agent leftward.
# (1, 0) will move the agent downward.
const arrow_characters = Dict((0, 1) => ">", (-1, 0) => "^", (0, -1) => "<", (1, 0) => "v", nothing => ".")

"""
    show_grid(gmdp::GridMarkovDecisionProcess, mapping::Dict)

Returns a string representation of the current grid using [<, >, v, ^] to indication the heading and "." to indicate no action.
"""
function show_grid(gmdp::GridMarkovDecisionProcess, policy::Dict)
    mapping = Dict(state => arrow_characters[action] for (state, action) in policy)
    io = IOBuffer()
    for i in 1:size(gmdp.grid, 1)
        for j in 1:size(gmdp.grid, 2)
            print(io, get(mapping, (i,j), "o"))
        end
        print(io,"\n")
    end
    return String(take!(io))
end

"""
    turn_heading

Given an input heading, return the heading in inc positions further in a clockwise manner

The headings can be specified with H. The heading h should occur in the collection of headings H.
"""
function turn_heading(h::Tuple{Int64,Int64}, inc::Int64; H=[(0,1),(1,0),(0,-1),(-1,0)])
    # get index of current heading
    i = findfirst(x->x==h, H)
    # return incremented heading
    return H[(i + inc + length(H) - 1) % length(H) + 1]
end

## --------------------------------------------- ##
## Solution methods
## --------------------------------------------- ##

"""
    value_iteration(mdp::T; epsilon::Float64=0.001) where {T <: AbstractMarkovDecisionProcess}

Return the utilities of the MDP's states as a Dict by applying the value iteration algorithm (Fig. 17.4)
on the given Markov decision process 'mdp' and a arbitarily small positive number 'epsilon'.
"""
function value_iteration(mdp::T; ϵ::Float64=0.001) where {T <: AbstractMarkovDecisionProcess}
    # initialise utility to zero for all states
    U_prime = Dict(collect(Pair(state, 0.0) for state in mdp.states))
    while true
        U = copy(U_prime)
        δ = 0.0
        for state in mdp.states
            U_prime[state] = reward(mdp, state) + mdp.gamma * maximum(map(x -> sum(p*U[newstate] for (p, newstate) in x),
                                                                          [transition_model(mdp, state, action) for action in actions(mdp, state)]))
            δ = max(δ, abs(U_prime[state] - U[state]));
        end
        if (δ < ((ϵ * (1 - mdp.gamma))/mdp.gamma))
            return U_prime;
        end
    end
end

function policy_extraction(mdp::T, U::Dict) where T<:AbstractMarkovDecisionProcess
	pi = Dict()
	for state in mdp.states
		A = collect(actions(mdp, state))
		expected_utility = map(x -> sum(p*U[newstate] for (p, newstate) in x),  [transition_model(mdp, state, action) for action in A])
		pi[state] = A[argmax(expected_utility)]
	end
	return pi
end

#function expected_utility(mdp::T, U::Dict, state::Tuple{Int64, Int64}, action::Nothing) where {T <: AbstractMarkovDecisionProcess}
#    return sum((p * U[state_prime] for (p, state_prime) in transition_model(mdp, state, action)));
#end

"""
    policy_evaluation(pi::Dict, U::Dict, mdp::T; k::Int64=20) where {T <: AbstractMarkovDecisionProcess}

Return the updated utilities of the MDP's states by applying the modified policy iteration
algorithm on the given Markov decision process 'mdp', utility function 'U', policy 'pi',
and number of Bellman updates to use 'k'.
"""
function policy_evaluation(pi::Dict, U::Dict, mdp::T; k::Int64=20) where {T <: AbstractMarkovDecisionProcess}
    for i in 1:k
        for state in mdp.states
            U[state] = reward(mdp, state) + mdp.gamma * sum(p * U[newstate] for (p, newstate) in transition_model(mdp, state, pi[state]));
        end
    end
    return U;
end

function policy_evaluation(pi::Dict, U::Dict, gmdp::GridMarkovDecisionProcess; k::Int64=20)
    for i in 1:k
        for state in gmdp.states
            U[state] = reward(gmdp, state) + gmdp.gamma * sum(p * U[newstate] for (p, newstate) in transition_model(gmdp, state, pi[state]));
        end
    end
    return U;
end

"""
    policy_iteration(mdp::T) where {T <: AbstractMarkovDecisionProcess}

Return a policy using the policy iteration algorithm (Fig. 17.7) given the Markov decision process 'mdp'.

it initialises the policy at random
"""
function policy_iteration(mdp::T) where {T <: AbstractMarkovDecisionProcess}
    U = Dict(Pair(state, 0.0) for state in mdp.states)
    pi = Dict(Pair(state, rand(actions(mdp, state))) for state in mdp.states)
    while true
        U = policy_evaluation(pi, U, mdp)
        unchanged = true
        for state in mdp.states
            A = collect(actions(mdp, state))
            expected_utility = map(x -> sum(p*U[newstate] for (p, newstate) in x),  [transition_model(mdp, state, action) for action in A])
            action = A[argmax(expected_utility)]
            if action != pi[state]
                pi[state] = action
                unchanged = false
            end
        end
        if unchanged
            return pi, U;
        end
    end
end


text:
using Meganet
using Base.Test

# include("../../src/optimization/sgd.jl")

# network
nf = 18
nc = 3
nclass=10
nex = 40
for TYPE=[Float64,Float32]
    K = getDenseKernel(TYPE,[nf,nf])
    Bout = randn(TYPE,nf,3)
    nLayer = getBatchNormLayer(TYPE,[div(nf,nc),nc],isTrainable=true)
    L     = getSingleLayer(TYPE,K,nLayer,Bout=Bout)
    net = getResNN(TYPE,L,8,1.0)

    # regularizers
    pRegTh = getTikhonovReg(TYPE,alpha=1e-5)
    pRegW = getTikhonovReg(TYPE,alpha=1e-5)
    # loss
    pLoss = getSoftMaxLoss(TYPE)
    # data
    objFun = dnnObjFctn(net,pLoss,pRegTh,pRegW)

    Y   = randn(TYPE,nf,nex)/10
    nw2 = (pLoss.addBias) ? nf+1 : nf
    C   = full(sparse(full(rand(1:nclass,nex)),collect(1:nex),ones(nex),nclass,nex))
    C  = convert.(TYPE,C)
    #initialize
    theta = initTheta(net)
    W     = randn(TYPE,nclass,nw2)/10

    opt = getSGDsolver(TYPE,learningRate=1e-4,maxEpochs=10,miniBatch=2)
    display(opt)

solve(opt,objFun::dnnObjFctn,[vec(theta);vec(W)],Y,C,Y,C)
end


text:
using ObjectPools
using Test, BenchmarkTools

@testset "Creating objects through the allocate" begin
    @test allocate(Int, 4) === 4
    @test allocate(Dict) isa Dict
    @test allocate(Dict{Int, Int}, [1 => 2])[1] === 2
end

@testset "Reusing objects" begin
    d1 = allocate(Dict{Int, Int})
    release(d1)

    d2 = allocate(Dict{Int, Int})
    @test d2 === d1

    d3 = allocate(Dict{Int, Int})
    @test d3 !== d1
end

function fill_pool(type, args...)
    ds = [allocate(type, args...) for i = 1:ObjectPools.MAX_POOL_SIZE]
    @test map(release, ds) == [true for i = 1:ObjectPools.MAX_POOL_SIZE]
end

const NUM_TYPES = 2

function getvalpools()
    for i = 1:NUM_TYPES
        ObjectPools.getpool(Val{i})
    end
end

function getvalpools_static()
    for i = 1:NUM_TYPES
        ObjectPools.getpool_static(Val{i})
    end
end

@testset "Reusing multiple objects of multiple types" begin
    fill_pool(Dict{String, Any})
    fill_pool(UInt, 1)
    fill_pool(Float16, 1.0)
    fill_pool(Float16, 1.0)
end

@testset "Filling up the pool" begin
    type = Dict{UInt, Int}
    extra_obj = allocate(type)
    ds = [allocate(type) for i = 1:ObjectPools.MAX_POOL_SIZE]
    @test map(release, ds) == [true for i = 1:ObjectPools.MAX_POOL_SIZE]
    @test release(extra_obj) == false
end

@testset "Benchmarks" begin
    @info "Running benchmarks:"
    @info "Dict constructor: Dict{Int, Int}()"
    @btime Dict{Int, Int}()

    @info "Creating new objects when the pool is empty: allocate(Dict{Int, Int})"
    @test (@btime allocate(Dict{Int, Int})) isa Dict{Int, Int}

    @info "Acquire-Release: begin a=allocate(Dict{Int, Int}); release(a) end"
    @btime begin a=allocate(Dict{Int, Int}); release(a) end

    @info "Acquire-Release from locally cached pool: begin pool=ObjectPools.getpool(Dict{Int, Int}); a=allocate(pool, Dict{Int, Int}); release(pool, a) end"
    @btime begin pool=ObjectPools.getpool(Dict{Int, Int}); a=allocate(pool, Dict{Int, Int}); release(pool, a) end

    pool = ObjectPools.getpool(Dict{Int, Int})
    @info "Acquire-Release from a constant pool: begin a=allocate(pool, Dict{Int, Int}); release(pool, a) end"
    @btime begin a=allocate($pool, Dict{Int, Int}); release($pool, a) end

    @info "getpool vs getpool_static (respectively) for $(NUM_TYPES)+ types and dynamic dispatch"
    getvalpools()
    getvalpools_static()
    @btime getvalpools()
    @btime getvalpools_static()
end


text:
# This file was generated, do not modify it. # hide
pipe.model.linear_regressor = RidgeRegressor()
fit!(pipe, rows=train)
ŷ = predict(pipe, rows=test)
round(rms(ŷ, y[test])^2, sigdigits=4)


text:
module CollegeStratMakie

using DocStringExtensions, FileIO, StatsBase
using CairoMakie, Makie, ColorSchemes
using FilesLH
using CollegeStratBase

# Helpers
export blank_plot, make_axis, make_legend, add_title!
# Settings
export default_theme, plot_defaults, get_colors
# Saving
export figsave
# Bar graphs
export bar_graph, bar_graph!
# Grouped bar graphs
export grouped_bar_graph, grouped_bar_graph!, grouped_bar_legend,
    grouped_bar_xyz
# Line graphs
export line_plot, add_line!, add_error_band!
# Scatter plots
export scatter_plot, add_scatter!
# Other graphs
export histogram_plot, contour_plot
export plot_regressions, plot_regression
# Subplots
export subplot_layout, subplot_pos, link_axes, make_axes!, subplots

include("helpers.jl");
include("settings.jl");
include("saving.jl");
include("bar_graphs.jl");
include("grouped_bar_graphs.jl");
include("line_graphs.jl");
include("scatter_graphs.jl");
include("subplots.jl");
include("regressions.jl");
include("xyz.jl");
include("other_plots.jl");

end # module


text:
using TcpInstruments
using Test

pwr = initialize(AgilentE36312A)
@info "Successfully connected $(pwr.model) at $(pwr.address)"

"""
Spec:
enable_output()
disable_output()
set_voltage(x)
get_voltage()
set_current_limit(x)
get_current_limit()
"""

@testset "Output" begin
    @info get_output(pwr)
    enable_output(pwr)

    @info get_output(pwr), "OUTPUT"
    disable_output(pwr)
    set_channel(pwr, 3)
    @test get_output(pwr) == false
    set_channel(pwr, 2)
    disable_output(pwr)
    @test get_output(pwr) == false
    set_channel(pwr, 1)
    disable_output(pwr)
    @test get_output(pwr) == false

    enable_output(pwr)

    @test get_output(pwr) == true
    set_channel(pwr, 2)
    @test get_output(pwr) == false
    set_channel(pwr, 3)
    @test get_output(pwr) == false

    set_channel(pwr, 2)
    enable_output(pwr)
    @test get_output(pwr) == true
end

@testset "Current" begin
    @info get_current_limit(pwr)
    set_current_limit(pwr, 1)
    @info get_current_limit(pwr)
    set_current_limit(pwr, 2)
    @info get_current_limit(pwr)
    set_current_limit(pwr, 3)
    @info get_current_limit(pwr)
    set_current_limit(pwr, 4)
    @info get_current_limit(pwr)
    set_current_limit(pwr, 0)
    @info get_current_limit(pwr)
    set_current_limit(pwr, "MAX")
    @info get_current_limit(pwr)
    set_current_limit(pwr, "MIN")
    @info get_current_limit(pwr)
    set_current_limit(pwr, 5)
    @info get_current_limit(pwr)
end

@testset "Voltage" begin
    function test_volt(v, i)
        set_voltage(pwr, v; chan=i)
        ans = get_voltage(pwr; chan=i)
        @info ans
        if ans == 6 && v != 6
            @info "This channel correctly limits to 6V"
        else
            @test ans == v
        end
    end
    for i in 1:3
        @info "CHANNEL: $i"
        test_volt(2, i)
        test_volt(4, i)
        test_volt(6, i)

        test_volt(8, i)
        test_volt(20, i)
        test_volt(25, i)
        @info get_voltage(pwr; chan=i)
        set_voltage(pwr, 0; chan=i)
        @info get_voltage(pwr; chan=i)
        set_voltage(pwr, "MAX"; chan=i)
        @info get_voltage(pwr; chan=i)
        set_voltage(pwr, "MIN"; chan=i)
        @info get_voltage(pwr; chan=i)
    end
end

@info query(pwr, "INST:NSEL?")
@info query(pwr, "APPLy? CH1"), 1
@info query(pwr, "APPLy? CH2"), 2
@info query(pwr, "APPLy? CH3"), 3

set_channel(pwr, 3)
set_voltage(pwr, 7.7; chan=2)
@test get_channel(pwr) == "3"
set_voltage(pwr, 7.7)
@info get_voltage(pwr)
@info get_channel(pwr)

terminate(pwr)

@info "Successfully disconnected"
@info "Goodbye"


text:
"""
This module supplies some functions for testing of the implementation.
Look at `?` for

* `fdtest` : general finite-difference test
* `fdtest_R2R`: fd test for F : ℝ → ℝ
"""
module Testing

using Test
using JuLIP: AbstractCalculator, AbstractAtoms, energy, gradient, forces,
         calculator, set_positions!, dofs,
         mat, vecs, positions, rattle!, set_dofs!, set_calculator!
using JuLIP.Potentials: PairPotential, evaluate, evaluate_d, @D
using Printf
using LinearAlgebra: norm

import ACEbase.Testing: print_tf, test_fio, h0, h1, h2, h3

export fdtest, fdtest_hessian, h0, h1, h2, h3, print_tf, test_fio

"""
first-order finite-difference test for scalar F

* `fdtest(F::Function, dF::Function, x)`
* `fdtest(V::PairPotential, r::Vector{Float64})`
* `fdtest(calc::AbstractCalculator, at::AbstractAtoms)`
"""
function fdtest(F::Function, dF::Function, x; verbose=true)
   errors = Float64[]
   E = F(x)
   dE = dF(x)
   # loop through finite-difference step-lengths
   verbose && @printf("---------|----------- \n")
   verbose && @printf("    h    | error \n")
   verbose && @printf("---------|----------- \n")
   for p = 2:11
      h = 0.1^p
      dEh = copy(dE)
      for n = 1:length(dE)
         x[n] += h
         dEh[n] = (F(x) - E) / h
         x[n] -= h
      end
      push!(errors, norm(dE - dEh, Inf))
      verbose && @printf(" %1.1e | %4.2e  \n", h, errors[end])
   end
   verbose && @printf("---------|----------- \n")
   if minimum(errors) <= 1e-3 * maximum(errors)
      verbose && println("passed")
      return true
   else
      @warn("""It seems the finite-difference test has failed, which indicates
      that there is an inconsistency between the function and gradient
      evaluation. Please double-check this manually / visually. (It is
      also possible that the function being tested is poorly scaled.)""")
      return false
   end
end

function fdtest_hessian(F::Function, dF::Function, x; verbose=true)
   errors = Float64[]
   F0 = F(x)
   dF0 = dF(x)
   dFh = copy(Matrix(dF0))
   @assert size(dFh) == (length(F0), length(x))
   # loop through finite-difference step-lengths
   verbose &&  @printf("---------|----------- \n")
   verbose &&  @printf("    h    | error \n")
   verbose &&  @printf("---------|----------- \n")
   for p = 2:11
      h = 0.1^p
      for n = 1:length(x)
         x[n] += h
         dFh[:, n] = (F(x) - F0) / h
         x[n] -= h
      end
      push!(errors, norm(dFh - dF0, Inf))
      verbose &&  @printf(" %1.1e | %4.2e  \n", h, errors[end])
   end
   @printf("---------|----------- \n")
   if minimum(errors) <= 1e-3 * maximum(errors)
      verbose &&  println("passed")
      return true
   else
      @warn("""It seems the finite-difference test has failed, which indicates
            that there is an inconsistency between the function and gradient
            evaluation. Please double-check this manually / visually. (It is
            also possible that the function being tested is poorly scaled.)""")
      return false
   end
end

"finite-difference test for a function V : ℝ → ℝ"
function fdtest_R2R(F::Function, dF::Function, x::Vector{Float64};
                     verbose=true)
   errors = Float64[]
   E = [ F(t) for t in x ]
   dE = [ dF(t) for t in x ]
   # loop through finite-difference step-lengths
   if verbose
      @printf("---------|----------- \n")
      @printf("    h    | error \n")
      @printf("---------|----------- \n")
   end
   for p = 2:11
      h = 0.1^p
      dEh = ([F(t+h) for t in x ] - E) / h
      push!(errors, norm(dE - dEh, Inf))
      if verbose
         @printf(" %1.1e | %4.2e  \n", h, errors[end])
      end
   end
   if verbose
      @printf("---------|----------- \n")
   end
   if minimum(errors) <= 1e-3 * maximum(errors[1:2])
      println("passed")
      return true
   else
      @warn("""is seems the finite-difference test has failed, which indicates
            that there is an inconsistency between the function and gradient
            evaluation. Please double-check this manually / visually. (It is
            also possible that the function being tested is poorly scaled.)""")
      return false
   end
end

fdtest(V::PairPotential, r::AbstractVector; kwargs...) =
               fdtest_R2R(s -> V(s), s -> (@D V(s)), collect(r); kwargs...)

function fdtest(calc::AbstractCalculator, at::AbstractAtoms;
                verbose=true, rattle=0.01)
   X0 = copy(positions(at))
   calc0 = calculator(at)
   set_calculator!(at, calc)
   # random perturbation to positions (and cell for VariableCell)
   # perturb atom positions a bit to get out of equilibrium states
   # Don't use `rattle!` here which screws up the `VariableCell` constraint
   # test!!! (but why?!?!?)
   x = dofs(at)
   x += rattle * rand(length(x))
   # call the actual FD test
   result = fdtest( x -> energy(at, x), x -> gradient(at, x), x;
                    verbose = verbose )
   # restore original atom positions
   set_positions!(at, X0)
   set_calculator!(at, calc0)
   return result
end

end


text:
module GHWT

include("utils.jl")

using ..GraphSignal, ..GraphPartition, ..BasisSpecification, LinearAlgebra

include("common.jl")

export ghwt_core!, ghwt_analysis!, fine2coarse!, ghwt_synthesis, ghwt_c2f_bestbasis, ghwt_f2c_bestbasis, ghwt_bestbasis, GHWT_jkl

"""
    ghwt_core!(GP::GraphPart)

performs the forward GHWT transform, which computes tag and compinfo of GP, but not the expansion coefficients.

### Input Arguments
* `GP::GraphPart`: a GraphPart object (with or without tag and compinfo data); note that tag and compinfo will be modified after this function is executed
"""
function ghwt_core!(GP::GraphPart)

    #
    # 0. Preliminaries
    #
    rs = GP.rs                  # for simple notation!
    (N, jmax) = Base.size(rs)
    N = N - 1

    # tag -- when expressed in binary, tag indicates the sequence of
    # average (0's) and difference (1's) operations used to obtain the given
    # basis vector
    if isempty(GP.tag)
        GP.tag = zeros(Int, N, jmax)
    end
    tag = GP.tag                # for simple notation!

    # allocate space for compinfo (used for synthesis from f2c dictionary)
    # if tag == 0 && n == 1:                        compinfo = 0
    # if tag == 0 && n >= 2:                        compinfo = n1
    # if tag == 1:                                  compinfo = n2
    # if tag >= 2 && coeff is formed from 2 coeffs: compinfo = 1
    # if tag >= 2 && coeff is formed from 1 coeff:  compinfo = 0
    if isempty(GP.compinfo)
        GP.compinfo = zeros(Int, N, jmax)
    end
    compinfo = GP.compinfo      # for simple notation!

    #
    # 1. Perform the transform
    #
    for j = (jmax - 1):-1:1 # # Bottom-1 (jmax-1) to up to the coarsest (j=1)
        regioncount = count(!iszero,rs[:, j]) - 1
        for r = 1:regioncount
            rs1 = rs[r, j]      # the start of the 1st subregion
            rs3 = rs[r + 1, j]  # 1 + the end of the 2nd subregion
            n = rs3 - rs1       # # of points in the current region
            if n > 1            # single node region passes this part
                rs2 = rs1 + 1   # the start of the 2nd subregion
                while rs2 < rs3 && tag[rs2, j + 1] != 0 ### && rs2 < N+1
                    rs2 += 1
                end
                if rs2 == rs3 # the parent region is a copy of the subregion
                    tag[rs1:(rs3 - 1), j] = tag[rs1:(rs3 - 1), j + 1]
                else           # the parent region has 2 child regions
                    n1 = rs2 - rs1 # # of pts in the 1st subregion
                    n2 = rs3 - rs2 # # of pts in the 2nd subregion
                    # SCALING COEFFICIENT (n > 1)
                    compinfo[rs1, j] = n1;
                    # HAAR COEFFICIENT
                    compinfo[rs1 + 1, j] = n2; tag[rs1 + 1, j] = 1
                    # WALSH COEFFICIENTS
                    # sweep through the coefficients in subregions 1 & 2
                    parent = rs1 + 2 # ind of the new coeffs to be created on level j
                    child1 = rs1 + 1 # ind of the current coeffs in subregion 1
                    child2 = rs2 + 1 # ind of the current coeffs in subregion 2
                    while parent < rs3
                        # no matching coefficient (use subregion 1)
                        if child1 < rs2 && (child2 == rs3 || tag[child1, j + 1] < tag[child2, j + 1])
                            tag[parent, j] = 2 * tag[child1, j + 1]
                            child1 += 1; parent += 1
                        # no matching coefficient (use subregion 2)
                        elseif child2 < rs3 && (child1 == rs2 || tag[child2, j + 1] < tag[child1, j + 1])
                            tag[parent, j] = 2 * tag[child2, j + 1]
                            child2 += 1; parent += 1
                        # matching coefficients
                        else
                            tag[parent, j] = 2 * tag[child1, j + 1]
                            #tag[parent + 1, j] = 2 * tag[child1, j + 1] + 1
                            tag[parent + 1, j] = tag[parent, j] + 1
                            compinfo[parent, j] = 1; compinfo[parent + 1, j]  = 1
                            child1 += 1; child2 += 1; parent += 2
                        end # of if child1 < rs2 ...
                    end # of while
                end # of if rs2 == rs3 ...
            elseif n <= 0 # n must be positive
                error("n must be positive: n = ", n)
            end # of if n > 1 ... elseif n <= 0 ...; note: do nothing for n==1
        end # of for r = 1:regioncount
    end # of for j = (jmax - 1):-1:1
end # of function ghwt_core!

"""
    ghwt_core!(GP::GraphPart, dmatrix::Array{Float64,3})

performs the forward GHWT transform, which generates expansion coefficients, tag, and compinfo.

### Input Arguments
* `GP::GraphPart`: a GraphPart object (with or without tag and compinfo data); note that tag and compinfo will be modified after this function is executed
* `dmatrix::Array{Float64,3}`: a matrix of expansion coefficients with only the last column filled in as the original input signal(s) `f`; after this function is executed, this matrix contains whole expansion coefficients of the input signal(s) relative to the GHWT dictionary
"""
function ghwt_core!(GP::GraphPart, dmatrix::Array{Float64,3})

    #
    # 0. Preliminaries
    #
    rs = GP.rs                  # for simple notation!
    (N, jmax) = Base.size(rs)
    N = N - 1

    # tag -- when expressed in binary, tag indicates the sequence of
    # average (0's) and difference (1's) operations used to obtain the given
    # basis vector
    if isempty(GP.tag)
        GP.tag = zeros(Int, N, jmax)
    end
    tag = GP.tag                # for simple notation!

    # allocate space for compinfo (used for synthesis from fine-to-coarse dictionary)
    # if tag == 0 && n == 1:                        compinfo = 0
    # if tag == 0 && n >= 2:                        compinfo = n1
    # if tag == 1:                                  compinfo = n2
    # if tag >= 2 && coeff is formed from 2 coeffs: compinfo = 1
    # if tag >= 2 && coeff is formed from 1 coeff:  compinfo = 0
    if isempty(GP.compinfo)
        GP.compinfo = zeros(Int, N, jmax)
    end
    compinfo = GP.compinfo      # for simple notation!

    #
    # 1. Perform the transform
    #
    for j = (jmax - 1):-1:1 # Bottom-1 (jmax-1) to up to the coarsest (j=1)
        regioncount = count(!iszero,rs[:, j]) - 1
        for r = 1:regioncount
            rs1 = rs[r, j]      # the start of the 1st subregion
            rs3 = rs[r + 1, j]  # 1 + the end of the 2nd subregion
            n = rs3 - rs1       # # of points in the current region
            if n == 1           # single node region
                # SCALING COEFFICIENT (n == 1)
                dmatrix[rs1, j, :] = dmatrix[rs1, j + 1, :]
            elseif n > 1
                rs2 = rs1 + 1 # the start of the 2nd subregion
                while rs2 < rs3 && tag[rs2, j + 1] != 0 ### && rs2 < N+1
                    rs2 += 1
                end
                if rs2 == rs3 # the parent region is a copy of the subregion
                    dmatrix[rs1:(rs3 - 1), j, :] = dmatrix[rs1:(rs3 - 1), j + 1, :]
                    tag[rs1:(rs3 - 1), j] = tag[rs1:(rs3 - 1), j + 1]
                else       # the parent region has 2 child regions
                    n1 = rs2 - rs1 # # of pts in the 1st subregion
                    n2 = rs3 - rs2 # # of pts in the 2nd subregion

                    # SCALING COEFFICIENT (n > 1)
                    dmatrix[rs1, j, :] =
                        ( sqrt(n1) * dmatrix[rs1, j + 1, :] +
                          sqrt(n2) * dmatrix[rs2, j + 1, :] ) / sqrt(n)
                    compinfo[rs1, j] = n1

                    # HAAR COEFFICIENT
                    dmatrix[rs1 + 1, j, :] =
                        ( sqrt(n2) * dmatrix[rs1, j + 1, :] -
                          sqrt(n1) * dmatrix[rs2, j + 1, :] ) / sqrt(n)
                    # Note these sqrt's arguments are not typos!!
                    compinfo[rs1 + 1, j] = n2; tag[rs1 + 1, j] = 1

                    # WALSH COEFFICIENTS
                    # sweep through the coefficients in subregions 1 & 2
                    parent = rs1 + 2 # ind of the new coeffs to be created on level j
                    child1 = rs1 + 1 # ind of the current coeffs in subregion 1
                    child2 = rs2 + 1 # ind of the current coeffs in subregion 2
                    while parent < rs3
                        # no matching coefficient (use subregion 1)
                        if child1 < rs2 && (child2 == rs3 || tag[child1, j + 1] < tag[child2, j + 1])
                            dmatrix[parent, j, :] = dmatrix[child1, j + 1, :]
                            tag[parent, j] = 2 * tag[child1, j + 1]
                            child1 += 1; parent += 1
                        # no matching coefficient (use subregion 2)
                        elseif child2 < rs3 && (child1 == rs2 || tag[child2, j + 1] < tag[child1, j + 1])
                            dmatrix[parent, j, :] = dmatrix[child2, j + 1, :]
                            tag[parent, j] = 2 * tag[child2, j + 1]
                            child2 += 1; parent += 1
                        # matching coefficients
                        else
                            dmatrix[parent, j, :] =
                                ( dmatrix[child1, j + 1, :] +
                                  dmatrix[child2, j + 1, :] ) / sqrt(2)
                            tag[parent, j] = 2 * tag[child1, j + 1]
                            compinfo[parent, j] = 1

                            dmatrix[parent + 1, j, :] =
                                ( dmatrix[child1, j + 1, :] -
                                  dmatrix[child2, j + 1, :] ) / sqrt(2)
                            #tag[parent + 1, j] = 2 * tag[child1, j + 1] + 1
                            tag[parent + 1, j] = tag[parent, j] + 1
                            compinfo[parent + 1, j]  = 1
                            child1 += 1; child2 += 1; parent += 2
                        end # of if child1 < rs2 ...
                    end # of while
                end # of if rs2 == rs3 ...
            else # n must be positive
                error("n must be positive: n = ", n)
            end # of if n == 1 ... elseif n > 1 ...
        end # of for r = 1:regioncount
    end # of for j = (jmax - 1):-1:1
end # of function ghwt_core!

"""
    dmatrix = ghwt_analysis!(G::GraphSig, GP::GraphPart = nothing, c2f::Bool = true)

For a GraphSig object `G`, generate the matrix of GHWT expansion coefficients

### Input Arguments
* `G::GraphSig`: an input GraphSig object
* `GP::GraphPart`: an input GraphPart object (optional); after this function is run, `GP`'s `compinfo`, `tag`, etc. are filled

### Output Argument
* `dmatrix::Array{Float64,3}`: the 3D array of expansion coefficients (i.e., for each input signal vector, the matrix of coefficients; hence, for multiple input signals, the coefficients are organized as a 3D array)
"""
function ghwt_analysis!(G::GraphSig; GP::GraphPart = nothing, c2f::Bool = true)

    #
    # 0. Preliminaries
    #
    if isnothing(GP)
        GP = partition_tree_fiedler(G)
    end

    # Simplify the subsequent notation without overhead. We can do this
    # since `ind`, `rs`, `f` share the same memory as `GP.ind`, `GP.rs`, `G.f`
    ind = GP.ind
    rs = GP.rs
    f = G.f
    (N, jmax) = Base.size(rs)
    N = N - 1

    # allocate space for the expansion coefficients
    fcols = size(f, 2)
    dmatrix = zeros(N, jmax, fcols)
    dmatrix[:, jmax, :] = f[ind, :]

    ## 1. Perform the transform

    # generate expansion coefficients, tag, and compinfo
    ghwt_core!(GP, dmatrix)

    # fill in the remaining (i.e. fine-to-coarse) GraphPart fields
    if c2f
        return dmatrix
    else
        return fine2coarse!(GP, dmatrix=dmatrix, coefp = true)
    end

end # of ghwt_analysis!

"""
    (dmatrixf2c, IX) = fine2coarse!(GP::GraphPart;
    dmatrix::Array{Float64,3} = zeros(0, 0, 0),
    coefp::Bool = false, indp::Bool = false)

Fill in the fine-to-coarse info (rs2f2c, tagf2c, and compinfof2c) in a
GraphPart object.  Also, rearrange a matrix of expansion coefficients.

### Input Arguments
* `GP::GraphPart`: an input GraphPart object without fine-to-coarse info (rsf2c, tagf2c, compinfof2c); after this function, rsf2c, tagf2c, compinfof2c are filled. Note that rs, tag, compinfo are intact.
* `dmatrix::Array{Float64,3}`: a matrix of expansion coefficients in coarse-to-fine arrangement (default: null matrix)
* `coefp::Bool`: a flag to return the rearranged f2c coefficients (default: false)
* `indp::Bool`: a flag to return the reordering index (default: false)

### Output Arguments
* `dmatrixf2c::Array{Float64,3}`: a matrix of expansion coefficients in fine-to-coarse arrangement (if requested)
* `IX::Vector{Unsigned}`: the reordering index for all levels
"""
function fine2coarse!(GP::GraphPart;
                      dmatrix::Array{Float64,3} = zeros(0, 0, 0),
                      coefp::Bool = false, indp::Bool = false)
    #
    # 0. Preliminaries
    #

    # get constants
    (N, jmax) = Base.size(GP.rs)
    N = N - 1

    # allocate spaces
    if isempty(GP.rsf2c)
        GP.rsf2c = zeros(Int, N + 1, jmax)
    end
    GP.rsf2c[1, :] .= 1
    GP.rsf2c[2, 1] = N + 1
    if isempty(GP.tagf2c)
        GP.tagf2c = zeros(Int, N, jmax)
    end
    if isempty(GP.compinfof2c)
        GP.compinfof2c = zeros(Int, N, jmax)
    end
    IX = zeros(Int, N, jmax)
    if coefp
        dmatrixf2c = zeros(Base.size(dmatrix))
    end

    # make sure the coarse-to-fine tag field is filled in
    if isempty(GP.tag) || isempty(GP.compinfo)
        ghwt_core!(GP)
    end

    # for simple notation!
    ind = GP.ind; rs = GP. rs
    tag = GP.tag; compinfo = GP.compinfo
    rsf2c = GP.rsf2c; tagf2c = GP.tagf2c; compinfof2c = GP.compinfof2c

    #
    # 1. Generate the fine-to-coarse dictionary
    #

    # put the coefficients into fine-to-coarse order
    for j = 1:jmax
        # put the basis into tag order
        # MATLAB: [GP.tagf2c(:,j),IX(:,j)] = sort(GP.tag(:,jmax+1-j));
        IX[:, j] = sortperm(tag[:, jmax + 1 - j])
        tagf2c[:, j] = tag[IX[:, j], jmax + 1 - j]
        compinfof2c[:, j] = compinfo[IX[:, j], jmax + 1 - j]

        if coefp
            dmatrixf2c[:, j, :] = dmatrix[IX[:, j], jmax + 1 - j, :]
        end

        # fill in the fine-to-coarse regionstarts
        if j > 1
            rsf2crow = 2
            for row = 2:N
                if tagf2c[row, j] > tagf2c[row - 1, j]
                    rsf2c[rsf2crow, j] = row
                    rsf2crow += 1
                end
            end
            rsf2c[rsf2crow, j] = N + 1
        end
    end # of for j = 1:jmax

    #
    # 2. Return the appropriate variables if requested
    #
    if coefp && indp
        return dmatrixf2c, IX
    elseif coefp && !indp
        return dmatrixf2c
    elseif !coefp && indp
        return IX
    end
end # of function fine2coarse!

"""
    f = ghwt_synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec)

Given a vector of GHWT expansion coefficients and info about the graph
partitioning and the choice of basis, reconstruct the signal

### Input Arguments
* `dvec::Matrix{Float64}`: the expansion coefficients corresponding to the chosen basis
* `GP::GraphPart`: an input GraphPart object
* `BS::BasisSpec`: an input  BasisSpec object

### Output Arguments
* `f::Matrix{Float64}`: the reconstructed signal(s)
"""
function ghwt_synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec)

    #
    # 0. Preliminaries
    #
    # if necessary, fill in the tag info in GP
    if isempty(GP.tag)
        ghwt_core!(GP)
    end

    # figure out which dictionary is used: coarse-to-fine or fine-to-coarse
    if !BS.c2f && isempty(GP.tagf2c)
        fine2coarse!(GP)
    end

    # constants
    N = Base.length(GP.ind)
    jmax = Base.size(GP.rs, 2)
    fcols = Base.size(dvec, 2)

    # fill in the appropriate entries of dmatrix
    dmatrix = dvec2dmatrix(dvec, GP, BS)

    #
    # 1. Synthesis from the basis dictionary
    #
    if BS.c2f                   # for coarse-to-fine case
        # for simpler notation!
        dvec_loc = zeros(Int, size(GP.tag))
        for i = 1:length(BS.levlist)
            dvec_loc[BS.levlist[i][1], BS.levlist[i][2]] = 1
        end

        rs = GP.rs
        tag = GP.tag
        for j = 1:(jmax - 1)    # from top to bottom-1
            regioncount = count(!iszero, rs[:, j]) - 1
            for r = 1:regioncount
                rs1 = rs[r, j]      # the start of the 1st subregion
                rs3 = rs[r + 1, j]  # 1 + the end of the 2nd subregion
                n = rs3 - rs1       # # of points in the current region
                # only proceed forward if coefficients do not exist
                #if count(!iszero, dmatrix[rs1:(rs3 - 1), j + 1, :]) == 0 &&
                #    count(!iszero, dmatrix[rs1:(rs3 - 1), j, :]) > 0
                if count(!iszero, dvec_loc[rs1:(rs3-1),j]) > 0
                    if n == 1   # single node region
                        # SCALING COEFFICIENT (n == 1)
                        if dvec_loc[rs1,j] == 1
                            dmatrix[rs1, j + 1, :] = dmatrix[rs1, j, :]
                            dvec_loc[rs1,j+1] = 1
                        end
                    elseif n > 1
                        rs2 = rs1 + 1 # the start of the 2nd subregion
                        while rs2 < rs3 && tag[rs2, j + 1] != 0 ### && rs2 < N+1
                            rs2 += 1
                        end
                        if rs2 == rs3 # the parent is a copy of the subregion
                            if dvec_loc[rs1:rs3-1,j] == 1
                                dmatrix[rs1:(rs3 - 1), j + 1, :] = dmatrix[rs1:(rs3 - 1), j, :]
                                dvec_loc[rs1:rs3-1,j+1] = 1
                            end
                        else   # the parent region has 2 child regions
                            n1 = rs2 - rs1 # # of pts in the 1st subregion
                            n2 = rs3 - rs2 # # of pts in the 2nd subregion

                            # SCALING COEFFICIENTS (n > 1)
                            if dvec_loc[rs1,j] == 1 && dvec_loc[rs1+1,j] == 1
                                dmatrix[rs1, j + 1, :] =
                                    ( sqrt(n1) * dmatrix[rs1, j, :] +
                                      sqrt(n2) * dmatrix[rs1 + 1, j, :] ) / sqrt(n)

                            # HAAR COEFFICIENTS
                                dmatrix[rs2, j + 1, :] =
                                    ( sqrt(n2) * dmatrix[rs1, j, :] -
                                      sqrt(n1) * dmatrix[rs1 + 1, j, :] ) / sqrt(n)

                                dvec_loc[rs1,j+1] = 1
                                dvec_loc[rs2,j+1] = 1
                            end

                            # WALSH COEFFICIENTS
                            # search through the remaining coefs in each subregion
                            parent = rs1 + 2
                            child1 = rs1 + 1
                            child2 = rs2 + 1
                            while child1 < rs2 || child2 < rs3
                                # subregion 1 has the smaller tag
                                if child2 == rs3 ||
                                    (tag[child1, j + 1] < tag[child2, j + 1] &&
                                     child1 < rs2)

                                    if dvec_loc[parent,j]==1
                                        dmatrix[child1, j + 1, :] =
                                            dmatrix[parent, j, :]
                                        dvec_loc[child1, j+1] =1
                                    end
                                    child1 += 1; parent += 1

                                # subregion 2 has the smaller tag
                                elseif child1 == rs2 ||
                                    (tag[child2, j + 1] < tag[child1, j + 1] &&
                                     child2 < rs3)

                                    if dvec_loc[parent, j] == 1
                                        dmatrix[child2, j + 1, :] =
                                            dmatrix[parent, j, :]
                                        dvec_loc[child2, j+1] = 1
                                    end
                                    child2 += 1; parent += 1

                                # both subregions have the same tag
                                else
                                    if dvec_loc[parent,j] == 1 && dvec_loc[parent+1, j] == 1
                                        dmatrix[child1, j + 1, :] =
                                            ( dmatrix[parent, j, :] +
                                              dmatrix[parent + 1, j, :] ) / sqrt(2)
                                        dmatrix[child2, j + 1, :] =
                                            ( dmatrix[parent, j, :] -
                                              dmatrix[parent + 1, j, :] ) / sqrt(2)
                                        dvec_loc[child1,j+1]=1
                                        dvec_loc[child2,j+1]=1
                                    end
                                    child1 += 1; child2 += 1; parent += 2
                                end # of if child2 == r3 ...
                            end # of while child1 < rs2 ...
                        end # of if rs2 == rs3 ... else
                    end # of if n == 1 ... elseif n > 1 ...
                end # of if countnz(...)
            end # of for r = 1:regioncount
        end # of for j = 1:(jmax - 1)
        # MATLAB: if fcols == 1
        #            f = dmatrix[:, end]
        #         else
        #            f = reshape(dmatrix[:, end, :], N, fcols)
        #         end
        # GS.f = dmatrix[:, end, :] # Do this later for clarity
    else                        # for fine-to-coarse case
        # for simpler notation!
        rsf2c = GP.rsf2c
        tagf2c = GP.tagf2c
        compinfof2c = GP.compinfof2c
        for j = (jmax - 1):-1:1 # from bottom-1 to top
            regioncount = count(!iszero, rsf2c[:, j]) - 1
            for r = 1:regioncount
                rs1 = rsf2c[r, j] # the start of the 1st subregin
                rs3 = rsf2c[r + 1, j] # 1 + the end of the 2nd subregion
                # only proceed forward if coefficients do not exist
                if count(!iszero, dmatrix[rs1:(rs3 - 1), j, :]) == 0 &&
                    count(!iszero, dmatrix[rs1:(rs3 - 1), j + 1, :]) > 0
                    # one subregion ==> copy the coefficients
                    if tagf2c[rs1, j + 1] == tagf2c[rs3 - 1, j + 1]
                        dmatrix[rs1:(rs3 - 1), j, :] =
                            dmatrix[rs1:(rs3 - 1), j + 1, :]
                    # two subregions ==> compute the coefficients
                    else
                        rs2 = rs1 + 1 # the start of the 2nd subregion
                        while rs2 < rs3 &&
                            tagf2c[rs1, j + 1] == tagf2c[rs2, j + 1] ### && rs2 < N+1
                            rs2 += 1
                        end
                        if rs2 == rs3 # the parent is a copy of the subregion
                            # no need to do the following
                            # dmatrix[rs1:(rs3 - 1), (j + 1) , :] =
                            #  dmatrix[rs1:(rs3 - 1) , j, :]
                            nothing
                        else   # the parent region has 2 child regions
                            parent = rs1 # the current coef in the parent
                            child1 = rs1 # the current coefs in the children
                            child2 = rs2

                            # SCALING & HAAR COEFFICIENTS
                            if tagf2c[rs1, j] == 0
                                while parent < rs3
                                    # 1 coef formed from 1 coef
                                    if compinfof2c[child1, j + 1] == 0
                                        dmatrix[parent, j, :] =
                                            dmatrix[child1, j + 1, :]
                                        parent += 1; child1 += 1
                                    # 2 coefs formed from 2 coefs
                                    else
                                        n1 = compinfof2c[child1, j + 1]
                                        n2 = compinfof2c[child2, j + 1]
                                        n = n1 + n2 # # of points in the parent
                                        # SCALING COEFFICIENT
                                        dmatrix[parent, j, :] =
                                            ( sqrt(n1) * dmatrix[child1, j + 1, :] +
                                              sqrt(n2) * dmatrix[child2, j + 1, :] ) / sqrt(n)
                                        # HAAR COEFFICIENT
                                        dmatrix[parent + 1, j, :] =
                                            ( sqrt(n2) * dmatrix[child1, j + 1, :] -
                                              sqrt(n1) * dmatrix[child2, j + 1, :] ) / sqrt(n)
                                        parent += 2; child1 += 1; child2 += 1
                                    end
                                end # of while
                            # WALSH COEFFICIENTS
                            else
                                while parent < rs3
                                    # 1 coef formed from 1 coef
                                    if compinfof2c[child1, j + 1] == 0
                                        dmatrix[parent, j, :] =
                                            dmatrix[child1, j + 1, :]
                                        parent += 1; child1 += 1

                                    # 2 coefs formed from 2 coefs
                                    else
                                        dmatrix[parent, j, :] =
                                            ( dmatrix[child1, j + 1, :] +
                                              dmatrix[child2, j + 1, :] ) / sqrt(2)
                                        dmatrix[parent + 1, j, :] =
                                            ( dmatrix[child1, j + 1, :] -
                                              dmatrix[child2, j + 1, :] ) / sqrt(2)
                                        parent += 2; child1 += 1; child2 += 1
                                    end # of if compinfof2c[...]
                                end # of while parent < rs3
                            end # of if tagf2c[rs1, j] == 0 ... else ...
                        end # of if rs2 == rs3 ... else ...
                    end # of if tagf2c[rs1, j + 1] == tagf2c[rs3 - 1, j + 1] ...
                end # of if countnz(dmatrix[rs1:(rs3 - 1), j, :]) == 0 ...
            end # for r = 1:regioncount
        end # for j=(jmax - 1):-1:1
        # GS.f = dmatrix[:, 1, :]
    end # of if BS.c2f ... else ...

    #
    # 2. Final prepartion for returning variables
    #

    ftemp = dmatrix[:, 1, :] # for fine-to-coarse and make f available
    if BS.c2f
        ftemp = dmatrix[:, end, :] # for coarse-to-fine
    end
    # put the reconstructed values in the correct order and return it.
    f = zeros(size(ftemp))
    f[GP.ind, :] = ftemp

    return f
end # of function ghwt_synthesis (without G input)

"""
    (f, GS) = ghwt_synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec, G::GraphSig)

Given a vector of GHWT expansion coefficients and info about the graph
partitioning and the choice of basis, reconstruct the signal

### Input Arguments
* `dvec::Matrix{Float64}`: the expansion coefficients corresponding to the chosen basis
* `GP::GraphPart`: an input GraphPart object
* `BS::BasisSpec`: an input  BasisSpec object
* `G::GraphSig`: an input  GraphSig object

### Output Arguments
* `f::Matrix{Float64}`: the reconstructed signal(s)
* `GS::GraphSig`: the reconstructed GraphSig object
"""
function ghwt_synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec, G::GraphSig)

    # call the normal ghwt_synthesis
    f = ghwt_synthesis(dvec, GP, BS)

    # create a GraphSig object with the reconstructed data and return it with f.
    GS = deepcopy(G)        # generate a new copy of G as GS
    replace_data!(GS, f)    # then replace the graph signal field
    return f, GS
end # of ghwt_synthesis (with G input)

"""
    (dvecc2f, BSc2f) = ghwt_c2f_bestbasis(dmatrix::Array{Float64,3}, GP::GraphPart;
                                          cfspec::Any = 1.0, flatten::Any = 1.0,
                                          j_start::Int = 1, j_end::Int = size(dmatrix,2),
                                          useParent::Bool = true)

Select the coarse-to-fine best basis from the matrix of GHWT expansion coefficients

### Input Arguments
* `dmatrix::Array{Float64,3}`: the matrix of expansion coefficients
* `GP::GraphPart`: an input GraphPart object
* `cfspec::Any`: the specification of cost functional to be used (default = 1.0, i.e., 1-norm)
* `flatten::Any`: the method for flattening vector-valued data to scalar-valued data (default = 1.0, i.e, 1-norm)
* `useParent::Bool`: the flag to indicate if we update the selected best basis
    subspace to the parent when parent and child have the same cost (default = true)

###  Output Arguments
* `dvecc2f::Matrix{Float64}`: the vector of expansion coefficients corresponding to the coarse-to-fine best basis
* `BSc2f::BasisSpec`: a BasisSpec object which specifies the coarse-to-fine best basis
"""
function ghwt_c2f_bestbasis(dmatrix::Array{Float64,3}, GP::GraphPart;
                            cfspec::Any = 1.0, flatten::Any = 1.0,
                            j_start::Int = 1, j_end::Int = size(dmatrix,2),
                            useParent::Bool = true)

    # determine the cost functional to be used
    costfun = cost_functional(cfspec)

    # constants and dmatrix cleanup
    (N, jmax, fcols) = Base.size(dmatrix)
    if ~(1 <= j_start <= j_end <= jmax)
        error("j_start and j_end should satisfy the inequality 1 <= j_start <= j_end <= "*string(jmax))
    end
    dmatrix[ abs.(dmatrix) .< 10^2 * eps() ] .= 0

    # "flatten" dmatrix
    if fcols > 1
        dmatrix0 = deepcopy(dmatrix)      # keep the original dmatrix as dmatrix0
        if !isnothing(flatten)
            dmatrix = dmatrix_flatten(dmatrix, flatten)
        end
    end

    ## Find the best-basis from the GHWT coarse-to-fine dictionary
    # allocate/initialize
    dvecc2f = dmatrix[:, j_end, 1]
    levlistc2f = j_end * ones(Int, N)

    # set the tolerance with sign
    tol = 10^4 * eps() * (1 - useParent * 2)

    # perform the basis search
    for j = (j_end - 1):-1:j_start
        regioncount = count(!iszero, GP.rs[:, j]) - 1
        for r = 1:regioncount
            indr = GP.rs[r, j]:(GP.rs[r + 1, j]-1)
            ##### compute the cost of the current best basis
            costBB = costfun(dvecc2f[indr])
            costNEW = costfun(dmatrix[indr, j, 1])
            if costBB >= costNEW + tol
                #(dvecc2f[indr], levlistc2f[indr]) = bbchange(dmatrix[indr, j, 1], j)
                dvecc2f[indr] = dmatrix[indr, j, 1]
                levlistc2f[indr] .= j
            end
        end
    end
    #levlistc2f = Array{Int}(levlistc2f[ levlistc2f .!= 0 ])
    levlistc2f = collect(enumerate(levlistc2f))
    BSc2f = BasisSpec(levlistc2f, c2f = true, description = "GHWT c2f Best Basis")
    #levlist2levlengths!(GP, BSc2f)

    # if we flattened dmatrix, then "unflatten" the expansion coefficients
    if fcols > 1
        dvecc2f = dmatrix2dvec(dmatrix0, GP, BSc2f)
    else
        dvecc2f = reshape(dvecc2f,(length(dvecc2f),1))
    end

    # Return things
    return dvecc2f, BSc2f
end # of function ghwt_c2f_bestbasis

"""
    (dvecf2c, BSf2c) = ghwt_f2c_bestbasis(dmatrix::Array{Float64,3}, GP::GraphPart;
                                          cfspec::Any = 1.0, flatten::Any = 1.0,
                                          j_start::Int = 1, j_end::Int = size(dmatrix,2),
                                          useParent::Bool = true)

Select the fine-to-coarse best basis from the matrix of GHWT expansion coefficients

### Input Arguments
* `dmatrix::Array{Float64,3}`: the matrix of expansion coefficients
* `GP::GraphPart`: an input GraphPart object
* `cfspec::Any`: the specification of cost functional to be used (default: 1.0, i.e., 1-norm)
* `flatten::Any`: the method for flattening vector-valued data to scalar-valued data (default: 1.0, i.e., 1-norm)
* `useParent::Bool`: the flag to indicate if we update the selected best basis
    subspace to the parent when parent and child have the same cost (default = true)

###  Output Arguments
* `dvecf2c::Matrix{Float64}`: the vector of expansion coefficients corresponding to the fine-to-coarse best basis
*  `BSf2c::BasisSpec`: a BasisSpec object which specifies the fine-to-coarse best basis
"""
function ghwt_f2c_bestbasis(dmatrix::Array{Float64,3}, GP::GraphPart;
                            cfspec::Any = 1.0, flatten::Any = 1.0,
                            j_start::Int = 1, j_end::Int = size(dmatrix,2),
                            useParent::Bool = true)

    # determine the cost functional to be used
    costfun = cost_functional(cfspec)

    # constants and dmatrix cleanup
    (N, jmax, fcols) = Base.size(dmatrix)
    if ~(1 <= j_start <= j_end <= jmax)
        error("j_start and j_end should satisfy the inequality 1 <= j_start <= j_end <= "*string(jmax))
    end

    dmatrix[ abs.(dmatrix) .< 10^2 * eps() ] .= 0

    # "flatten" dmatrix
    if fcols > 1
        dmatrix0 = deepcopy(dmatrix)      # keep the original dmatrix as dmatrix0
        if !isnothing(flatten)
            dmatrix = dmatrix_flatten(dmatrix, flatten)
        end
    end

    ## Find the best-basis from the GHWT fine-to-coarse dictionary

    # generate the fine-to-coarse GraphPart fields and coefficient matrix
    dmatrixf2c = fine2coarse!(GP, dmatrix = dmatrix, coefp = true)

    # allocate/initialize
    dvecf2c = dmatrixf2c[:, j_end, 1]
    levlistf2c = j_end * ones(Int, N)

    # set the tolerance with sign
    tol = 10^4 * eps() * (1 - useParent * 2)

    # perform the basis search
    for j = (j_end - 1):-1:j_start
        regioncount = count(!iszero, GP.rsf2c[:, j]) - 1
        for r = 1:regioncount
            indr = GP.rsf2c[r, j]:(GP.rsf2c[r + 1, j] - 1)
            ##### compute the cost of the current best basis
            costBB = costfun(dvecf2c[indr])
            costNEW = costfun(dmatrixf2c[indr, j, 1])
            if costBB >= costNEW + tol
                #(dvecf2c[indr], levlistf2c[indr]) = bbchange(dmatrixf2c[indr, j, 1], j)
                dvecf2c[indr] = dmatrixf2c[indr, j, 1]
                levlistf2c[indr] .= j
            end
        end
    end

    #levlistf2c = Array{Int}(levlistf2c[ levlistf2c .!= 0 ])
    levlistf2c = collect(enumerate(levlistf2c))
    BSf2c = BasisSpec(levlistf2c, c2f = false, description = "GHWT f2c Best Basis")
    #levlist2levlengths!(GP, BSf2c)
    # costf2c = costfun(dvecf2c)

    # if we flattened dmatrix, then "unflatten" the expansion coefficients
    if fcols > 1
        dvecf2c = dmatrix2dvec(dmatrix0, GP, BSf2c)
    else
        dvecf2c = reshape(dvecf2c,(length(dvecf2c),1))
    end

    # Return things
    return dvecf2c, BSf2c
end # of function ghwt_f2c_bestbasis

"""
    (dvec, BS) = ghwt_bestbasis(dmatrix::Array{Float64,3}, GP::GraphPart; cfspec::Any, flatten::Any = 1.0)

Select the overall best basis among the c2f and f2c best bases

### Input Arguments
* `dmatrix::Array{Float64,3}`: the matrix of expansion coefficients
* `GP::GraphPart`: an input GraphPart object
* `cfspec::Any`: the specification of cost functional to be used (default: 1.0, i.e., 1-norm)
* `flatten::Any`: the method for flattening vector-valued data to scalar-valued data (default: 1.0, i.e., 1-norm)

###  Output Arguments
* `dvec::Matrix{Float64}`: the vector of expansion coefficients corresponding to the best basis
* `BS::BasisSpec`: a BasisSpec object which specifies the best basis
"""
function ghwt_bestbasis(dmatrix::Array{Float64,3}, GP::GraphPart;
                        cfspec::Any = 1.0, flatten::Any = 1.0, j_start::Int = 1, j_end::Int = size(dmatrix,2))

    costfun = cost_functional(cfspec)
    jmax = size(dmatrix,2)
    if ~(1 <= j_start <= j_end <= jmax)
        error("j_start and j_end should satisfy the inequality 1 <= j_start <= j_end <= "*string(jmax))
    end

    # Select the c2f best basis
    (dvecc2f, BSc2f) = ghwt_c2f_bestbasis(dmatrix, GP, cfspec = cfspec, flatten = flatten, j_start = j_start, j_end = j_end)
    costc2f = costfun(dvecc2f)

    # Select the f2c best basis
    (dvecf2c, BSf2c) = ghwt_f2c_bestbasis(dmatrix, GP, cfspec = cfspec, flatten = flatten, j_start = jmax + 1 - j_end, j_end = jmax + 1 - j_start)
    costf2c = costfun(dvecf2c)

    # Compare the coarse-to-fine and fine-to-coarse best-bases
    if costc2f >= costf2c
        dvec = dvecf2c
        BS = BSf2c
    else
        dvec = dvecc2f
        BS = BSc2f
    end

    # Return things
    return dvec, BS
end # of function ghwt_bestbasis

"""
    function GHWT_jkl(GP::GraphPart, drow::Int, dcol::Int; c2f::Bool = true)

Generate the (j,k,l) indices for the GHWT basis vector corresponding to the coefficient dmatrix(drow,dcol)

### Input Arguments
* `GP`: a GraphPart object
* `drow`: the row of the expansion coefficient
* `dcol`: the column of the expansion coefficient

### Output Argument
* `j`: the level index of the expansion coefficient
* `k`: the subregion index of the expansion coefficient
* `l`: the tag of the expansion coefficient
"""
function GHWT_jkl(GP::GraphPart, drow::Int, dcol::Int; c2f::Bool = true)
    if c2f
        j = dcol - 1

        k = findfirst(GP.rs[:,dcol] .> drow)
        k -= 2

        if isempty(GP.tag)
            ghwt_core!(GP)
        end
        l = GP.tag[drow,dcol]
    else
        jmax = size(GP.rs, 2)
        j = jmax - dcol

        IX = fine2coarse!(GP; indp = true)
        k = findfirst(GP.rs[:,j+1] .> IX[drow,dcol])
        k -= 2

        l = GP.tagf2c[drow,dcol]
    end
    return j,k,l
end # of function GHWT_jkl

end # of module GHWT


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Structures related to the gravity models.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export ICGEM, GravityModel_Coefs

"""
    ICGEM

Structure to store the information contained in ICGEM files.

"""
@with_kw struct ICGEM
    # Mandatory keywords.
    product_type::String
    modelname::String
    gravity_constant::Float64
    radius::Float64
    max_degree::Int
    errors::Symbol

    # Optional keywords.
    tide_system::Symbol
    norm::Symbol

    # Coefficients.
    Clm::Matrix{Float64}           = Matrix{Float64}(undef,0,0)
    Slm::Matrix{Float64}           = Matrix{Float64}(undef,0,0)
    sigmaC::Matrix{Float64}        = Matrix{Float64}(undef,0,0)
    sigmaC_formal::Matrix{Float64} = Matrix{Float64}(undef,0,0)
    sigmaS::Matrix{Float64}        = Matrix{Float64}(undef,0,0)
    sigmaS_formal::Matrix{Float64} = Matrix{Float64}(undef,0,0)
end

"""
    GravityModel_Coefs{T}

Structure to store the information about a gravity model.

"""
@with_kw struct GravityModel_Coefs{T}
    name::String
    μ::T
    R0::T
    legendre_norm::Symbol
    n_max::Int
    C::Matrix{T}
    S::Matrix{T}
end