text:
"""
`FlipX <: ImageTransformation`

Description
============

Reverses the pixel order along the x-axis. In other words it
mirrors the Image along the y-axis (i.e. horizontally).

If created using the parameter `chance`, the operation will be
lifted into a `Either`. See the documentation of `Either` for more
information.

Usage
======

    FlipX()

    FlipX(chance)

Arguments
==========

- **`chance`** : The probability of the operation to be applied to
the image if the operation is used in `transform`. Must be within
the set [0, 1].

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Images
    using Augmentor

    # Create a test image with 4 pixel.
    img = grayim(UInt8[200 150; 50 1])

    # Apply the image operation to the image
    img_new = transform(FlipX(), img)

    # Check manually if the result is as expected
    @assert img_new == grayim(UInt8[50 1; 200 150])

see also
=========

`ImageTransformation`, `Either`, `transform`
"""
immutable FlipX <: ImageTransformation
end

FlipX(chance) = ProbableOperation(FlipX(); chance = chance)

Base.showcompact(io::IO, op::FlipX) = print(io, "Flip x-axis.")
multiplier(::FlipX) = 1

function transform{T<:AbstractImage}(op::FlipX, img::T)
    result = copyproperties(img, flipx(img))::T
    _log_operation!(result, op)
end


text:
using JuMP, Cbc, CSV, PrettyTables

#Importar datos de Excel
nodes_import = CSV.read("Node_Coordinates.csv", header = true)

#Definir que hay en cada columna del excel y cuantos nodos hay en total
num_nodes = size(nodes_import,1)
idcol = 1
Xcol = 2
Ycol = 3

#Generar la Matriz de distancias a partir del numero de nodos
distance_matrix = Array{Float64}(undef, (num_nodes, num_nodes))
#Llenar la matriz de distancias entre todos los puntos usando Pitagoras
for n in 1:num_nodes
    for s in 1:num_nodes
        distance_matrix[n,s] = sqrt((nodes_import[n,Ycol] - nodes_import[s,Ycol])^2 +
        (nodes_import[n,Xcol] - nodes_import[s,Xcol])^2)
    end
end

println(pretty_table(nodes_import))
println(distance_matrix)

mtzModel = Model(Cbc.Optimizer)
@variable(mtzModel, rutas[1:num_nodes,1:num_nodes], binary=true)
@variable(mtzModel, 2 <= u[2:num_nodes] <= num_nodes)
@objective(mtzModel, Min, sum(rutas[:,:].*distance_matrix[:,:]))
for i = 1:num_nodes
    @constraint(mtzModel, sum(rutas[i,:]) == 1)
    @constraint(mtzModel, sum(rutas[:,i]) == 1)
    @constraint(mtzModel, rutas[i,i] == 0)
end
for i = 2:num_nodes
    for j = 2:num_nodes
        @constraint(mtzModel, u[i]-u[j]+1 <= (num_nodes-1)*(1-rutas[i,j]))
    end
end
stats = JuMP.optimize!(mtzModel)
print(pretty_table(JuMP.value.(rutas)))

ext: jl

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


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

"""
    calculate_mortar_assembly(elements::Dict{Int, Vector{Int}},
                              element_types::Dict{Int, Symbol},
                              coords::Dict{Int, Vector{Float64}},
                              slave_element_ids::Vector{Int},
                              master_element_ids::Vector{Int})

Given data, calculate projection matrices `D` and `M`. This is the main
function of package. Relation between matrices is ``D u_s = M u_m``, where
``u_s`` is slave nodes and ``u_m`` master nodes.

# Example

Calculate mortar matrices for simple problem in README.md

```jldoctest ex1
Xs = Dict(1 => [0.0, 1.0], 2 => [5/4, 1.0], 3 => [2.0, 1.0])
Xm = Dict(4 => [0.0, 1.0], 5 => [1.0, 1.0], 6 => [2.0, 1.0])
coords = merge(Xm , Xs)
Es = Dict(1 => [1, 2], 2 => [2, 3])
Em = Dict(3 => [4, 5], 4 => [5, 6])
elements = merge(Es, Em)
element_types = Dict(1 => :Seg2, 2 => :Seg2, 3 => :Seg2, 4 => :Seg2)
slave_element_ids = [1, 2]
master_element_ids = [3, 4]
s, m, D, M = calculate_mortar_assembly(elements, element_types, coords,
                                       slave_element_ids, master_element_ids)

# output

([1, 2, 3], [4, 5, 6],
  [1, 1]  =  0.416667
  [2, 1]  =  0.208333
  [1, 2]  =  0.208333
  [2, 2]  =  0.666667
  [3, 2]  =  0.125
  [2, 3]  =  0.125
  [3, 3]  =  0.25,
  [1, 4]  =  0.366667
  [2, 4]  =  0.133333
  [1, 5]  =  0.25625
  [2, 5]  =  0.65
  [3, 5]  =  0.09375
  [1, 6]  =  0.00208333
  [2, 6]  =  0.216667
  [3, 6]  =  0.28125)

```

`s` and `m` contains slave and master dofs:
```jldoctest ex1
julia> s, m
([1, 2, 3], [4, 5, 6])
```

`D` is slave side mortar matrix:
```jldoctest ex1
julia> full(D[s,s])
3×3 Array{Float64,2}:
 0.416667  0.208333  0.0
 0.208333  0.666667  0.125
 0.0       0.125     0.25
```

`M` is master side mortar matrix:
```jldoctest ex1
julia> full(M[s,m])
3×3 Array{Float64,2}:
 0.366667  0.25625  0.00208333
 0.133333  0.65     0.216667
 0.0       0.09375  0.28125
```

"""
function calculate_mortar_assembly(elements, element_types, coords,
                                   slave_element_ids, master_element_ids)
    S = Int[]
    M = Int[]
    for sid in slave_element_ids
        push!(S, elements[sid]...)
    end
    for mid in master_element_ids
        push!(M, elements[mid]...)
    end
    S = sort(unique(S))
    M = sort(unique(M))
    ns = length(S)
    nm = length(M)
    normals = calculate_normals(elements, element_types, coords)
    segmentation = calculate_segments(slave_element_ids, master_element_ids,
                                      elements, element_types, coords, normals)
    D_I = Int[]
    D_J = Int[]
    D_V = Float64[]
    M_I = Int[]
    M_J = Int[]
    M_V = Float64[]
    for sid in slave_element_ids
        mids = [mid for (mid, xi) in segmentation[sid]]
        De, Me = calculate_mortar_matrices(sid, elements, element_types,
                                           coords, normals, segmentation)
        for (i,I) in enumerate(elements[sid])
            for (j,J) in enumerate(elements[sid])
                push!(D_I, I)
                push!(D_J, J)
                push!(D_V, De[i,j])
            end
        end
        for mid in mids
            for (i,I) in enumerate(elements[sid])
                for (j,J) in enumerate(elements[mid])
                   push!(M_I, I)
                   push!(M_J, J)
                   push!(M_V, Me[mid][i,j])
                end
            end
        end
    end

    # return global matrices + slave and master dofs
    Dg = sparse(D_I, D_J, D_V)
    Mg = sparse(M_I, M_J, M_V)
    return S, M, Dg, Mg
end


text:
using ConstructionBase: setproperties
import Unitful

export ControlNamelist,
    SystemNamelist,
    ElectronsNamelist,
    IonsNamelist,
    CellNamelist,
    DosNamelist,
    BandsNamelist,
    ElectronicTemperatureSetter,
    ElecTempSetter
export getxmldir, wfcfiles, getpseudodir

# From https://discourse.julialang.org/t/aliases-for-union-t-nothing-and-union-t-missing/15402/4
const Maybe{T} = Union{T,Nothing}

# The default values are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90.
"""
    ControlNamelist(calculation, title, verbosity, restart_mode, wf_collect, nstep, iprint, tstress, tprnfor, dt, outdir, wfcdir, prefix, lkpoint_dir, max_seconds, etot_conv_thr, forc_conv_thr, disk_io, pseudo_dir, tefield, dipfield, lelfield, nberrycyc, lorbm, lberry, gdir, nppstr, lfcpopt, gate)
    ControlNamelist(; kwargs...)
    ControlNamelist(::ControlNamelist; kwargs...)

Represent the `CONTROL` namelist of `pw.x`.
"""
struct ControlNamelist <: Namelist
    calculation::String
    title::String
    verbosity::String
    restart_mode::String
    wf_collect::Bool
    nstep::UInt
    iprint::UInt
    tstress::Bool
    tprnfor::Bool
    dt::Float64
    outdir::String
    wfcdir::String
    prefix::String
    lkpoint_dir::Bool
    max_seconds::Float64
    etot_conv_thr::Float64
    forc_conv_thr::Float64
    disk_io::String
    pseudo_dir::String
    tefield::Bool
    dipfield::Bool
    lelfield::Bool
    nberrycyc::UInt
    lorbm::Bool
    lberry::Bool
    gdir::UInt8
    nppstr::UInt
    lfcpopt::Bool
    gate::Bool
end # struct ControlNamelist
function ControlNamelist(;
    calculation = "scf",
    title = " ",
    verbosity = "low",
    restart_mode = "from_scratch",
    wf_collect = true,
    nstep = 50,
    iprint = 100000,
    tstress = false,
    tprnfor = false,
    dt = 20.0,
    outdir = "./",
    wfcdir = "./",
    prefix = "pwscf",
    lkpoint_dir = true,
    max_seconds = 10000000.0,
    etot_conv_thr = 1e-4,
    forc_conv_thr = 1e-3,
    disk_io = ifelse(calculation == "scf", "low", "medium"),
    pseudo_dir = raw"$HOME/espresso/pseudo/",
    tefield = false,
    dipfield = false,
    lelfield = false,
    nberrycyc = 1,
    lorbm = false,
    lberry = false,
    gdir = 1,  # The QE default value is `0`!
    nppstr = 0,
    lfcpopt = false,
    gate = false,
)
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1282-L1369.
    @assert calculation in ("scf", "nscf", "bands", "relax", "md", "vc-relax", "vc-md")
    @assert verbosity in ("high", "low", "debug", "medium", "default", "minimal")
    @assert restart_mode in ("from_scratch", "restart")
    @assert iprint >= 1
    @assert disk_io in ("high", "medium", "low", "none", "default")
    @assert dt >= 0
    # @assert !(lkpoint_dir && wf_collect) "`lkpoint_dir` currently doesn't work together with `wf_collect`!"
    @assert max_seconds >= 0
    @assert etot_conv_thr >= 0
    @assert forc_conv_thr >= 0
    @assert gdir in 1:3
    @assert !all((gate, tefield, !dipfield)) "`gate` cannot be used with `tefield` if dipole correction is not active!"
    @assert !all((gate, dipfield, !tefield)) "dipole correction is not active if `tefield = false`!"
    return ControlNamelist(
        calculation,
        title,
        verbosity,
        restart_mode,
        wf_collect,
        nstep,
        iprint,
        tstress,
        tprnfor,
        dt,
        outdir,
        wfcdir,
        prefix,
        lkpoint_dir,
        max_seconds,
        etot_conv_thr,
        forc_conv_thr,
        disk_io,
        pseudo_dir,
        tefield,
        dipfield,
        lelfield,
        nberrycyc,
        lorbm,
        lberry,
        gdir,
        nppstr,
        lfcpopt,
        gate,
    )
end
ControlNamelist(nml::ControlNamelist; kwargs...) = setproperties(nml; kwargs...)

getxmldir(nml::ControlNamelist) = expanduser(joinpath(nml.outdir, nml.prefix * ".save"))

wfcfiles(nml::ControlNamelist, n = 1) =
    [joinpath(getxmldir(nml), nml.prefix * ".wfc$i") for i in 1:n]

"""
    getpseudodir(nml::ControlNamelist)

Get the directory storing the pseudopotential files.
"""
getpseudodir(nml::ControlNamelist) = abspath(expanduser(nml.pseudo_dir))

"""
    SystemNamelist(ibrav, celldm, A, B, C, cosAB, cosAC, cosBC, nat, ntyp, nbnd, tot_charge, starting_charge, tot_magnetization, starting_magnetization, ecutwfc, ecutrho, ecutfock, nr1, nr2, nr3, nr1s, nr2s, nr3s, nosym, nosym_evc, noinv, no_t_rev, force_symmorphic, use_all_frac, occupations, one_atom_occupations, starting_spin_angle, degauss, smearing, nspin, noncolin, ecfixed, qcutz, q2sigma, input_dft, exx_fraction, screening_parameter, exxdiv_treatment, x_gamma_extrapolation, ecutvcut, nqx1, nqx2, nqx3, localization_thr, lda_plus_u, lda_plus_u_kind, Hubbard_U, Hubbard_J0, Hubbard_alpha, Hubbard_beta, starting_ns_eigenvalue, U_projection_type, edir, emaxpos, eopreg, eamp, angle1, angle2, constrained_magnetization, fixed_magnetization, lambda, report, lspinorb, assume_isolated, esm_bc, esm_w, esm_efield, esm_nfit, fcp_mu, vdw_corr, london, london_s6, london_c6, london_rvdw, london_rcut, ts_vdw_econv_thr, ts_vdw_isolated, xdm, xdm_a1, xdm_a2, space_group, uniqueb, origin_choice, rhombohedral, zgate, relaxz, block, block_1, block_2, block_height)
    SystemNamelist(; kwargs...)
    SystemNamelist(::SystemNamelist; kwargs...)

Represent the `SYSTEM` namelist of `pw.x`.
"""
struct SystemNamelist <: Namelist
    ibrav::Int8
    celldm::Vector{Maybe{Float64}}
    A::Float64
    B::Float64
    C::Float64
    cosAB::Float64
    cosAC::Float64
    cosBC::Float64
    nat::UInt
    ntyp::UInt8
    nbnd::UInt
    tot_charge::Float64
    starting_charge::Vector{Maybe{Float64}}
    tot_magnetization::Float64
    starting_magnetization::Vector{Maybe{Float64}}
    ecutwfc::Float64
    ecutrho::Float64
    ecutfock::Float64
    nr1::UInt
    nr2::UInt
    nr3::UInt
    nr1s::UInt
    nr2s::UInt
    nr3s::UInt
    nosym::Bool
    nosym_evc::Bool
    noinv::Bool
    no_t_rev::Bool
    force_symmorphic::Bool
    use_all_frac::Bool
    occupations::String
    one_atom_occupations::Bool
    starting_spin_angle::Bool
    degauss::Float64
    smearing::String
    nspin::UInt8
    noncolin::Bool
    ecfixed::Float64
    qcutz::Float64
    q2sigma::Float64
    input_dft::String
    exx_fraction::Float64
    screening_parameter::Float64
    exxdiv_treatment::String
    x_gamma_extrapolation::Bool
    ecutvcut::Float64
    nqx1::UInt
    nqx2::UInt
    nqx3::UInt
    localization_thr::Float64
    lda_plus_u::Bool
    lda_plus_u_kind::UInt8
    Hubbard_U::Vector{Maybe{Float64}}
    Hubbard_J0::Vector{Maybe{Float64}}
    Hubbard_alpha::Vector{Maybe{Float64}}
    Hubbard_beta::Vector{Maybe{Float64}}
    # Hubbard_J::Vector{Vector{Maybe{Float64}}}
    starting_ns_eigenvalue::Float64
    U_projection_type::String
    edir::UInt8
    emaxpos::Float64
    eopreg::Float64
    eamp::Float64
    angle1::Vector{Maybe{Float64}}
    angle2::Vector{Maybe{Float64}}
    constrained_magnetization::String
    fixed_magnetization::Vector{Maybe{Float64}}
    lambda::Float64
    report::UInt
    lspinorb::Bool
    assume_isolated::String
    esm_bc::String
    esm_w::Float64
    esm_efield::Float64
    esm_nfit::UInt
    fcp_mu::Float64
    vdw_corr::String
    london::Bool
    london_s6::Float64
    london_c6::Vector{Maybe{Float64}}
    london_rvdw::Vector{Maybe{Float64}}
    london_rcut::Float64
    ts_vdw_econv_thr::Float64
    ts_vdw_isolated::Bool
    xdm::Bool
    xdm_a1::Float64
    xdm_a2::Float64
    space_group::UInt8
    uniqueb::Bool
    origin_choice::UInt8
    rhombohedral::Bool
    zgate::Float64
    relaxz::Bool
    block::Bool
    block_1::Float64
    block_2::Float64
    block_height::Float64
end # struct SystemNamelist
function SystemNamelist(;
    ibrav = 127,
    celldm = zeros(6),  # Must specify
    A = 0.0,
    B = 0.0,
    C = 0.0,
    cosAB = 0.0,
    cosAC = 0.0,
    cosBC = 0.0,
    nat = 0,
    ntyp = 0,
    nbnd = 0,
    tot_charge = 0.0,
    starting_charge = [],
    tot_magnetization = -1.0,
    starting_magnetization = [],
    ecutwfc = 0.0,
    ecutrho = 4ecutwfc,
    ecutfock = ecutrho,
    nr1 = 0,
    nr2 = 0,
    nr3 = 0,
    nr1s = 0,
    nr2s = 0,
    nr3s = 0,
    nosym = false,
    nosym_evc = false,
    noinv = false,
    no_t_rev = false,
    force_symmorphic = false,
    use_all_frac = false,
    occupations = "fixed",
    one_atom_occupations = false,
    starting_spin_angle = false,
    degauss = 0.0,
    smearing = "gaussian",
    nspin = 1,
    noncolin = false,
    ecfixed = 0.0,
    qcutz = 0.0,
    q2sigma = 0.1,  # The default value in QE's source code is 0.01
    input_dft = "none",
    exx_fraction = 0.25,
    screening_parameter = 0.106,
    exxdiv_treatment = "gygi-baldereschi",
    x_gamma_extrapolation = true,
    ecutvcut = 0.0,
    nqx1 = 1,
    nqx2 = 1,
    nqx3 = 1,
    localization_thr = 0.0,  # This is only for QE 6.4
    lda_plus_u = false,
    lda_plus_u_kind = 0,
    Hubbard_U = [],
    Hubbard_J0 = [],
    Hubbard_alpha = [],
    Hubbard_beta = [],
    # Hubbard_J = [zeros(ntyp)]  ,  # The default value in QE's source code is just one 0.0
    starting_ns_eigenvalue = -1.0,  # It's actually a multidimensional array.
    U_projection_type = "atomic",
    edir = 1,
    emaxpos = 0.5,
    eopreg = 0.1,
    eamp = 0.001,  # The default value in QE's source code is 0.0
    angle1 = [],
    angle2 = [],
    constrained_magnetization = "none",
    fixed_magnetization = zeros(3),  # The default value in QE's source code is just one 0.0
    lambda = 1.0,
    report = 100,
    lspinorb = false,
    assume_isolated = "none",
    esm_bc = "pbc",
    esm_w = 0.0,
    esm_efield = 0.0,
    esm_nfit = 4,
    fcp_mu = 0.0,
    vdw_corr = "none",
    london = false,
    london_s6 = 0.75,
    london_c6 = [],
    london_rvdw = [],
    london_rcut = 200.0,
    ts_vdw_econv_thr = 1e-06,
    ts_vdw_isolated = false,
    xdm = false,
    xdm_a1 = 0.6836,  # The default value in QE's source code is 0.0
    xdm_a2 = 1.5045,  # The default value in QE's source code is 0.0
    space_group = 0,
    uniqueb = false,
    origin_choice = 1,
    rhombohedral = true,
    zgate = 0.5,
    relaxz = false,
    block = false,
    block_1 = 0.45,
    block_2 = 0.55,
    block_height = 0.1,  # The default value in QE's source code is 0.0
)
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1378-L1499.
    @assert ibrav in union(0:1:14, (-3, -5, -9, 91, -12, -13), 127)
    @assert ntyp <= 10 && ntyp <= nat
    @assert smearing in (
        "gaussian",
        "gauss",
        "methfessel-paxton",
        "m-p",
        "mp",
        "marzari-vanderbilt",
        "cold",
        "m-v",
        "mv",
        "fermi-dirac",
        "f-d",
        "fd",
    )
    @assert nspin in (1, 2, 4)
    @assert ecutwfc >= 0  # From `set_cutoff` in https://github.com/QEF/q-e/blob/7573301/PW/src/input.f90#L1597-L1639
    @assert ecutrho >= ecutwfc # From `set_cutoff`
    @assert ecfixed >= 0
    @assert qcutz >= 0
    @assert q2sigma >= 0
    @assert lda_plus_u_kind in 0:1
    @assert edir in 1:3
    @assert space_group in 0:230
    @assert origin_choice in 1:2
    @assert length(starting_charge) <= ntyp
    @assert length(starting_magnetization) <= ntyp
    @assert length(Hubbard_U) <= ntyp
    @assert length(Hubbard_J0) <= ntyp
    @assert length(Hubbard_alpha) <= ntyp
    @assert length(Hubbard_beta) <= ntyp
    # @assert all(length(x) <= ntyp for x in Hubbard_J)
    @assert length(angle1) <= ntyp
    @assert length(angle2) <= ntyp
    @assert length(fixed_magnetization) <= 3
    @assert length(london_c6) <= ntyp
    @assert length(london_rvdw) <= ntyp
    @assert exxdiv_treatment in
            ("gygi-baldereschi", "gygi-bald", "g-b", "vcut_ws", "vcut_spherical", "none")
    @assert !(x_gamma_extrapolation && exxdiv_treatment in ("vcut_ws", "vcut_spherical")) "`x_gamma_extrapolation` cannot be used with `vcut`!"
    return SystemNamelist(
        ibrav,
        celldm,
        A,
        B,
        C,
        cosAB,
        cosAC,
        cosBC,
        nat,
        ntyp,
        nbnd,
        tot_charge,
        starting_charge,
        tot_magnetization,
        starting_magnetization,
        ecutwfc,
        ecutrho,
        ecutfock,
        nr1,
        nr2,
        nr3,
        nr1s,
        nr2s,
        nr3s,
        nosym,
        nosym_evc,
        noinv,
        no_t_rev,
        force_symmorphic,
        use_all_frac,
        occupations,
        one_atom_occupations,
        starting_spin_angle,
        degauss,
        smearing,
        nspin,
        noncolin,
        ecfixed,
        qcutz,
        q2sigma,
        input_dft,
        exx_fraction,
        screening_parameter,
        exxdiv_treatment,
        x_gamma_extrapolation,
        ecutvcut,
        nqx1,
        nqx2,
        nqx3,
        localization_thr,
        lda_plus_u,
        lda_plus_u_kind,
        Hubbard_U,
        Hubbard_J0,
        Hubbard_alpha,
        Hubbard_beta,
        starting_ns_eigenvalue,
        U_projection_type,
        edir,
        emaxpos,
        eopreg,
        eamp,
        angle1,
        angle2,
        constrained_magnetization,
        fixed_magnetization,
        lambda,
        report,
        lspinorb,
        assume_isolated,
        esm_bc,
        esm_w,
        esm_efield,
        esm_nfit,
        fcp_mu,
        vdw_corr,
        london,
        london_s6,
        london_c6,
        london_rvdw,
        london_rcut,
        ts_vdw_econv_thr,
        ts_vdw_isolated,
        xdm,
        xdm_a1,
        xdm_a2,
        space_group,
        uniqueb,
        origin_choice,
        rhombohedral,
        zgate,
        relaxz,
        block,
        block_1,
        block_2,
        block_height,
    )
end
SystemNamelist(nml::SystemNamelist; kwargs...) = setproperties(nml; kwargs...)

"""
    ElectronsNamelist(electron_maxstep, scf_must_converge, conv_thr, adaptive_thr, conv_thr_init, conv_thr_multi, mixing_mode, mixing_beta, mixing_ndim, mixing_fixed_ns, diagonalization, ortho_para, diago_thr_init, diago_cg_maxiter, diago_david_ndim, diago_full_acc, efield, efield_cart, efield_phase, startingpot, startingwfc, tqr)
    ElectronsNamelist(; kwargs...)
    ElectronsNamelist(::ElectronsNamelist; kwargs...)

Represent the `ELECTRONS` namelist of `pw.x`.
"""
struct ElectronsNamelist <: Namelist
    electron_maxstep::UInt
    scf_must_converge::Bool
    conv_thr::Float64
    adaptive_thr::Bool
    conv_thr_init::Float64
    conv_thr_multi::Float64
    mixing_mode::String
    mixing_beta::Float64
    mixing_ndim::UInt
    mixing_fixed_ns::UInt
    diagonalization::String
    ortho_para::UInt
    diago_thr_init::Float64
    diago_cg_maxiter::UInt
    diago_david_ndim::UInt
    diago_full_acc::Bool
    efield::Float64
    efield_cart::Vector{Maybe{Float64}}
    efield_phase::String
    startingpot::String  # This depends on `calculation`
    startingwfc::String  # This depends on `calculation`
    tqr::Bool
end # struct ElectronsNamelist
function ElectronsNamelist(;
    electron_maxstep = 100,
    scf_must_converge = true,
    conv_thr = 1e-6,
    adaptive_thr = false,
    conv_thr_init = 1e-3,
    conv_thr_multi = 0.1,
    mixing_mode = "plain",
    mixing_beta = 0.7,
    mixing_ndim = 8,
    mixing_fixed_ns = 0,
    diagonalization = "david",
    ortho_para = 0,
    diago_thr_init = 0.0,
    diago_cg_maxiter = 20,
    diago_david_ndim = 4,
    diago_full_acc = false,
    efield = 0.0,
    efield_cart = zeros(3),
    efield_phase = "none",
    startingpot = "atomic",  # This depends on `calculation`
    startingwfc = "atomic+random",  # This depends on `calculation`
    tqr = false,
)
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1508-L1543.
    @assert mixing_mode in ("plain", "TF", "local-TF")
    @assert diagonalization in ("david", "cg", "cg-serial", "david-serial", "ppcg")  # Different from docs
    @assert efield_phase in ("read", "write", "none")
    @assert startingpot in ("atomic", "file")
    @assert startingwfc in ("atomic", "atomic+random", "random", "file")
    return ElectronsNamelist(
        electron_maxstep,
        scf_must_converge,
        conv_thr,
        adaptive_thr,
        conv_thr_init,
        conv_thr_multi,
        mixing_mode,
        mixing_beta,
        mixing_ndim,
        mixing_fixed_ns,
        diagonalization,
        ortho_para,
        diago_thr_init,
        diago_cg_maxiter,
        diago_david_ndim,
        diago_full_acc,
        efield,
        efield_cart,
        efield_phase,
        startingpot,
        startingwfc,
        tqr,
    )
end
ElectronsNamelist(nml::ElectronsNamelist; kwargs...) = setproperties(nml; kwargs...)

"""
    IonsNamelist(ion_dynamics, ion_positions, pot_extrapolation, wfc_extrapolation, remove_rigid_rot, ion_temperature, tempw, tolp, delta_t, nraise, refold_pos, upscale, bfgs_ndim, trust_radius_max, trust_radius_min, trust_radius_ini, w_1, w_2)
    IonsNamelist(; kwargs...)
    IonsNamelist(::IonsNamelist; kwargs...)

Represent the `IONS` namelist of `pw.x`.

Input this namelist only if `calculation` is `"relax"`, `"md"`, `"vc-relax"`, or `"vc-md"`.
"""
struct IonsNamelist <: Namelist
    ion_dynamics::String
    ion_positions::String
    pot_extrapolation::String
    wfc_extrapolation::String
    remove_rigid_rot::Bool
    ion_temperature::String
    tempw::Float64
    tolp::Float64
    delta_t::Float64
    nraise::UInt
    refold_pos::Bool
    upscale::Float64
    bfgs_ndim::UInt
    trust_radius_max::Float64
    trust_radius_min::Float64  # The default value in QE's source code is 0.0001
    trust_radius_ini::Float64
    w_1::Float64
    w_2::Float64
end # struct IonsNamelist
function IonsNamelist(;
    ion_dynamics = "none",
    ion_positions = "default",
    pot_extrapolation = "atomic",
    wfc_extrapolation = "none",
    remove_rigid_rot = false,
    ion_temperature = "not_controlled",
    tempw = 300.0,
    tolp = 100.0,
    delta_t = 1.0,
    nraise = 1,
    refold_pos = false,
    upscale = 100.0,
    bfgs_ndim = 1,
    trust_radius_max = 0.8,
    trust_radius_min = 1e-3,  # The default value in QE's source code is 0.0001
    trust_radius_ini = 0.5,
    w_1 = 0.01,
    w_2 = 0.5,
)
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1552-L1585.
    @assert ion_dynamics in
            ("none", "bfgs", "damp", "verlet", "langevin", "langevin-smc", "beeman")
    @assert ion_positions in ("default", "from_input")
    @assert pot_extrapolation in ("none", "atomic", "first_order", "second_order")
    @assert wfc_extrapolation in ("none", "first_order", "second_order")
    @assert ion_temperature in (
        "rescaling",
        "rescale-v",
        "rescale-T",
        "reduce-T",
        "berendsen",
        "andersen",
        "initial",
        "not_controlled",
    )
    @assert tempw > 0
    return IonsNamelist(
        ion_dynamics,
        ion_positions,
        pot_extrapolation,
        wfc_extrapolation,
        remove_rigid_rot,
        ion_temperature,
        tempw,
        tolp,
        delta_t,
        nraise,
        refold_pos,
        upscale,
        bfgs_ndim,
        trust_radius_max,
        trust_radius_min,
        trust_radius_ini,
        w_1,
        w_2,
    )
end
IonsNamelist(nml::IonsNamelist; kwargs...) = setproperties(nml; kwargs...)

"""
    CellNamelist(cell_dynamics, press, wmass, cell_factor, press_conv_thr, cell_dofree)
    CellNamelist(; kwargs...)
    CellNamelist(::CellNamelist; kwargs...)

Represent the `CELL` namelist of `pw.x`.

Input this namelist only if `calculation` is `"vc-relax"` or `"vc-md"`.
"""
struct CellNamelist <: Namelist
    cell_dynamics::String
    press::Float64
    wmass::Float64
    cell_factor::Float64
    press_conv_thr::Float64
    cell_dofree::String
end # struct CellNamelist
function CellNamelist(;
    cell_dynamics = "none",
    press = 0.0,
    wmass = 0.0,
    cell_factor = 0.0,
    press_conv_thr = 0.5,
    cell_dofree = "all",
)
    # These checks are from https://github.com/QEF/q-e/blob/4132a64/Modules/read_namelists.f90#L1596-L1625.
    @assert cell_dynamics in ("none", "sd", "damp-pr", "damp-w", "bfgs", "pr", "w")
    @assert wmass >= 0
    @assert cell_dofree in (
        "all",
        "ibrav",
        "x",
        "y",
        "z",
        "xy",
        "xz",
        "yz",
        "xyz",
        "shape",
        "volume",
        "2Dxy",
        "2Dshape",
        "epitaxial_ab",  # New in 6.4
        "epitaxial_ac",  # New in 6.4
        "epitaxial_bc",  # New in 6.4
    )
    return CellNamelist(
        cell_dynamics,
        press,
        wmass,
        cell_factor,
        press_conv_thr,
        cell_dofree,
    )
end
CellNamelist(nml::CellNamelist; kwargs...) = setproperties(nml; kwargs...)

# The following default values are picked from `<QE source>/PP/src/dos.f90`
"""
    DosNamelist(prefix, outdir, ngauss, degauss, Emin, Emax, DeltaE, fildos)
    DosNamelist(; kwargs...)
    DosNamelist(::DosNamelist; kwargs...)

Represent the `DOS` namelist of `dos.x`.
"""
struct DosNamelist <: Namelist
    prefix::String
    outdir::String
    ngauss::Int
    degauss::Float64
    Emin::Float64
    Emax::Float64
    DeltaE::Float64
    fildos::String
end # struct DosNamelist
function DosNamelist(;
    prefix = "pwscf",
    outdir = "./",
    ngauss = 0,
    degauss = 0.0,
    Emin = -1000000.0,
    Emax = 1000000.0,
    DeltaE = 0.01,
    fildos = "$(prefix).dos",
)
    @assert ngauss in (0, 1, -1, -99)
    return DosNamelist(prefix, outdir, ngauss, degauss, Emin, Emax, DeltaE, fildos)
end
DosNamelist(nml::DosNamelist; kwargs...) = setproperties(nml; kwargs...)

# The following default values are picked from `<QE source>/PP/src/bands.f90`
"""
    BandsNamelist(prefix, outdir, filband, spin_component, lsigma, lp, filp, lsym, no_overlap, plot_2d, firstk, lastk)
    BandsNamelist(; kwargs...)
    BandsNamelist(::BandsNamelist; kwargs...)

Represent the `BANDS` namelist of `bands.x`.
"""
struct BandsNamelist <: Namelist
    prefix::String
    outdir::String
    filband::String
    spin_component::UInt8
    lsigma::Vector{Maybe{Bool}}  # The default value in QE's source code is just one `false`
    lp::Bool
    filp::String
    lsym::Bool
    no_overlap::Bool
    plot_2d::Bool
    firstk::UInt
    lastk::UInt
end # struct BandsNamelist
function BandsNamelist(;
    prefix = "pwscf",
    outdir = "./",
    filband = "bands.out",
    spin_component = 1,
    lsigma = falses(3),  # The default value in QE's source code is just one `false`
    lp = false,
    filp = "p_avg.dat",
    lsym = true,
    no_overlap = true,
    plot_2d = false,
    firstk = 0,
    lastk = 10000000,
)
    @assert spin_component in 1:2
    return BandsNamelist(
        prefix,
        outdir,
        filband,
        spin_component,
        lsigma,
        lp,
        filp,
        lsym,
        no_overlap,
        plot_2d,
        firstk,
        lastk,
    )
end
BandsNamelist(nml::BandsNamelist; kwargs...) = setproperties(nml; kwargs...)

@batteries ControlNamelist eq = true hash = true
@batteries SystemNamelist eq = true hash = true
@batteries ElectronsNamelist eq = true hash = true
@batteries IonsNamelist eq = true hash = true
@batteries CellNamelist eq = true hash = true
@batteries DosNamelist eq = true hash = true
@batteries BandsNamelist eq = true hash = true

function (x::VerbositySetter)(control::ControlNamelist)
    if x.v == "high"
        @set! control.verbosity = "high"
        @set! control.wf_collect = true
        @set! control.tstress = true
        @set! control.tprnfor = true
        @set! control.disk_io = "high"
    else
        @set! control.verbosity = "low"
        @set! control.wf_collect = false
        @set! control.tstress = false
        @set! control.tprnfor = false
        @set! control.disk_io = "low"
    end
    return control
end

struct ElectronicTemperatureSetter{T<:Number} <: Setter
    temp::T
end
function (x::ElectronicTemperatureSetter)(system::SystemNamelist)
    @set! system.occupations = "smearing"
    @set! system.smearing = "fermi-dirac"
    @set! system.degauss = degauss(x.temp)
    return system
end
degauss(temp::Real) = temp
degauss(temp::Unitful.Temperature) = ustrip(u"K", temp) / 315775.02480407 * 2
degauss(temp::Unitful.Energy) = ustrip(u"Ry", temp)
degauss(temp::Unitful.Frequency) = ustrip(u"Hz", temp) / 6579683920502000.0 * 2
degauss(temp::Unitful.Wavenumber) = ustrip(u"m^-1", temp) / 21947463.13632 * 2

const ElecTempSetter = ElectronicTemperatureSetter

# _coupledargs(::Type{ControlNamelist}) = (:calculation => :disk_io,)
# _coupledargs(::Type{SystemNamelist}) = (:ecutwfc => :ecutrho, :ecutrho => :ecutfock)
# _coupledargs(::Type{DosNamelist}) = (:prefix => :fildos,)

groupname(::Type{ControlNamelist}) = "CONTROL"
groupname(::Type{SystemNamelist}) = "SYSTEM"
groupname(::Type{ElectronsNamelist}) = "ELECTRONS"
groupname(::Type{IonsNamelist}) = "IONS"
groupname(::Type{CellNamelist}) = "CELL"


text:
"""
---
Constant Modulus Zero Autocorrelation (CAZAC) sequence generation. Sequence can be used as preamble sequences for OFDM systems. The function generate the ZC sequence in frequency domain.
Generates a Zadoff-Chu  sequence on allocated subcarrier mapped on nFFT vector with generated kernel muPSS and power boost zcBoost
# --- Syntax
zcSeq	  = genZCSequence(nFFT,allocatedSubcarrier,muPSS=0,zcBoost=0)
# --- Input parameters
- nFFT				  : FFT size for output [Int]
- allocatedSubcarrier : vector of allocated subcarrier [Array{Int,L}].
- muPSS				  : Kernel for ZC sequence [Int] -- default 0
- zcBoost			  : Power boost (in dB) applied to sequence [Float32] -- default 0.
# --- Output parameters
- zcSeq				  : ZC sequence [Array{Complex{Float64}},nFFT]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function genZCSequence(nFFT::Int,allocatedSubcarrier::Array{Int},muPss=0::Int,zcBoost=0)
	# --- Getting ZC length
	nbSubcarriers	= length(allocatedSubcarrier);
	# --- Create a subcarrier index vector
	nI				= (1:nbSubcarriers);
	# --- Create frequency container
	seqZC	= zeros(Complex{Float64},nFFT);
	# --- Fill frequency location with ZC sequence
	seqZC[allocatedSubcarrier] = (10).^(zcBoost./20) .* exp.(im*pi*muPss*nI.*(nI.+mod(nbSubcarriers,2))/nbSubcarriers);
	# --- Return generated sequence
	return seqZC;
end


text:
## A domain representing a union of subdomains

canonicaldomain(d::UnionDomain) = d  # we could map all to canonical, but then there would be overlap

isambiguous(d::UnionDomain) = isempty(d.domains)
convert(::Type{UnionDomain{DD,T}},::AnyDomain) where {DD,T} =
    UnionDomain{DD,T}(map(D->D(AnyDomain()),DD.parameters))
convert(::Type{IT},::AnyDomain) where {IT<:UnionDomain} = UnionDomain(tuple())

#support tuple set
components(d::AbstractVector) = d

pieces(d::UnionDomain) = components(d)
piece(d::UnionDomain,k) = pieces(d)[k]
npieces(d::UnionDomain) = length(components(d))

union(::AnyDomain, d::UnionDomain) = d
union(d::UnionDomain, ::AnyDomain) = d

arclength(d::UnionDomain) = mapreduce(arclength,+,d.domains)

reverseorientation(d::UnionDomain) = UnionDomain(reverse(map(reverseorientation,d.domains)))

leftendpoint(d::UnionDomain) = leftendpoint(first(components(d)))
rightendpoint(d::UnionDomain) = rightendpoint(last(components(d)))

# determine the number of points per piece
function components_npoints(d, n::Int)
    N = ncomponents(d)
    k = n ÷ N
    r = n - N*k
    [fill(k+1, r); fill(k, N-r)]
end

pieces_npoints(d, n::Int) = components_npoints(d, n)

points(d::UnionDomain,n) = vcat(points.(pieces(d), pieces_npoints(d,n))...)

rand(d::UnionDomain) = rand(component(d,rand(1:length(d))))
checkpoints(d::UnionDomain) = mapreduce(checkpoints,union,d.domains)

## to move over
function Base.merge(d1::UnionDomain,m::IntervalOrSegment)
    T = float(promote_type(eltype(d1), eltype(m)))
    ret=d1.domains

    for k=length(ret):-1:1
        it=intersect(ret[k],m)
        if arclength(it) ≠ 0
            sa=setdiff(ret[k],it)
            m=setdiff(m,it)
            if arclength(sa) ≤ eps(T)
                ret = [ret[1:k-1]...; it; ret[k+1:end]...]
            else
                ret = [ret[1:k-1]...; sa; it; ret[k+1:end]...]
            end
            if arclength(m) == 0
                break
            end
        end
    end
    if !isempty(m)
        ret = [ret...; m]
    end

    UnionDomain(sort!(ret,by=leftendpoint))
end

function merge(d1::UnionDomain,d2::UnionDomain)
    ret=d1
    for m in d2.domains
        ret=merge(ret,m)
    end
    ret
end


text:
import JSON
using Sockets
include("header.jl")

#global table_list
table_list = Dict()

function join_table(p::Player)
    global table_list
    let key_output
        key_output = ""
        for key in keys(table_list)
            table = table_list[key]
            if round(Int64, time()) - table["created"] > 1200
                delete!(table_list, key)
            elseif table["game"] != p.current_game || table["minimum"] > p.balance
                continue
            else
                key_output *= "$key\n"
            end
        end
        if key_output == ""
            print_dict(p, "table_3")
        else
            print_dict(p, "table_1")
            write(p.socket, key_output)
        end
    end

    while true
        print_dict(p, "table_2")
        s = readline(p.socket)
        if s == ""
            print_dict(p, "repeat")
            continue
        elseif s == "l"
            return false
        elseif s == "c"
            table = create_table(p)
            write(p.socket, "You approach the $(p.current_game) table $(table["name"]). The dealer smiles at you and slightly nods his head as a greeting.\n")
            break
        elseif haskey(table_list, s) && table_list[s]["game"] == p.current_game
            table = table_list[s]
            if table["minimum"] <= p.balance
                write(p.socket, "You approach the $(p.current_game) table $(table["name"]). The dealer smiles at you and slightly nods his head as a greeting.\n")
                break
            else
                print_dict(p, "table_4")
                table = table_list[s]
                s = readline(p.socket)
                if table["passphrase"] == s
                    write(p.socket, "You approach the $(p.current_game) table $(table["name"]). The dealer smiles at you and slightly nods his head as a greeting.\n")
                    break;
                else
                    print_dict(p, "table_5")
                    continue
                end
            end
        else
            print_dict(p, "repeat")
            continue
        end
    end
    return true
end

function create_table(p::Player)
    global table_list

    print_dict(p, "table_12")
    let key, name, minimum, passphrase
        while true
            key = readline(p.socket)
            if length(key) > 30
                print_dict(p, "table_13")
                continue
            elseif haskey(table_list,key)
                if round(Int64, time()) - table_list[key]["created"] <= 1200
                    print_dict(p, "table_14")
                    continue
                else
                    break
                end
            else
                break
            end
        end

        print_dict(p, "table_6")
        while true
            name = readline(p.socket)

            if length(name) > 64
                print_dict(p, "table_7")
                continue
            else
                break
            end
        end
        print_dict(p, "table_8")
        while true
            s = readline(p.socket)
            minimum = tryparse(Int, s)
            if minimum == nothing || minimum < 1
                print_dict(p, "table_9")
                continue
            else
                break
            end
        end
        print_dict(p, "table_10")
        while true
            passphrase = readline(p.socket)
            if length(passphrase) > 24
                print_dict(p, "table_11")
            else
                break
            end
        end

        global table_list[key] = Dict("name" => name, "minimum" => minimum, "passphrase" => passphrase, "game" => p.current_game, "created" => round(Int64, time()))

        print_dict(p, "table_15")
        print_dict(p, "spacer")

        return table_list[key]
    end
end


text:
function load_data()
    path = pwd()
    XFA = XLSX.readxlsx(path*"/data/algorithm_data_ns.xlsx")
    XFI = XLSX.readxlsx(path*"/data/investment_data_ns.xlsx")
    XFO_Electricity = XLSX.readxlsx(path*"/data/operation_electricity_data_ns.xlsx")
    XFO_heat = XLSX.readxlsx(path*"/data/operation_heat_data_ns.xlsx")
    XFO_NG = XLSX.readxlsx(path*"/data/operation_NG_data_ns.xlsx")
    XFO_W = XLSX.readxlsx(path*"/data/operation_water_data_ns.xlsx")
    XFO_O = XLSX.readxlsx(path*"/data/operation_oil_data_ns.xlsx")
    XFO_hydrogen = XLSX.readxlsx(path*"/data/operation_hydrogen_data_ns.xlsx")
    XFU = XLSX.readxlsx(path*"/data/uncertainty_data_ns.xlsx")

    unc= u_type(0,Vector{Tuple{String,String}}(),String[],String[],zeros(0,0),zeros(0,0),Dict{String,Int64}(),Dict{String,Int64}())
    ms = ms_type(String[],String[],String[],String[],String[],String[],String[],String[],String[],0:0,0:0)
    mp = mp_type(0,
                zeros(0,0),zeros(0,0,0),zeros(0,0,0),zeros(0,0,0),
                zeros(0),
                zeros(0,0),zeros(0),
                zeros(0),zeros(0),
                Vector{Array{Int64,1}}(),Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),
                Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),
                Dict{String,Int64}(),Dict{String,Int64}())
    ps = ps_type(0:0,0:0,Array{Int64,1}(),Array{UnitRange{Int64},1}(),String[],String[],String[],String[],String[],
                String[],String[],String[],String[],String[],String[],String[],String[],
                String[],String[],String[],String[],String[],
                String[],
                String[],String[],
                String[],String[],String[])
    pp = pp_type(0.,0.,
                0.,zeros(0,0),zeros(0,0),0.,zeros(0),0.,0.,0.,zeros(0,0),0.,0.,0.,0.,0.,0.,0.,zeros(0,0),zeros(0),zeros(0),
                0.,0.,0.,0.,0.,
                zeros(0,0),zeros(0,0),0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
                zeros(0,0), 0.,
                zeros(0,0),zeros(0,0),zeros(0,0),0.,0.,
                0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
                0.,0.,0.,0.,
                Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}(),Dict{String,Int64}())

    # */ ------------------------- structure ----------------------------------- /* #
    sheet = XFU["structure"]

    ndict = Dict(:t0 =>"present",
                 :t5 =>"5 years")

    tarray = [:t0,:t5]

    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    t0__col = findfirst(col_names .== ndict[:t0 ])[2]
    t5__col = findfirst(col_names .== ndict[:t5 ])[2]
    t0__value = convert(Array{Int64,1},sheet[XLSX.encode_column_number(t0__col)*"2"*":"*XLSX.encode_column_number(t0__col)*"$stpr"][:,1])
    t5__value = convert(Array{Int64,1},sheet[XLSX.encode_column_number(t5__col)*"2"*":"*XLSX.encode_column_number(t5__col)*"$stpr"][:,1])
    nω1 = length(unique(t0__value))
    nω2 = length(unique(t5__value))

    d_dict = Dict{Symbol,Array{Tuple{Int64,UnitRange{Int64},Int64},1}}()
    for sym in keys(ndict)
        tc = findfirst(col_names .== ndict[sym])[2]
        tv = convert(Array{Int64,1},sheet[XLSX.encode_column_number(tc)*"2"*":"*XLSX.encode_column_number(tc)*"$stpr"][:,1])
        ay = [(tc,1:length(tv),tc)]
        d_dict[sym] = ay
    end
    d_array = [(1,2)]
    s_dict = Dict{Symbol,Array{Tuple{Int64,Int64,Int64},1}}()
    for sym in keys(ndict)
        tc = findfirst(col_names .== ndict[sym])[2]
        tv = convert(Array{Int64,1},sheet[XLSX.encode_column_number(tc)*"2"*":"*XLSX.encode_column_number(tc)*"$stpr"][:,1])
        tn = unique(tv)
        ay = Array{Tuple{Int64,Int64,Int64},1}(undef,length(tn))
        for i in 1:length(tn)
            tr = findfirst(tv .== tn[i])
            ay[i] = (tn[i],tr,tc)
        end
        s_dict[sym] = ay
    end
    s_array = Array{Tuple{Int64,Int64},1}(undef,length(t0__value))
    for i in 1:length(t0__value)
        s_array[i] = (t0__value[i],t5__value[i])
    end

    ################ parameter ##################
    sheet = XFU["parameters"]
    ndict = Dict(:shn =>"sheet_name",
                :rc  =>"location",
                :type=>"type")
    udict = Dict(:μD  =>"demand_scaling",
                :μE =>"co2_budget_scaling",
                :CCO2 =>"co2_price")
    uarray = [:μD,:μE,:CCO2]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    row_names = sheet[XLSX.encode_column_number(sttc)*"3"*":"*XLSX.encode_column_number(sttc)*"$stpr"]
    l_col = findfirst(col_names .== ndict[:rc  ])[2]
    t_col = findfirst(col_names .== ndict[:type])[2]
    l_ay = convert(Array{Int64,1},sheet[XLSX.encode_column_number(l_col)*"3"*":"*XLSX.encode_column_number(l_col)*"$stpr"][:,1])
    t_ay = convert(Array{Int64,1},sheet[XLSX.encode_column_number(t_col)*"3"*":"*XLSX.encode_column_number(t_col)*"$stpr"][:,1])
    NI1 = nω1^sum(t_ay)
    NI2 = nω2^sum(t_ay)
    NI0 = NI1
    NI = NI2
    th = String[];tc =String[];
    for i in 1:length(l_ay)
        if l_ay[i]==0
            push!(th,String(row_names[i]))
        elseif l_ay[i]==1
            push!(tc,String(row_names[i]))
    end
    end
    nh = length(l_ay)-sum(l_ay)
    nc = sum(l_ay)

    scol = findfirst(col_names .== ndict[:shn])[2]
    svalue = convert(Array{String,1},sheet[XLSX.encode_column_number(scol)*"3"*":"*XLSX.encode_column_number(scol)*"$stpr"][:,1])
    ddict = Dict{Symbol,Dict{Int64,Float64}}()
    adict = Dict{Symbol,Array{Tuple{Int64,Int64},1}}()

    for sym in uarray
        row = findfirst(row_names .== udict[sym])[1]
        sname = svalue[row]
        sheet = XFU[sname]
        stt,stp = sheet.dimension.start,sheet.dimension.stop
        sttr,sttc = stt.row_number,stt.column_number
        stpr,stpc = stp.row_number,stp.column_number
        values = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc)*"3"*":"*XLSX.encode_column_number(stpc)*"$stpr"])
        vdict = Dict{Int64,Float64}()
        if t_ay[row] == 0
        for tm in tarray
            dic_tm = d_dict[tm]
            for j in 1:length(dic_tm)
                v = values[dic_tm[j][2]]
                vdict[dic_tm[j][1]] = sum(v)/length(v)
            end
        end
        varray = d_array
        else
        for tm in tarray
            dic_tm = s_dict[tm]
            for j in 1:length(dic_tm)
                vdict[dic_tm[j][1]] = values[dic_tm[j][2]]
            end
        end
        varray = s_array
        end
        ddict[sym] = vdict
        adict[sym] = varray
    end

    glb_array = Array{typeof(ntuple(x->0,length(uarray))),2}(undef,0,2)
    for i1 in adict[uarray[1]] for i2 in adict[uarray[2]] for i3 in adict[uarray[3]]
        ay = Array{typeof(ntuple(x->0,length(uarray))),2}(undef,1,2)
        for j in 1:2 ay[j] = (i1[j],i2[j],i3[j]) end
        glb_array = vcat(glb_array,ay)
    end end end
    num_dict = Dict{Int64,typeof(glb_array[1,1])}()
    dict_num = Dict{typeof(glb_array[1,1]),Int64}()
    i_n = 1
    for j in 1:2
        for el in unique(glb_array[:,j])
            num_dict[i_n] = el
            dict_num[el] = i_n
            i_n += 1
        end
    end
    nr = i_n-2
    H,C = zeros(nr,nh),zeros(nr,nc)
    for i in 1:nr
        ih,ic=1,1
        dc = num_dict[i+1]
        for j in 1:length(dc)
            if l_ay[j] == 0
                H[i,ih] = ddict[uarray[j]][dc[j]]
                ih += 1
            end
            if l_ay[j] == 1
                C[i,ic] = ddict[uarray[j]][dc[j]]
                ic += 1
            end
        end
    end

    itoi0 = Vector{Array{Int64,1}}()
    for i in 1:NI1
        rw = findfirst(x->x==num_dict[i+NI1],glb_array[:,2])
        push!(itoi0,[dict_num[glb_array[rw,1]]])
    end

    setfield!(unc,:ni,NI)
    setfield!(unc,:th,th)
    setfield!(unc,:tc,tc)
    setfield!(unc,:h,H)
    setfield!(unc,:c,C)

    map_h=Dict{String,Int64}();map_c=Dict{String,Int64}();
    for i in 1:length(th)
        map_h[th[i]]=i
    end
    for i in 1:length(tc)
        map_c[tc[i]]=i
    end
    setfield!(unc,:map_h,map_h)
    setfield!(unc,:map_c,map_c)
    setfield!(ms,:I0,1:NI0)
    setfield!(ms,:I ,1:NI)
    setfield!(mp,:π0,vcat(ones(NI1)./NI1))
    setfield!(mp,:π ,vcat(ones(NI2)./NI2))
    setfield!(mp,:map,itoi0)

    ########## read nodes sets ##########
    sheet = XFI["nodes"]
    ndict = Dict(:node =>"node",
                :type   =>"type")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    fcol = findfirst(col_names .== ndict[:node])[2]
    tcol = findfirst(col_names .== ndict[:type])[2]
    fvalue = convert(Array{String,1},sheet[XLSX.encode_column_number(fcol)*"3"*":"*XLSX.encode_column_number(fcol)*"$stpr"][:,1])
    tvalue = sheet[XLSX.encode_column_number(tcol)*"3"*":"*XLSX.encode_column_number(tcol)*"$stpr"][:,1]
    ZP=String[];ZH=String[];ZO=String[];Z=String[];
    for i in 1:length(tvalue)
        if tvalue[i]==0
            push!(ZP,fvalue[i])
            push!(Z,fvalue[i])
        elseif tvalue[i]==1
            push!(ZH,fvalue[i])
            push!(Z,fvalue[i])
        elseif tvalue[i]==2
            push!(ZO,fvalue[i])
            push!(Z,fvalue[i])
        end
    end
    setfield!(ms,:ZP,ZP)
    setfield!(ms,:ZH,ZH)
    setfield!(ms,:ZO,ZO)
    setfield!(ms,:Z,Z)
    setfield!(ps,:ZP,ZP)
    setfield!(ps,:ZH,ZH)
    setfield!(ps,:ZO,ZO)
    setfield!(ps,:Z,Z)
    map_zp=Dict{String,Int64}();map_zh=Dict{String,Int64}();map_zo=Dict{String,Int64}();
    for i in 1:length(ms.ZP)
        map_zp[ms.ZP[i]]=i
    end
    for i in 1:length(ms.ZH)
        map_zh[ms.ZH[i]]=i
    end
    for i in 1:length(ms.ZO)
        map_zo[ms.ZO[i]]=i
    end
    map_zpz=Dict{String,Int64}();map_zhz=Dict{String,Int64}();map_zoz=Dict{String,Int64}();
    for i in 1:length(ps.Z)
        if tvalue[i]==0
            map_zpz[ps.Z[i]]=i
        elseif tvalue[i]==1
            map_zhz[ps.Z[i]]=i
        elseif tvalue[i]==2
            map_zoz[ps.Z[i]]=i
        end
    end
    setfield!(mp,:map_zp,map_zp)
    setfield!(mp,:map_zh,map_zh)
    setfield!(mp,:map_zo,map_zo)
    setfield!(pp,:map_zp,map_zp)
    setfield!(pp,:map_zh,map_zh)
    setfield!(pp,:map_zo,map_zo)
    setfield!(pp,:map_zpz,map_zpz)
    setfield!(pp,:map_zhz,map_zhz)
    setfield!(pp,:map_zoz,map_zoz)

    ########### read investment ################
    sheet = XFI["investment"]
    ndict = Dict(:device=>"device",
                :ci=>"CaPeX",
                :cf=>"FixOM",
                :lf=>"Life",
                :type=>"type",
                :xm=>"Max_capacity",
                :location=>"location")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    Pnames = convert(Array{String,1},sheet[XLSX.encode_column_number(sttc)*"3"*":"*XLSX.encode_column_number(sttc)*"$stpr"][:,1])
    locationcol = findfirst(col_names .== ndict[:location])[2]
    location = convert(Array{Int64,1},sheet[XLSX.encode_column_number(locationcol)*"3"*":"*XLSX.encode_column_number(locationcol)*"$stpr"][:,1])
    typecol = findfirst(col_names .== ndict[:type])[2]
    type = sheet[XLSX.encode_column_number(typecol)*"3"*":"*XLSX.encode_column_number(typecol)*"$stpr"][:,1]

    L=String[];G=String[];RS=String[];RW=String[];SE=String[];EB=String[];Dcov=String[];W=String[];C=String[];GB=String[];Dsep=String[];
    PO=String[];PWI=String[];PWL=String[];E=String[];F=String[];SHy=String[];B=String[];GI=String[];
    for i in 1:length(Pnames)
        if type[i]=="line"
            push!(L,Pnames[i])
        elseif type[i]=="OCGT"
            push!(G,Pnames[i])
        elseif type[i]=="solar"
            push!(RS,Pnames[i])
        elseif type[i]=="wind"
            push!(RW,Pnames[i])
        elseif type[i]=="battery"
            push!(SE,Pnames[i])
        elseif type[i]=="E_boiler"
            push!(EB,Pnames[i])
        elseif type[i]=="converter"
            push!(Dcov,Pnames[i])
        elseif type[i]=="NG_well"
            push!(W,Pnames[i])
        elseif type[i]=="compressor"
            push!(C,Pnames[i])
        elseif type[i]=="NG_boiler"
            push!(GB,Pnames[i])
        elseif type[i]=="separator"
            push!(Dsep,Pnames[i])
        elseif type[i]=="oil_exp_pump"
            push!(PO,Pnames[i])
        elseif type[i]=="water_inj_pump"
            push!(PWI,Pnames[i])
        elseif type[i]=="water_lift_pump"
            push!(PWL,Pnames[i])
        elseif type[i]=="gas_inj_compressor"
            push!(GI,Pnames[i])
        elseif type[i]=="electrolyser"
            push!(E,Pnames[i])
        elseif type[i]=="fuel_cell"
            push!(F,Pnames[i])
        elseif type[i]=="hydrogen_store"
            push!(SHy,Pnames[i])
        elseif type[i]=="bus"
            push!(B,Pnames[i])
    end
    end
    setfield!(ps,:L,L);setfield!(ps,:G,G);setfield!(ps,:RS,RS);setfield!(ps,:RW,RW);setfield!(ps,:SE,SE);setfield!(ps,:EB,EB);
    setfield!(ps,:Dcov,Dcov);setfield!(ps,:W,W);setfield!(ps,:C,C);setfield!(ps,:GB,GB);setfield!(ps,:Dsep,Dsep);
    setfield!(ps,:PO,PO);setfield!(ps,:PWI,PWI);setfield!(ps,:PWL,PWL);setfield!(ps,:E,E);setfield!(ps,:F,F);
    setfield!(ps,:SHy,SHy);setfield!(ps,:B,B);setfield!(ps,:GI,GI);

    PP=String[];PH=String[];PO=String[];L=String[];P=String[];
    for l in 1:length(location)
        if location[l]==0
            push!(PP,Pnames[l])
            push!(P,Pnames[l])
        elseif location[l]==1
            push!(PH,Pnames[l])
            push!(P,Pnames[l])
        elseif location[l]==2
            push!(PO,Pnames[l])
            push!(P,Pnames[l])
        elseif location[l]==3
            push!(L,Pnames[l])
            push!(P,Pnames[l])
        end
    end
    setfield!(ms,:PP,PP)
    setfield!(ms,:PH,PH)
    setfield!(ms,:PO,PO)
    setfield!(ms,:L,L)
    setfield!(ms,:P,P)
    setfield!(ps,:P,P)
    map_ppp=Dict{String,Int64}();map_php=Dict{String,Int64}();map_pop=Dict{String,Int64}();map_ll=Dict{String,Int64}()
    for pp in 1:length(PP)
        map_ppp[PP[pp]]=pp
    end
    for ph in 1:length(PH)
        map_php[PH[ph]]=length(PP)+ph
    end
    for po in 1:length(PO)
        map_pop[PO[po]]=length(PP)+length(PH)+po
    end
    for l in 1:length(L)
        map_ll[L[l]]=length(PP)+length(PH)+length(PO)+l
    end
    setfield!(mp,:map_ppp,map_ppp)
    setfield!(mp,:map_php,map_php)
    setfield!(mp,:map_pop,map_pop)
    setfield!(mp,:map_ll,map_ll)
    cicol = findfirst(col_names .== ndict[:ci])[2]
    cival = convert(Array{Float64,1},sheet[XLSX.encode_column_number(cicol)*"3"*":"*XLSX.encode_column_number(cicol)*"$stpr"][:,1])
    cfcol = findfirst(col_names .== ndict[:cf])[2]
    cfval = convert(Array{Float64,1},sheet[XLSX.encode_column_number(cfcol)*"3"*":"*XLSX.encode_column_number(cfcol)*"$stpr"][:,1])
    xmcol = findfirst(col_names .== ndict[:xm])[2]
    xmval = convert(Array{Float64,1},sheet[XLSX.encode_column_number(xmcol)*"3"*":"*XLSX.encode_column_number(xmcol)*"$stpr"][:,1])
    lfcol = findfirst(col_names .== ndict[:lf])[2]
    lfval = convert(Array{Float64,1},sheet[XLSX.encode_column_number(lfcol)*"3"*":"*XLSX.encode_column_number(lfcol)*"$stpr"][:,1])
    c_inv_0 = cival./lfval

    ci = hcat(repeat(c_inv_0,1,NI1)).*10^3
    cf = cfval.*10^3
    xm0 = xmval.*10^-3
    setfield!(mp,:ci,ci)
    setfield!(mp,:cf,cf)
    setfield!(mp,:xm0,xm0)

    udict = Dict(:zn =>"node",
                :type=>"node_type",
                :x0 =>"capacity_at_0",
                :x5 =>"capacity_at_5")

    map_pp=Dict{String,Int64}();map_ph=Dict{String,Int64}();map_po=Dict{String,Int64}();map_l=Dict{String,Int64}()
    for pp in 1:length(PP)
        map_pp[PP[pp]]=pp
    end
    for ph in 1:length(PH)
        map_ph[PH[ph]]=ph
    end
    for po in 1:length(PO)
        map_po[PO[po]]=po
    end
    for l in 1:length(L)
        map_l[L[l]]=l
    end
    setfield!(mp,:map_pp,map_pp)
    setfield!(mp,:map_ph,map_ph)
    setfield!(mp,:map_po,map_po)
    setfield!(mp,:map_l,map_l)
    setfield!(pp,:map_l,map_l)

    xhpp=zeros(length(ms.PP),length(ms.ZP),NI2)
    for p in PP
        sheet = XFI[p]
        stt,stp = sheet.dimension.start,sheet.dimension.stop
        sttr,sttc = stt.row_number,stt.column_number
        stpr,stpc = stp.row_number,stp.column_number
        col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
        x5_col = findfirst(col_names .== udict[:x5])[2]
        x5_val = convert(Array{Float64,1},sheet[XLSX.encode_column_number(x5_col)*"3"*":"*XLSX.encode_column_number(x5_col)*"$stpr"][:,1])
        xhpp[map_pp[p],:,:] = hcat(repeat(x5_val,1,NI2)).*10^-3
    end
    setfield!(mp,:xhpp,xhpp)

    xhph=zeros(length(ms.PH),length(ms.ZH),NI2)
    for p in PH
        sheet = XFI[p]
        stt,stp = sheet.dimension.start,sheet.dimension.stop
        sttr,sttc = stt.row_number,stt.column_number
        stpr,stpc = stp.row_number,stp.column_number
        col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
        x5_col = findfirst(col_names .== udict[:x5])[2]
        x5_val = convert(Array{Float64,1},sheet[XLSX.encode_column_number(x5_col)*"3"*":"*XLSX.encode_column_number(x5_col)*"$stpr"][:,1])
        xhph[map_ph[p],:,:] = hcat(repeat(x5_val,1,NI2)).*10^-3
    end
    setfield!(mp,:xhph,xhph)

    xhpo=zeros(length(ms.PO),length(ms.ZO),NI2)
    for p in PO
        sheet = XFI[p]
        stt,stp = sheet.dimension.start,sheet.dimension.stop
        sttr,sttc = stt.row_number,stt.column_number
        stpr,stpc = stp.row_number,stp.column_number
        col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
        x5_col = findfirst(col_names .== udict[:x5])[2]
        x5_val = convert(Array{Float64,1},sheet[XLSX.encode_column_number(x5_col)*"3"*":"*XLSX.encode_column_number(x5_col)*"$stpr"][:,1])
        xhpo[map_po[p],:,:] = hcat(repeat(x5_val,1,NI2)).*10^-3
    end
    setfield!(mp,:xhpo,xhpo)
    xhl=zeros(length(ms.L),NI2)

    ############ read lines #################
    sheet = XFI["lines"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    x5_col = findfirst(col_names .== udict[:x5])[2]
    x5_val = convert(Array{Float64,1},sheet[XLSX.encode_column_number(x5_col)*"3"*":"*XLSX.encode_column_number(x5_col)*"$stpr"][:,1])
    xhl = hcat(repeat(x5_val,1,NI2)).*10^-3
    setfield!(mp,:xhl,xhl)
    setfield!(mp,:κ,5.)

    tx=Vector{Tuple{String,String}}()
    for i in 1:length(PP)
        for j in 1:length(ZP)
            push!(tx,(PP[i],ZP[j]))
        end
    end
    for i in 1:length(PH)
        for j in 1:length(ZH)
            push!(tx,(PH[i],ZH[j]))
        end
    end
    for i in 1:length(PO)
        for j in 1:length(ZO)
            push!(tx,(PO[i],ZO[j]))
        end
    end
    for l in 1:length(L)
        push!(tx,(L[l],L[l]))
    end
    for i in 1:length(th)
        push!(tx,(th[i],"-"))
    end
    setfield!(unc,:tx,tx)

    ############## gas turbine ############
    sheet = XFO_Electricity["gas_turbine"]
    ndict = Dict(:CG=>"VarOM",
                :CFuel=>"FuelCost",
                :EG=>"CO2_emissions",
                :ηG=>"efficiency",
                :GGR=>"ramping")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    Gnames = sheet[XLSX.encode_column_number(sttc)*"3"*":"*XLSX.encode_column_number(sttc)*"$stpr"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(Float64,sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1])
        setfield!(pp,vname,vvalue)
    end

    ############ battery ############
    sheet = XFO_Electricity["storage"]

    ndict = Dict(:ηSE=>"efficiency",
                :ΗSE=>"power_ratio",
                :CSE=>"VarOM")

    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    Bnames = sheet[XLSX.encode_column_number(sttc)*"3"*":"*XLSX.encode_column_number(sttc)*"$stpr"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(Float64,sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1])
        setfield!(pp,vname,vvalue)
    end

    # ############ time series ############
    sheet = XFO_Electricity["time_series"]

    ndict = Dict(1=>"start",
                2=>"end")

    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    Snames = sheet[XLSX.encode_column_number(sttc)*"3"*":"*XLSX.encode_column_number(sttc)*"$stpr"]
    setfield!(ps,:S,1:length(Snames))
    Hs = Array{Int64,2}(undef,maximum(ps.S),2)
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        Hs[:,vname] = convert(Array{Int64,1},sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][:,1])
    end
    setfield!(ps,:Hs,[Hs[s,1]:Hs[s,2] for s in ps.S])
    setfield!(ps,:H,1:8760)
    sH=Int64[];
    for s in ps.S
        for i in ps.Hs[s]
            push!(sH,i)
        end
    end
    setfield!(ps,:sH,sH)

    ############ renewable ############
    sheet = XFO_Electricity["wind"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    wvalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])
    setfield!(pp,:RPW,wvalue)
    sheet = XFO_Electricity["solar"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    svalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])
    setfield!(pp,:RPS,svalue)
    ############ onshore electricity price ############
    sheet = XFO_Electricity["electricity_price"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    pricevalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])
    setfield!(pp,:CZO,pricevalue)
    ############ transmission line efficiency #########
    sheet = XFO_Electricity["line_loss"]
    ndict = Dict(:ηL=>"loss_factor")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    Gnames = sheet[XLSX.encode_column_number(sttc)*"3"*":"*XLSX.encode_column_number(sttc)*"$stpr"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(Array{Float64,1},sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][:,1])
        setfield!(pp,vname,vvalue)
    end
    ############ bus line incidence matrix ############
    sheet = XFO_Electricity["bus_line_matrix"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    A = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+1)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])
    setfield!(pp,:AE,A)

    ############ other electricity ############
    sheet = XFO_Electricity["other"]
    ndict = Dict(:CLShed => "load_shedding",
                :CGShed => "generation_shedding",
                :WS  =>"alpha_adjustment",
                :Ht=>"hour",
                :μD=>"demand_scaling",
                :μE=>"CO2_limit_scaling",
                :CCO2=>"CO2_tax",
                :ΜE=>"CO2_limit",
                :σRes=>"spinning_reserve")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end
    ############ abatement cost ############
    sheet = XFO_Electricity["abatement_cost"]
    ndict = Dict(:CZOE => "abatement_cost_emission",
                :CZOP => "abatement_cost_power")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(Array{Float64,1},sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][:,1])
        setfield!(pp,vname,vvalue)
    end

    ############ natural gas ############
    sheet = XFO_NG["production"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    Dvalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])*0.84
    setfield!(pp,:VD,Dvalue)
    sheet = XFO_NG["gas_injection"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    Dvalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])*0.84
    setfield!(pp,:VGI,Dvalue)
    ############ gas well ############
    sheet = XFO_NG["well"]
    ndict = Dict(:CW =>"VarOM")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end
    ############ compressor ############
    sheet = XFO_NG["compressor"]
    ndict = Dict(:γC =>"compression ratio",
                 :ηC => "efficiency",
                 :α => "polytropic exponent",
                 :CC => "VarOM")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end
    ############ separator ############
    sheet = XFO_NG["separator"]
    ndict = Dict(:CSEP =>"VarOM",
                :κSEP=>"electricity_demand")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end
    ############ NG_boiler ############
    sheet = XFO_NG["NG_boiler"]
    ndict = Dict(:CGB =>"VarOM",
                 :ηGB =>"efficiency")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end
    ############ other ############
    sheet = XFO_NG["other"]
    ndict = Dict(:CNGShed =>"load_shed_cost",
                 :θNG => "gas_energy_content")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end
    ###### oil data ##########
    sheet = XFO_O["production"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    Dvalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])*900
    setfield!(pp,:VOD,Dvalue)
    sheet = XFO_O["oil_pump"]
    ndict = Dict(:κPO=>"electricity_demand_po")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end

    ###### water data ##########
    sheet = XFO_W["water_bore"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    Dvalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])*1000
    setfield!(pp,:VWB,Dvalue)
    sheet = XFO_W["water_injection"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    Dvalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])*1000
    setfield!(pp,:VWI,Dvalue)
    sheet = XFO_W["water_lift"]
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    Dvalue = convert(Array{Float64,2},sheet[XLSX.encode_column_number(sttc+1)*"$(sttr+2)"*":"*XLSX.encode_column_number(stpc)*"$stpr"])*1000
    setfield!(pp,:VWL,Dvalue)
    sheet = XFO_W["pump_wi"]
    ndict = Dict(:κPWI=>"electricity_demand_pwi")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end
    setfield!(pp,:VWL,Dvalue)
    sheet = XFO_W["pump_wl"]
    ndict = Dict(:κPWL=>"electricity_demand_pwl")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end

    ############ heat data ############
    sheet = XFO_heat["efficiency"]
    ndict = Dict(:ηGh =>"heat_recovery_gas_turbine",
                :ηFh  =>"heat_recovery_fuel_cell",
                :ρSEP=>"heat_demand_separator",
                :ηEB=>"electricity_boiler_efficiency",
                :CHLShed=>"heat_shed_cost")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end

    ############ electrolyser ############
    sheet = XFO_hydrogen["electrolyser"]
    ndict = Dict(:CE =>"VarOM",
                :ηES =>"efficiency_ES",
                :ηEF =>"efficiency_EF",
                :σE => "min_power")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end

    ############ hydrogen storage ############
    sheet = XFO_hydrogen["hydrogen_store"]
    ndict = Dict(:CSHy =>"VarOM",
                :σSHy =>"cushion_gas")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end

    ############ fuel cell ############
    sheet = XFO_hydrogen["fuel_cell"]
    ndict = Dict(:CF =>"VarOM",
                 :ηF => "efficiency",
                 :FR =>"ramping")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end

    ############ other hydrogen data ############
    sheet = XFO_hydrogen["other"]
    ndict = Dict(:θHY=>"hydrogen_energy_content")
    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]
    for vname in keys(ndict)
        vcol = findfirst(col_names .== ndict[vname])[2]
        vvalue = convert(typeof(getfield(pp,vname)),sheet[XLSX.encode_column_number(vcol)*"3"*":"*XLSX.encode_column_number(vcol)*"$stpr"][1,1])
        setfield!(pp,vname,vvalue)
    end

    ########## settings #########
    sheet = XFA["settings"]

    ndict = Dict(:Imax=>"Max_Iteration",
                 :tol =>"Convergence_Tol",
                 :alg =>"Type",
                 :γs  =>"stabilization_factor",
                 :w   =>"Subproblems_each_iter")

    stt,stp = sheet.dimension.start,sheet.dimension.stop
    sttr,sttc = stt.row_number,stt.column_number
    stpr,stpc = stp.row_number,stp.column_number
    col_names = sheet[XLSX.encode_column_number(sttc)*"1"*":"*XLSX.encode_column_number(stpc)*"1"]

    return unc,ms,mp,ps,pp
end


text:
"""
Basic function to read JEMRIS phantoms.
"""
function read_phantom_jemris(filename)
	# A(:,:,:,1)=Sample.M0;
	# I=find(Sample.T1); R1 =zeros(size(Sample.T1));  R1(I) =1./Sample.T1(I);
	# I=find(Sample.T2); R2 =zeros(size(Sample.T2));  R2(I) =1./Sample.T2(I);
	# I=find(Sample.T2S);R2S=zeros(size(Sample.T2S)); R2S(I)=1./Sample.T2S(I);
	# A(:,:,:,2)=R1;	#1/T1
	# A(:,:,:,3)=R2;	#1/T2
	# A(:,:,:,4)=R2S;	#1/T2s
	# A(:,:,:,5)=Sample.DB;
	fid = HDF5.h5open(filename)
	data = read(fid["sample/data"])
	Δx = read(fid["sample/resolution"]) * 1e-3 #[m]
	offset = read(fid["sample/offset"]) * 1e-3 #[m]
	mask = data[1,:,:,:] .!= 0
	#Maps
	ρ =   data[1,:,:,:]
	T1 =  1e-3 ./ data[2,:,:,:]
	T2 =  1e-3 ./ data[3,:,:,:]
	T2s = 1e-3 ./ data[4,:,:,:]
	Δw = data[5,:,:,:]
	#Positions
	X, Y, Z = size(ρ)
	FOVx = (X-1)*Δx[1] #[m]
	FOVy = (Y-1)*Δx[2] #[m]
	FOVz = (Z-1)*Δx[3] #[m]
	xx = [(-FOVx/2:Δx[1]:FOVx/2)...;]
	yy = [(-FOVy/2:Δx[2]:FOVy/2)...;;]
	zz = [(-FOVz/2:Δx[3]:FOVz/2)...;;;]
	x = xx*1 .+ yy*0 .+ zz*0 .+ offset[1]	#spin x coordinates
	y = xx*0 .+ yy*1 .+ zz*0 .+ offset[2]	#spin y coordinates
	z = xx*0 .+ yy*0 .+ zz*1 .+ offset[3]	#spin z coordinates
	v = 0 # m/s

	phantom = Phantom(name=basename(filename),
					x=    x[mask],
					y=    y[mask],
					z=    z[mask],
					ρ=    ρ[mask],
					T1=	 T1[mask],
					T2=	 T2[mask],
					T2s=T2s[mask],
					Δw=	 Δw[mask],
					ux=(x,y,z,t)->v*t
					)
	phantom
end


text:
gamma1arg = [
(:gamma, :(y.*digamma(x)), (-Inf,Inf)),
(:lfact, :(sign(y).*digamma(x+1)), (-Inf,Inf)),
(:lgamma, :(digamma(x)), (-Inf,Inf)),
(:digamma, :(trigamma(x)), (-Inf,Inf)), # polygamma(0,x)
(:trigamma, :(polygamma(2,x)), (-Inf,Inf)), # polygamma(1,x)
(:invdigamma, :(1./trigamma(y)), (-Inf,Inf)),
# zeta: TODO. Riemann 1-arg zeta
# eta # TODO. related to zeta
]

for (f,g,r) in gamma1arg
    @eval @primitive $f(x),dy,y  (dy.*($g))
    addtest1(f,r)
end

gamma2arg = [
#:beta => :TODO,                          # gamma,operators
#:lbeta => :TODO,                         # gamma,operators
#:zeta => :TODO,                 # Hurwitz 2-arg zeta
]

# polygamma wants x1 to be a non-negative integer, x2 unrestricted
@primitive polygamma(x1,x2),dy,y  nothing  unbroadcast(x2,dy.*polygamma(x1+1,x2))
addtest2(polygamma, 0:5, (-Inf,Inf))


text:
#
println("Ad) Tests towards a SCF field: B-spline primitives and one-particle spectra in a local potential.")
#

if  false
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true), nuclearModel=Nuclear.Model(18., "point"),
                            configs=[Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],  ##
                            asfSettings=AsfSettings(AsfSettings(), scField=Basics.DFSField())  )

    wb = perform(wa)

elseif  true
    grid = Radial.Grid(Radial.Grid(true), rnt = 2.0e-6,h = 3.0e-2, hp = 6.0e-2, NoPoints = 900)
    grid = Radial.Grid(true)
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(16., "point"),
                            configs=[Configuration("[Ne] 3s^2 3p^2")],  ## , Basics.DFSField()
                            asfSettings=AsfSettings(AsfSettings(), scField=Basics.DFSField())  )

    wb = perform(wa)

elseif false
    #
    # Test for Björn, 7. Mai 2020
    ## grid=JAC.Radial.Grid(true)
    grid = Radial.Grid(Radial.Grid(true), rnt = 2.0e-5,h = 5.0e-2, hp = 5.0e-2, NoPoints = 2000)
    wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(8.0, "Fermi"),
                            configs=[Configuration("1s^2 2s"), Configuration("1s 2s^2"), Configuration("1s 2s 3p"), Configuration("1s 2s 4p")],
                            ## configs=[Configuration("1s^2 2s"), Configuration("1s^2 2p"), Configuration("1s 2s 3s"), Configuration("1s 2p 3s"),
                            ##          Configuration("1s 2s 2p")],  ##
                            asfSettings=AsfSettings(true, false, Basics.DFSField(), "hydrogenic", Dict{Subshell, Orbital}(), [1],    40, 1.0e-6, JAC.Subshell[], JAC.Subshell[],
                                                    true, false, NoneQed(), LSjjSettings(false),
                                                    false, [1,2,3,4], false, JAC.LevelSymmetry[] )  )

    wb = perform(wa)

end


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

"""
    StandardUvNormal{T<:Real} <: Distributions.Distribution{Univariate,Continuous}

A standard normal distribution with a mean of zero and a variance of one.

Constructor:
```
    StandardUvNormal()
    StandardUvNormal{T<:Real}()
```
"""
struct StandardUvNormal{T<:Real} <: Distributions.Distribution{Univariate,Continuous}
end

StandardUvNormal() = StandardUvNormal{Float64}()

Distributions.Normal(d::StandardUvNormal{T}) where T = Normal{T}(zero(T), one(T))
Base.convert(::Type{Distributions.Normal}, d::StandardUvNormal) = Normal(d)

Base.minimum(d::StandardUvNormal{T}) where T = convert(float(T), -Inf)
Base.maximum(d::StandardUvNormal{T}) where T = convert(float(T), +Inf)

StatsBase.params(d::StandardUvNormal) = (mean(d), var(d))
@inline Distributions.partype(d::StandardUvNormal{T}) where {T<:Real} = T

Distributions.location(d::StandardUvNormal) = mean(d)
Distributions.scale(d::StandardUvNormal) = var(d)

Base.eltype(::Type{StandardUvNormal{T}}) where T = T

Statistics.mean(d::StandardUvNormal{T}) where T = zero(T)
StatsBase.median(d::StandardUvNormal{T}) where T = zero(T)
StatsBase.mode(d::StandardUvNormal{T}) where T = zero(T)
StatsBase.modes(d::StandardUvNormal{T}) where T = FillArrays.Zeros{T}(1)

Statistics.var(d::StandardUvNormal{T}) where T = one(T)
StatsBase.std(d::StandardUvNormal{T}) where T = one(T)
StatsBase.skewness(d::StandardUvNormal{T}) where {T<:Real} = zero(T)
StatsBase.kurtosis(d::StandardUvNormal{T}) where {T<:Real} = zero(T)

StatsBase.entropy(d::StandardUvNormal) = entropy(Normal(d))

Distributions.logpdf(d::StandardUvNormal, x::Real) = logpdf(Normal(d), x)
Distributions.pdf(d::StandardUvNormal, x::Real) = pdf(Normal(d), x)
Distributions.logcdf(d::StandardUvNormal, x::Real) = logcdf(Normal(d), x)
Distributions.cdf(d::StandardUvNormal, x::Real) = cdf(Normal(d), x)
Distributions.logccdf(d::StandardUvNormal, x::Real) = logccdf(Normal(d), x)
Distributions.ccdf(d::StandardUvNormal, x::Real) = ccdf(Normal(d), x)
Distributions.quantile(d::StandardUvNormal, p::Real) = quantile(Normal(d), p)
Distributions.cquantile(d::StandardUvNormal, p::Real) = cquantile(Normal(d), p)
Distributions.mgf(d::StandardUvNormal, t::Real) = mgf(Normal(d), t)
Distributions.cf(d::StandardUvNormal, t::Real) = cf(Normal(d), t)

Distributions.gradlogpdf(d::StandardUvNormal, x::Real) = gradlogpdf(Normal(d), x)

Base.rand(rng::AbstractRNG, d::StandardUvNormal{T}) where T = randn(rng, float(T))

Distributions.truncated(d::StandardUvNormal, l::Real, u::Real) = truncated(Normal(d), l, u)

function Distributions.product_distribution(dists::AbstractVector{StandardUvNormal{T}}) where T
    StandardMvNormal{T}(length(eachindex(dists)))
end

"""
    StandardMvNormal{T<:Real} <: Distributions.AbstractMvNormal

A standard `n`-dimensional multivariate normal distribution with it's mean at
the origin and an identity covariance matrix.

Constructor:
```
    StandardMvNormal(n::Integer)
    StandardMvNormal{T<:Real}(n::Integer)
```
"""
struct StandardMvNormal{T<:Real} <: Distributions.AbstractMvNormal
    _dim::Int
end

StandardMvNormal(dim::Integer) = StandardMvNormal{Float64}(dim)

Distributions.Product(d::StandardMvNormal{T}) where T = Distributions.Product(Fill(StandardUvNormal{T}(), length(d)))
Base.convert(::Type{Distributions.Product}, d::StandardMvNormal) = Distributions.Product(d)

Distributions.MvNormal(d::StandardMvNormal{T}) where T = MvNormal(ScalMat(length(d), one(T)))
Base.convert(::Type{Distributions.MvNormal}, d::StandardMvNormal) = MvNormal(d)

function Distributions.insupport(d::StandardMvNormal, x::AbstractVector)
    length(d) == length(eachindex(x)) && all(isfinite, x)
end

Base.eltype(::Type{StandardMvNormal{T}}) where T = T

Base.length(d::StandardMvNormal{T}) where T = d._dim

@inline function Base.view(d::StandardMvNormal{T}, i::Integer) where T
    Base.@boundscheck Base.checkindex(Bool, Base.OneTo(length(d)), i) || Base.throw_boundserror(d, i)
    StandardUvNormal{T}()
end

function Base.view(d::StandardMvNormal{T}, I::AbstractArray{<:Integer}) where T
    Base.@boundscheck Base.checkindex(Bool, Base.OneTo(length(d)), I) || Base.throw_boundserror(d, I)
    StandardMvNormal{T}(length(eachindex(I)))
end

StatsBase.params(d::StandardMvNormal) = (mean(d), cov(d))
@inline Distributions.partype(d::StandardMvNormal{T}) where {T<:Real} = T

Statistics.mean(d::StandardMvNormal{T}) where T = FillArrays.Zeros{T}(length(d))
Statistics.var(d::StandardMvNormal{T}) where T = FillArrays.Ones{T}(length(d))
Statistics.cov(d::StandardMvNormal{T}) where T = Diagonal(var(d))

StatsBase.mode(d::StandardMvNormal) = mean(d)
StatsBase.modes(d::StandardMvNormal) = Fill(mean(d), 1)

Distributions.invcov(d::StandardMvNormal{T}) where T = cov(d)
Distributions.logdetcov(d::StandardMvNormal{T}) where T = zero(T)

StatsBase.entropy(d::StandardMvNormal) = entropy(MvNormal(d))

Distributions._logpdf(d::StandardMvNormal, x::AbstractVector{<:Real}) = logpdf(MvNormal(d), x)
Distributions._pdf(d::StandardMvNormal, x::AbstractVector{<:Real}) = pdf(MvNormal(d), x)
Distributions.sqmahal(d::StandardMvNormal, x::AbstractVector{<:Real}) = sqmahal(MvNormal(d), x)
# Distributions.sqmahal!(r::AbstractVector{<:Real}, d::StandardMvNormal, x::AbstractMatrix{<:Real}) = sqmahal!(r, MvNormal(d), x)

Distributions.gradlogpdf(d::StandardMvNormal, x::AbstractVector{<:Real}) = Distributions.gradlogpdf(MvNormal(d), x)

function Distributions._rand!(rng::AbstractRNG, d::StandardMvNormal, A::AbstractVector{T}) where {T<:Real}
    broadcast!(x -> randn(rng, T), A, A)
end

function Distributions._rand!(rng::AbstractRNG, d::StandardMvNormal, A::AbstractMatrix{T}) where {T<:Real}
    broadcast!(x -> randn(rng, T), A, A)
end


text:
using Test, Downloads
using JSONStat

io = IOBuffer()

resp = Downloads.download("https://json-stat.org/samples/canada.json", io) |> take!
@test JSONStat.name(JSONStat.read(resp)) == "Population by sex and age group. Canada. 2012"

resp = Downloads.download("https://json-stat.org/samples/galicia.json", io) |> take!
@test JSONStat.Tables.columnnames(JSONStat.read(resp)) |> length == 7

max_stars_repo_name: lsaenzt/JSONStat

max_issues_repo_name: lsaenzt/JSONStat

max_forks_repo_name: lsaenzt/JSONStat


text:
using Documenter, Fermi

makedocs(sitename="Fermi Documentation",format = Documenter.HTML(prettyurls=false))
deploydocs(
           repo = "github.com/FermiQC/Fermi.jl.git",
          )


text:
using ComponentArrays, RecursiveArrayTools

function plotdemo()

    x = ComponentVector(a = ones(2), b = ComponentVector(b1 = fill(2, 5), b2 = zeros(3)))
    v = Vector{typeof(x)}()

    #make a Vector of ComponentVectors to simulate a simulation log
    for i in 1:10
        push!(v, copy(x))
    end

    #this returns a 10-element Vector of ComponentVector which can be indexed as
    #a matrix, while preserving the axis structure of its ComponentVector
    #elements
    v = VectorOfArray(v)

    #this produces an actual ComponentMatrix of size [length(x), 10]. the
    #awesome part is that its row dimension preserves the Axes metadata in the
    #original x. its column dimension has a Flat Axis.
    m = convert(Array, v)

    #we can now easily extract and plot the log of any specific block in the
    #original x
    log_b = m[:b, :]
    @show log_b2 = log_b[:b2, :]

end


text:
"""
# Refs
[1] T. Bian and Z.-P. Jiang,
“Value Iteration, Adaptive Dynamic Programming, and Optimal Control of Nonlinear Systems,”
in 2016 IEEE 55th Conference on Decision and Control (CDC), Las Vegas, NV, USA, Dec. 2016, pp. 3375–3380.
doi: 10.1109/CDC.2016.7798777.
[2] Forked repo of the code implementation by T. Bian, https://github.com/JinraeKim/nonlinear-VI/blob/main/nonlinear_sys.R
# Variables
V̂ ∈ R: the estimate of (state) value function
    - V̂.basis: Φ
    - V̂.param: w
dV̂ ∈ R: the estimate of time derivative of (state) value function
    - dV̂.basis: Ψ
    - V̂.param: c (BE CAREFUL; not used)
n: state dim.
m: input dim.
d: polynomial degree
"""
mutable struct CTValueIterationADP
    n::Int
    m::Int
    N_Ψ::Int
    V̂::LinearApproximator
    dV̂::LinearApproximator
    data
    ΣΦᵀΦ_inv
    Θ
    running_cost::Function
    u_norm_max::Real
    exact_integration::Bool
    function CTValueIterationADP(n::Int, m::Int, running_cost, u_norm_max,
                                 exact_integration, d_value::Int=2, d_controller::Int=4;
                                 V̂=LinearApproximator(:tao_bian_nonlinear_VI_V̂),
                                 dV̂=LinearApproximator(:tao_bian_nonlinear_VI_dV̂),
                                )
        @assert u_norm_max > 0.0
        # TODO: test with changing the function approximator, e.g., basis.
        data = DataFrame()
        ΣΦᵀΦ_inv = nothing
        Θ = nothing
        N_Ψ = length(dV̂.basis(rand(n+m)))
        new(n, m, N_Ψ, V̂, dV̂, data, ΣΦᵀΦ_inv, Θ, running_cost, u_norm_max, exact_integration)
    end
end

"""
# Notes
adp.data contains M+1 data
j = 0, 1, ..., M-1
"""
function set_data!(adp::CTValueIterationADP, data_new)
    adp.data = data_new  # write over
    ts = adp.data.time  # length: M+1
    xs = adp.data.sol |> Map(datum -> datum.state) |> collect  # length: M+1
    us = adp.data.sol |> Map(datum -> datum.input) |> collect  # length: M+1
    x_js = xs[1:end-1]  # length: M
    # ΣΦᵀΦ_inv (Eq. 16)
    Φs = xs |> Map(Φ(adp)) |> collect  # length: M+1
    Φ_js = Φs[1:end-1]  # length: M
    adp.ΣΦᵀΦ_inv = Φ_js |> Map(Φ -> Φ' * Φ) |> collect |> sum |> inv
    # Θ (Eq. 12)
    Θ_js = nothing
    # if adp.exact_integration
    if 1 == 0  # TODO it makes the algorithm diverge
        ∫Ψs = adp.data.sol |> Map(datum -> datum.∫Ψ) |> collect
        @show size(∫Ψs)
        Θ_js = diff(∫Ψs)
    else
        # numerical integration
        t_intervals = ts |> Partition(2; step=1) |> Map(copy) |> collect
        x_intervals = xs |> Partition(2; step=1) |> Map(copy) |> collect
        u_intervals = us |> Partition(2; step=1) |> Map(copy) |> collect
        Ψ_intervals = zip(xs, us) |> MapSplat(Ψ(adp)) |> Partition(2; step=1) |> Map(copy) |> collect
        Θ_js = integrate.(t_intervals, Ψ_intervals)
    end
    ΣΘᵀΘ_inv = Θ_js |> Map(Θ_j -> Θ_j' * Θ_j) |> collect |> sum |> inv
    Φ_diff_js = Φs |> diff  # length: M
    ΣΘᵀΦ_diff = zip(Θ_js, Φ_diff_js) |> MapSplat((Θ_j, Φ_diff_j) -> Θ_j' * Φ_diff_j) |> collect |> sum
    adp.Θ = ΣΘᵀΘ_inv * ΣΘᵀΦ_diff
    nothing
end

Φ(adp::CTValueIterationADP) = (x) -> reshape(adp.V̂.basis(x), 1, :)
Ψ(adp::CTValueIterationADP) = (x, u) -> reshape(adp.dV̂.basis(vcat(x, u)), 1, :)  # 1 × N_Ψ
function Ĥ(adp::CTValueIterationADP)
    r = adp.running_cost
    return (x, u, c) -> dot(Ψ(adp)(x, u), c) + r(x, u)
end

"""
Minimise approximate Hamiltonian.
"""
function min_Ĥ(adp::CTValueIterationADP)
    # TODO: Add inexact way of obtaining minimums;
    # see https://github.com/JinraeKim/nonlinear-VI/blob/main/nonlinear_sys.R#L98
    @unpack m, u_norm_max = adp
    return function (x, c)
        opt = NLopt.Opt(:LN_COBYLA, m)
        opt.min_objective = (u, grad) -> Ĥ(adp)(x, u, c)
        opt.xtol_rel = 1e-3
        u_norm_const = function (u, grad)
            norm(u) - u_norm_max
        end
        inequality_constraint!(opt, u_norm_const)
        (minf, minx, ret) = NLopt.optimize(opt, rand(m))
    end
end

function update!(adp::CTValueIterationADP, lr)
    @unpack m, ΣΦᵀΦ_inv, data, Θ = adp
    xs = adp.data.sol |> Map(datum -> datum.state) |> collect
    term2 = xs |> Map(x -> Φ(adp)(x)' * min_Ĥ(adp)(x, Θ*adp.V̂.param)[1]) |> collect |> sum
    adp.V̂.param += lr * ΣΦᵀΦ_inv * term2
    nothing
end

"""
Infer the approximate optimal input.
"""
function approximate_optimal_input(adp::CTValueIterationADP)
    # TODO: Add inexact way of obtaining approximate optimal input;
    # see https://github.com/JinraeKim/nonlinear-VI/blob/main/nonlinear_sys.R#L182
    ŵ_f = adp.V̂.param
    ĉ_f = adp.Θ * ŵ_f
    return function (x, p, t)
        û = min_Ĥ(adp)(x, ĉ_f)[2]
    end
end


text:
using Test

@testset "Beam" begin

    include(string(@__DIR__) * "/Beam/" * "BeamTest1.jl")
    include(string(@__DIR__) * "/Beam/" * "BeamTest2.jl")
    include(string(@__DIR__) * "/Beam/" * "BeamTest3.jl")

end


text:
"""
Description

        The foo() method is for making foo

Usage

        foo(12) = "foo"

Method Signature(s)

      foo(x::Int; y::UTF8String="foo")

Details

      The foo function makes foo

References

      http://www.duckduckgo.com

See Also

      bar(x) which makes bars

"""

function foo()
    # awesome code here
end

max_stars_repo_name: milktrader/Tuxedo

max_issues_repo_name: milktrader/Tuxedo

max_forks_repo_name: milktrader/Tuxedo