text:
"""
`module  JAC.PeriodicTable`
    ... a submodel of JAC that contains methods and data from the periodic table of elements. Although we here collect
        data for various elements, no attempt will be made to set-up a full compilation. Instead, useful information
        and some semi-empirical data are compiled to simplify the use of the JAC package in the future.
"""
module PeriodicTable

    using  Printf, ..Basics, ..ManyElectron

    """
    `PeriodicTable.getAtomicNumber(symbol::Symbol)`
        ... returns the atomic number (nuclear charge Z) for the given element symbol.
    """
    function getAtomicNumber(symbol::Symbol)
        for  (k,v) in PeriodicTable.data
            if  v.symbol == symbol    return( v.nuclearCharge )    end
        end

        error("Undefined element symbol in the periodic table.")
    end

    """
    `PeriodicTable.getData()`
        ... provides various data from the periodic table of elements.

    + `("mass", Z::Int64)`  or  `("mass", symbol::Symbol)`
        ... to get mean mass value of the element with nuclear charge Z or symbol.

    + `("density", Z::Int64)`  or  `("density", symbol::Symbol)`
        ... to get the density [g/cm^3] of the element with nuclear charge Z or symbol.

    + `("1st IP", Z::Int64)`  or  `("1st IP", symbol::Symbol)`
        ... to get the first ionization potential [eV] of the element with nuclear charge Z or symbol.

    + `("ground configuration", Z::Int64)`  or  `("ground configuration", symbol::Symbol)`
        ... to get the ground configuration of the neutral atom with nuclear charge Z or symbol.

    + `("atomic radius", Z::Int64)`  or  `("atomic radius", symbol::Symbol)`
        ... to get the atomic radius of the neutral atom with nuclear charge Z or symbol.

    + `("polarizibility", Z::Int64)`  or  `("polarizibility", symbol::Symbol)`
        ... to get the polarizibility of the neutral atom with nuclear charge Z or symbol.

    + `("ground configuration", Z::Int64)`  or  `("ground configuration", symbol::Symbol)`
        ... to get the ground configuration of the neutral atom with nuclear charge Z or symbol.

    + `("chemical hardness", Z::Int64)`  or  `("chemical hardness", symbol::Symbol)`
        ... to get the chemical hardness of the neutral atom with nuclear charge Z or symbol.
    """
    function getData(sa::String, zsymbol::Union{Int64,Symbol})
        if  typeof(zsymbol) == Symbol   Z = PeriodicTable.getAtomicNumber(zsymbol)    else    Z = zsymbol     end

        if        sa == "mass"              return ( PeriodicTable.data[Z].mass )
        elseif    sa == "density"           println("\nDensities are given [g/cm^3]: ...")
                                            return ( PeriodicTable.data[Z].density )
        elseif    sa == "1st IP"            println("\nIonization potentials are given [eV]: ...")
                                            return ( PeriodicTable.data[Z].firstIP )
        elseif    sa == "ground configuration"
                                            return ( Configuration(PeriodicTable.data[Z].configuration) )
        elseif    sa == "atomic radius"     println("\nAtomic radius is given [Angstroem]: ...")
                                            return ( PeriodicTable.data[Z].atomicRadius )
        elseif    sa == "polarizibility"    println("\nAtomic polarizibility is given [10^-24 c.c]: ...")
                                            return ( PeriodicTable.data[Z].polarizibility )
        elseif    sa == "chemical hardness" println("\nChemical hardness is given [eV]: ...")
                                            return ( PeriodicTable.data[Z].chemHardness )
        else      error("Unsupported keystring:: $sa")
        end
    end

    """
    `PeriodicTable.getIsotopeData()`
        ... provides various data from the periodic table of elements. Not yet implemented.

    + `("mass", Z::Int64)`  or  `("mass", symbol::Symbol)`
        ... to get mean mass value of the element with nuclear charge Z or symbol

    """
    function getIsotopeData(sa::String, zsymbol::Union{Int64,Symbol})

        if        sa in ["alpha", "fine-structure constant alpha"]          return (FINE_STRUCTURE_CONSTANT)
        else      error("Unsupported keystring:: $sa")
        end
    end

    ########################################################################################################################
    ########################################################################################################################
    ########################################################################################################################

    """
    `struct  PeriodicTable.Isotope`  ... defines a type for collecting data for an individual isotope.

        + massNumber         ::Float64           ... Mass number of the isotope.
        + isotopeMass        ::Float64           ... Mass of the isotope.
        + spin               ::AngularJ64        ... Nuclear spin of the isotope
        + mu                 ::Float64           ... Magnetic dipole moment of the isotope in nuclear magnetons.
        + Q                  ::Float64           ... Electric quadrupole moment of the isotope.
    """
    struct  Isotope
        massNumber           ::Float64
        isotopeMass          ::Float64
        spin                 ::AngularJ64
        mu                   ::Float64
        Q                    ::Float64
    end

    """
    `PeriodicTable.Isotope()`  ... simple constructur for an empty instance of an isotope.
    """
    function Isotope(i::Int64)
        Isotope( i, NaN, AngularJ64(0), NaN, NaN)
    end

    """
    `struct  PeriodicTable.Element`  ... defines a type for collecting data for each element.

        + nuclearCharge     ::Int64                 ... Nuclear charge
        + symbol            ::Symbol                ... Symbol of the element.
        + name              ::String                ... Name of the element.
        + mass              ::Float64               ... Mean mass number of the element.
        + density           ::Float64               ... Density of the element [g/cm^3].
        + firstIP           ::Float64               ... First ionization potential [eV].
        + configuration     ::String                ... Ground configuration of the element.
        + atomicRadius      ::Float64               ... Mean atomic radius of the neutral atom (Ghosh & Biswas, 2002).
        + polarizibility    ::Float64               ... Mean polarizibility of the neutral atom (Ghosh & Biswas, 2002).
        + chemHardness      ::Float64               ... Mean chemical hardness of the element (Ghosh & Biswas, 2002).
    """
    struct Element
        nuclearCharge       ::Int64
        symbol              ::Symbol
        name                ::String
        mass                ::Float64
        density             ::Float64
        firstIP             ::Float64
        configuration       ::String
        atomicRadius        ::Float64
        polarizibility      ::Float64
        chemHardness        ::Float64
    end

    """
    `PeriodicTable.Element(charge::Int64, symbol::Symbol)`  ... simple constructur for defining an element
    """
    function Element(charge::Int64, symbol::Symbol)
        Element(charge, symbol, "not defined", 999.99, Isotope[])
    end

    # `Base.show(io::IO, element::PeriodicTable.Element)`  ... prepares a proper printout of the variable element::PeriodicTable.Element.
    function Base.show(io::IO, element::PeriodicTable.Element)
        println(io, "Element  $(element.name) ($(element.symbol), Z=$(element.nuclearCharge), mass=$(element.mass))  ")
        println(io, "configuration:            $(element.configuration)  ")
    end

    """
    `PeriodicTable.data::Dict{Int64,PeriodicTable.Element}`
        ... provides all available data for an element with given nuclear charge Z::Int64 by PeriodicTable.data[Z];
	    cf. PeriodicTable.Element
    """
    const data = Dict{Int64,PeriodicTable.Element}(
                1 => Element(  1, :H,  "hydrogen",         1.008,   0.09,  13.60,  "1s",			            0.5292,     0.6669,   13.588    ),
                2 => Element(  2, :He, "helium",           4.003,   0.18,  24.59,  "1s^2",			            0.3113,     0.1358,   22.383    ),
                3 => Element(  3, :Li, "lithum",           6.941,   0.53,   5.39,  "[He] 2s",			        1.6282,    19.4239,    4.4164   ),
                4 => Element(  4, :Be, "beryllium",        9.012,   1.85,   9.32,  "[He] 2s^2", 		        1.0855,     5.7558,    6.6244   ),
                5 => Element(  5, :B,  "boron",           10.811,   2.34,   8.30,  "[He] 2s^2 2p",		        0.8141,     2.428 ,    8.8328   ),
                6 => Element(  6, :C,  "carbon",          12.011,   2.26,  11.26,  "[He] 2s^2 2p^2",		    0.6513,     1.2432,   11.0407   ),
                7 => Element(  7, :N,  "nitrogen",        14.007,   1.25,  14.53,  "[He] 2s^2 2p^3",		    0.5427,     0.7193,   13.25     ),
                8 => Element(  8, :O,  "oxygen",          15.999,   1.43,  13.62,  "[He] 2s^2 2p^4",		    0.4652,     0.453 ,   15.4574   ),
                9 => Element(  9, :F,  "fluorine",        18.998,   1.70,  17.42,  "[He] 2s^2 2p^5",		    0.4071,     0.3034,   17.6634   ),
               10 => Element( 10, :Ne, "neon",            20.180,   0.90,  21.56,  "[Ne]",			            0.3618,     0.2131,   19.875    ),
               11 => Element( 11, :Na, "sodium",          22.990,   0.97,   5.14,  "[Ne] 3s",			        2.1649,    45.659 ,    3.3215   ),
               12 => Element( 12, :Mg, "magnesium",       24.305,   1.74,   7.65,  "[Ne] 3s^2", 		        1.6711,    21.0   ,    4.303    ),
               13 => Element( 13, :Al, "aluminium",       26.982,   2.70,   5.99,  "[Ne] 3s^2 3p",		        1.3607,    11.337 ,    5.2846   ),
               14 => Element( 14, :Si, "silicon",         28.086,   2.33,   8.15,  "[Ne] 3s^2 3p^2",		    1.1476,     6.8011,    6.2659   ),
               15 => Element( 15, :P,  "phosphorus",      30.974,   1.82,  10.49,  "[Ne] 3s^2 3p^3",		    0.9922,     4.3955,    7.2473   ),
               16 => Element( 16, :S,  "sulfur",          32.065,   2.07,  10.36,  "[Ne] 3s^2 3p^4",		    0.8738,     3.0023,    8.2293   ),
               17 => Element( 17, :Cl, "chlorine",        35.453,   3.21,  12.97,  "[Ne] 3s^2 3p^5",		    0.7807,     2.1412,    9.2107   ),
               18 => Element( 18, :Ar, "argon",           39.948,   1.78,  15.76,  "[Ar]",			            0.7056,     1.5808,    10.191   ),
               19 => Element( 19, :K,  "potassium",       39.098,   0.86,   4.34,  "[Ar] 4s",			        3.5598,   202.9969,     2.02    ),
               20 => Element( 20, :Ca, "calcium",         40.078,   1.55,   6.11,  "[Ar] 4s^2", 		        2.7479,    93.3717,     2.6162  ),
               21 => Element( 21, :Sc, "scandium",        44.956,   2.99,   6.56,  "[Ar] 3d 4s^2",		        2.6106,    80.0633,     2.7545  ),
               22 => Element( 22, :Ti, "titanium",        47.867,   4.54,   6.83,  "[Ar] 3d^2 4s^2",		    2.4861,    69.1712,     2.8924  ),
               23 => Element( 23, :V,  "vanadium",        50.942,   6.11,   6.75,  "[Ar] 3d^3 4s^2",		    2.3732,    60.1472,     3.03    ),
               24 => Element( 24, :Cr, "chromium",        51.996,   7.19,   6.77,  "[Ar] 3d^5 4s",		        2.2701,    52.6438,     3.1676  ),
               25 => Element( 25, :Mn, "manganese",       54.938,   7.43,   7.43,  "[Ar] 3d^5 4s^2",		    2.1754,    46.3265,     3.3055  ),
               26 => Element( 26, :Fe, "iron",            55.845,   7.87,   7.90,  "[Ar] 3d^6 4s^2",		    2.0885,    40.9939,     3.443   ),
               27 => Element( 27, :Co, "cobalt",          58.933,   8.90,   7.88,  "[Ar] 3d^7 4s^2",		    2.008 ,    36.4555,     3.5811  ),
               28 => Element( 28, :Ni, "nickel",          58.693,   8.90,   7.64,  "[Ar] 3d^8 4s^2",		    1.9337,    32.5372,     3.7187  ),
               29 => Element( 29, :Cu, "copper",          63.546,   8.96,   7.73,  "[Ar] 3d^10 4s",		        1.8648,    29.1769,     3.8561  ),
               30 => Element( 30, :Zn, "zinc",            65.390,   7.13,   9.39,  "[Ar] 3d^10 4s^2",		    1.8004,    26.2615,     3.994   ),
               31 => Element( 31, :Ga, "gallium",         69.723,   5.91,   6.00,  "[Ar] 3d^10 4s^2 4p",	    1.5663,    17.2917,     4.5909  ),
               32 => Element( 32, :Ge, "germanium",       72.640,   5.32,   7.90,  "[Ar] 3d^10 4s^2 4p^2",	    1.3862,    11.9864,     5.1874  ),
               33 => Element( 33, :As, "arsenic",         74.922,   5.72,   9.79,  "[Ar] 3d^10 4s^2 4p^3",	    1.2431,     8.6443,     5.7846  ),
               34 => Element( 34, :Se, "selenium",        78.960,   4.79,   9.75,  "[Ar] 3d^10 4s^2 4p^4",	    1.1269,     6.4397,     6.381   ),
               35 => Element( 35, :Br, "bromine",         79.904,   3.12,  11.81,  "[Ar] 3d^10 4s^2 4p^5",	    1.0305,     4.9244,     6.978   ),
               36 => Element( 36, :Kr, "krypton",         83.800,   3.75,  14.00,  "[Kr]",			            0.9493,     3.8497,     7.5748  ),
               37 => Element( 37, :Rb, "rubidium",        85.468,   1.63,   4.18,  "[Kr] 5s",			        4.8106,   500.9683,     1.4948  ),
               38 => Element( 38, :Sr, "strontium",       87.620,   2.54,   5.69,  "[Kr] 5s^2", 		        3.7135,   230.4462,     1.9364  ),
               39 => Element( 39, :Y,  "yttrium",         88.906,   4.47,   6.22,  "[Kr] 4d 5s^2",		        3.5278,   197.5715,     2.0383  ),
               40 => Element( 40, :Zr, "zirconium",       91.224,   6.51,   6.63,  "[Kr] 4d^2 5s^2",		    3.3598,   170.6683,     2.1402  ),
               41 => Element( 41, :Nb, "niobium",         92.906,   8.57,   6.76,  "[Kr] 4d^4 5s",		        3.2071,   148.4397,     2.2422  ),
               42 => Element( 42, :Mo, "molybdenum",      95.940,  10.22,   7.09,  "[Kr] 4d^5 5s",		        3.0677,   129.9126,     2.344   ),
               43 => Element( 43, :Tc, "technetium",      98.000,  11.50,   7.28,  "[Kr] 4d^5 5s^2",		    2.9398,   114.3315,     2.446   ),
               44 => Element( 44, :Ru, "ruthenium",      101.070,  12.37,   7.36,  "[Kr] 4d^7 5s",		        2.8222,   101.1523,     2.5479  ),
               45 => Element( 45, :Rh, "rhodium",        102.906,  12.41,   7.46,  "[Kr] 4d^8 5s",		        2.7137,    89.9286,     2.6498  ),
               46 => Element( 46, :Pd, "palladium",      106.420,  12.02,   8.34,  "[Kr] 4d^10",		        2.6132,    80.3028,     2.7517  ),
               47 => Element( 47, :Ag, "silver",         107.868,  10.50,   7.58,  "[Kr] 4d^10 5s",		        2.5199,    72.005 ,     2.8536  ),
               48 => Element( 48, :Cd, "cadmium",        112.411,   8.65,   8.99,  "[Kr] 4d^10 5s^2",		    2.433 ,    64.8095,     2.9555  ),
               49 => Element( 49, :In, "indium",         114.818,   7.31,   5.79,  "[Kr] 4d^10 5s^2 5p",	    2.1167,    42.6767,     3.3972  ),
               50 => Element( 50, :Sn, "tin",            118.710,   7.31,   7.34,  "[Kr] 4d^10 5s^2 5p^2",	    1.8732,    29.5777,     3.8388  ),
               51 => Element( 51, :Sb, "antimony",       121.760,   6.68,   8.61,  "[Kr] 4d^10 5s^2 5p^3",	    1.6799,    21.3335,     4.2805  ),
               52 => Element( 52, :Te, "tellurium",      127.600,   6.24,   9.01,  "[Kr] 4d^10 5s^2 5p^4",	    1.5228,    15.8906,     4.7221  ),
               53 => Element( 53, :I,  "iodine",         126.905,   4.93,  10.45,  "[Kr] 4d^10 5s^2 5p^5",	    1.3926,    12.1532,     5.0636  ),
               54 => Element( 54, :Xe, "xenon",          131.293,   5.90,  12.13,  "[Xe]",			            1.2828,     9.4993,     5.6056  ),
               55 => Element( 55, :Cs, "caesium",        132.906,   1.87,   3.89,  "[Xe] 6s",			        6.0615,  1002.1964,     1.1863  ),
               56 => Element( 56, :Ba, "barium",         137.327,   3.59,   5.21,  "[Xe] 6s^2", 		        4.6788,   460.9098,     1.5369  ),
               57 => Element( 57, :La, "lanthanum",      138.906,   6.15,   5.58,  "[Xe] 5d 6s^2",		        3.8102,   248.9177,     1.8873  ),
               58 => Element( 58, :Ce, "cerium",         140.116,   6.77,   5.54,  "[Xe] 4f 5d 6s^2",		    3.2133,   149.2883,     2.2378  ),
               59 => Element( 59, :Pr, "praseodymium",   140.908,   6.77,   5.47,  "[Xe] 4f^3 6s^2",		    2.778 ,    96.4738,     2.5885  ),
               60 => Element( 60, :Nd, "neodymium",      144.240,   7.01,   5.53,  "[Xe] 4f^4 6s^2",		    2.4468,    65.9186,     2.9389  ),
               61 => Element( 61, :Pm, "promethium",     145.000,   7.30,   5.58,  "[Xe] 4f^5 6s^2",		    2.1861,    47.0135,     3.2893  ),
               62 => Element( 62, :Sm, "samarium",       150.360,   7.52,   5.64,  "[Xe] 4f^6 6s^2",		    1.9756,    34.6984,     3.6398  ),
               63 => Element( 63, :Eu, "europium",       151.964,   5.24,   5.67,  "[Xe] 4f^7 6s^2",		    1.802 ,    26.3316,     3.9905  ),
               64 => Element( 64, :Gd, "gadolinium",     157.250,   7.90,   6.15,  "[Xe] 4f^7 5d 6s^2", 	    1.6565,    20.4544,     4.341   ),
               65 => Element( 65, :Tb, "terbium",        158.925,   8.23,   5.86,  "[Xe] 4f^9 6s^2",		    1.5328,    16.2057,     4.6913  ),
               66 => Element( 66, :Dy, "dysprosium",     162.500,   8.55,   5.94,  "[Xe] 4f^10 6s^2",		    1.4262,    13.0543,     5.0419  ),
               67 => Element( 67, :Ho, "holmium",        164.930,   8.80,   6.02,  "[Xe] 4f^11 6s^2",		    1.3335,    10.6707,     5.3924  ),
               68 => Element( 68, :Er, "erbium",         167.259,   9.07,   6.11,  "[Xe] 4f^12 6s^2",		    1.2521,     8.8334,     5.743   ),
               69 => Element( 69, :Tm, "thulium",        168.934,   9.32,   6.18,  "[Xe] 4f^13 6s^2",		    1.1801,     7.3955,     6.0934  ),
               70 => Element( 70, :Yb, "ytterbium",      173.040,   6.90,   6.25,  "[Xe] 4f^14 6s^2",		    1.1159,     6.253 ,     6.4439  ),
               71 => Element( 71, :Lu, "lutetium",       174.967,   9.84,   5.43,  "[Xe] 4f^14 5d 6s^2",	    1.0583,     5.3338,     6.7947  ),
               72 => Element( 72, :Hf, "hafnium",        178.490,  13.31,   6.83,  "[Xe] 4f^14 5d^2 6s^2",	    1.0079,     4.6075,     7.1344  ),
               73 => Element( 73, :Ta, "tantalum",       180.948,  16.65,   7.55,  "[Xe] 4f^14 5d^3 6s^2",	    0.9594,     3.9739,     7.4951  ),
               74 => Element( 74, :W,  "tungsten",       183.840,  19.35,   7.86,  "[Xe] 4f^14 5d^4 6s^2",	    0.9165,     3.4643,     7.8459  ),
               75 => Element( 75, :Re, "rhenium",        186.207,  21.04,   7.83,  "[Xe] 4f^14 5d^5 6s^2",	    0.8773,     3.0385,     8.1965  ),
               76 => Element( 76, :Os, "osmium",         190.230,  22.60,   8.44,  "[Xe] 4f^14 5d^6 6s^2",	    0.8413,     2.6796,     8.5472  ),
               77 => Element( 77, :Ir, "iridium",        192.217,  22.40,   8.97,  "[Xe] 4f^14 5d^7 6s^2",	    0.8182,     2.3756,     8.8973  ),
               78 => Element( 78, :Pt, "platinum",       195.078,  21.45,   8.96,  "[Xe] 4f^14 5d^9 6s",	    0.7776,     2.1175,     9.2474  ),
               79 => Element( 79, :Au, "gold",           196.967,  19.32,   9.23,  "[Xe] 4f^14 5d^10 6s",	    0.7492,     1.8924,     9.598   ),
               80 => Element( 80, :Hg, "mercury",        200.590,  13.55,  10.44,  "[Xe] 4f^14 5d^10 6s^2",	    3.0636,   129.646 ,     2.3456  ),
               81 => Element( 81, :Tl, "thallium",       204.383,  11.85,   6.11,  "[Xe] 4f^14 5d^10 6s^2 6p",  2.667 ,    85.3653,     2.6962  ),
               82 => Element( 82, :Pb, "lead",           207.200,  11.35,   7.42,  "[Xe] 4f^14 5d^10 6s^2 6p^2",2.3603,    59.1717,     3.0466  ),
               83 => Element( 83, :Bi, "bismuth",        208.980,   9.75,   7.29,  "[Xe] 4f^14 5d^10 6s^2 6p^3",2.1167,    42.6767,     3.3972  ),
               84 => Element( 84, :Po, "polonium",       209.000,   9.30,   8.42,  "[Xe] 4f^14 5d^10 6s^2 6p^4",1.9187,    31.7858,     3.7477  ),
               85 => Element( 85, :At, "astatine",       210.000,    NaN,   9.30,  "[Xe] 4f^14 5d^10 6s^2 6p^5",1.7546,    24.3079,     4.0983  ),
               86 => Element( 86, :Rn, "radon",          222.000,   9.73,  10.75,  "[Rn]",			            1.6164,    19.0046,     4.4487  ),
               87 => Element( 87, :Fr, "francium",       223.000,    NaN,   4.07,  "[Rn] 7s",			        7.2404,  1705.0   ,     0.9913  ),
               88 => Element( 88, :Ra, "radium",         226.000,   5.50,   5.28,  "[Rn] 7s^2", 		        5.5887,   785.4977,     1.2867  ),
               89 => Element( 89, :Ac, "actinium",       227.000,  10.07,   5.17,  "[Rn] 6d 7s^2",		        5.3091,   673.4033,     1.3544  ),
               90 => Element( 90, :Th, "thorium",        232.038,  11.72,   6.31,  "[Rn] 6d^2 7s^2",		    5.0569,   581.6815,     1.4222  ),
               91 => Element( 91, :Pa, "protactinium",   231.036,  15.40,   5.89,  "[Rn] 5f^2 6d 7s^2", 	    3.7042,   228.7156,     1.9413  ),
               92 => Element( 92, :U,  "uranium",        238.029,  18.95,   6.19,  "[Rn] 5f^3 6d 7s^2", 	    3.2177,   149.9164,     2.2348  ),
               93 => Element( 93, :Np, "neptunium",      237.000,  20.20,   6.27,  "[Rn] 5f^4 6d 7s^2", 	    2.8443,   103.5473,     2.5281  ),
               94 => Element( 94, :Pu, "plutonium",      244.000,  19.84,   6.03,  "[Rn] 5f^6 7s^2",		    2.3596,    59.1266,     3.0473  ),
               95 => Element( 95, :Am, "americium",      243.000,  13.67,   5.97,  "[Rn] 5f^7 7s^2",		    2.1525,    44.8789,     3.3407  ),
               96 => Element( 96, :Cm, "curium",         247.000,  13.50,   5.99,  "[Rn] 5f^7 6d 7s^2", 	    2.1097,    42.2547,     3.4084  ),
               97 => Element( 97, :Bk, "berkelium",      247.000,  14.78,   6.20,  "[Rn] 5f^9 7s^2",		    1.8308,    27.5918,     3.9277  ),
               98 => Element( 98, :Cf, "californium",    251.000,  15.10,   6.28,  "[Rn] 5f^10 7s^2",		    1.7035,    22.2453,     4.2212  ),
               99 => Element( 99, :Es, "einsteinium",    252.000,    NaN,   6.42,  "[Rn] 5f^11 7s^2",		    1.5928,    18.1843,     4.5146  ),
              100 => Element(100, :Fm, "fermium",        257.000,    NaN,   6.50,  "[Rn] 5f^12 7s^2",		    1.4956,    15.0542,     4.808   ),
              101 => Element(101, :Md, "mendelevium",    258.000,    NaN,   6.58,  "[Rn] 5f^13 7s^2",		    1.4096,    12.6038,     5.1013  ),
              102 => Element(102, :No, "nobelium",       259.000,    NaN,   6.65,  "[Rn] 5f^14 7s^2",		    1.3329,    10.6563,     5.3949  ),
              103 => Element(103, :Lr, "lawrencium",     262.000,    NaN,   4.90,  "[Rn] 5f^14 7s^2 7p",	    1.3164,    10.2631,     5.4629  ),
              104 => Element(104, :Rf, "rutherfordium",  261.000,    NaN,    NaN,  "[Rn] 5f^14 6d^2 7s^2",	       NaN,        NaN,        NaN  ),
              105 => Element(105, :Db, "dubnium",        262.000,    NaN,    NaN,  "[Rn] 5f^14 6d^3 7s^2",	       NaN,        NaN,        NaN  ),
              106 => Element(106, :Sg, "seaborgium",     266.000,    NaN,    NaN,  "[Rn] 5f^14 6d^4 7s^2",	       NaN,        NaN,        NaN  ),
              107 => Element(107, :Bh, "bohrium",        264.000,    NaN,    NaN,  "[Rn] 5f^14 6d^5 7s^2",	       NaN,        NaN,        NaN  ),
              108 => Element(108, :Hs, "hassium",        277.000,    NaN,    NaN,  "[Rn] 5f^14 6d^6 7s^2",	       NaN,        NaN,        NaN  ),
              109 => Element(109, :Mt, "meitnerium",     268.000,    NaN,    NaN,  "[Rn] 5f^14 6d^7 7s^2",	       NaN,        NaN,        NaN  ),
              110 => Element(110, :Ds, "darmstadtium",       NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^8 7s^2",	       NaN,        NaN,        NaN  ),
              111 => Element(111, :Rg, "roentgenium",        NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^9 7s^2",	       NaN,        NaN,        NaN  ),
              112 => Element(112, :Cn, "copernicium",        NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^10 7s^2",	       NaN,        NaN,        NaN  ),
              113 => Element(113, :Nh, "nihonium",           NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^8 7s^2 7p",	   NaN,        NaN,        NaN  ),
              114 => Element(114, :Fl, "flerovium",          NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^2 7s^2 7p^2",    NaN,        NaN,        NaN  ),
              115 => Element(115, :Mc, "moscovium",          NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^3 7s^2 7p^3",    NaN,        NaN,        NaN  ),
              116 => Element(116, :Lv, "livermorium",        NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^4 7s^2 7p^4",    NaN,        NaN,        NaN  ),
              117 => Element(117, :Ts, "tennessine",         NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^5 7s^2 7p^5",    NaN,        NaN,        NaN  ),
              118 => Element(118, :Og, "oganesson",          NaN,    NaN,    NaN,  "[Rn] 5f^14 6d^6 7s^2 7p^6",    NaN,        NaN,        NaN  )	 )

    """
    `PeriodicTable.isotopes::Dict{Int64,Array{PeriodicTable.Isotope,1}}`
        ... provides isotope data for an element with given nuclear charge Z::Int64 from PeriodicTable.isotopes[Z];
	        cf. PeriodicTable.Isotope
    """
    const isotopes = Dict{Int64,Array{PeriodicTable.Isotope,1}}(
                1 => Isotope[],
                2 => Isotope[ Isotope( 2, 2.015894,      AngularJ64(0), NaN, NaN),    Isotope( 2, 3.01602932265, AngularJ64(1//2), NaN, NaN),
                              Isotope( 4, 4.00260325413, AngularJ64(0), NaN, NaN)],
                3 => Isotope[],
                4 => Isotope[],
                5 => Isotope[],
                6 => Isotope[],
                7 => Isotope[],
                8 => Isotope[],
                9 => Isotope[],
               10 => Isotope[],
               11 => Isotope[],
               12 => Isotope[],
               13 => Isotope[],
               14 => Isotope[],
               15 => Isotope[],
               16 => Isotope[],
               17 => Isotope[],
               18 => Isotope[],
               19 => Isotope[],
               20 => Isotope[],
               21 => Isotope[],
               22 => Isotope[],
               23 => Isotope[],
               24 => Isotope[],
               25 => Isotope[],
               26 => Isotope[],
               27 => Isotope[],
               28 => Isotope[],
               29 => Isotope[],
               30 => Isotope[],
               31 => Isotope[],
               32 => Isotope[],
               33 => Isotope[],
               34 => Isotope[],
               35 => Isotope[],
               36 => Isotope[],
               37 => Isotope[],
               38 => Isotope[],
               39 => Isotope[],
               40 => Isotope[],
               41 => Isotope[],
               42 => Isotope[],
               43 => Isotope[],
               44 => Isotope[],
               45 => Isotope[],
               46 => Isotope[],
               47 => Isotope[],
               48 => Isotope[],
               49 => Isotope[],
               50 => Isotope[],
               51 => Isotope[],
               52 => Isotope[],
               53 => Isotope[],
               54 => Isotope[],
               55 => Isotope[],
               56 => Isotope[],
               57 => Isotope[],
               58 => Isotope[],
               59 => Isotope[],
               60 => Isotope[],
               61 => Isotope[],
               62 => Isotope[],
               63 => Isotope[],
               64 => Isotope[],
               65 => Isotope[],
               66 => Isotope[],
               67 => Isotope[],
               68 => Isotope[],
               69 => Isotope[],
               70 => Isotope[],
               71 => Isotope[],
               72 => Isotope[],
               73 => Isotope[],
               74 => Isotope[],
               75 => Isotope[],
               76 => Isotope[],
               77 => Isotope[],
               78 => Isotope[],
               79 => Isotope[],
               80 => Isotope[],
               81 => Isotope[],
               82 => Isotope[],
               83 => Isotope[],
               84 => Isotope[],
               85 => Isotope[],
               86 => Isotope[],
               87 => Isotope[],
               88 => Isotope[],
               89 => Isotope[],
               90 => Isotope[],
               91 => Isotope[],
               92 => Isotope[],
               93 => Isotope[],
               94 => Isotope[],
               95 => Isotope[],
               96 => Isotope[],
               97 => Isotope[],
               98 => Isotope[],
               99 => Isotope[],
              100 => Isotope[],
              101 => Isotope[],
              102 => Isotope[],
              103 => Isotope[],
              104 => Isotope[],
              105 => Isotope[],
              106 => Isotope[],
              107 => Isotope[],
              108 => Isotope[],
              109 => Isotope[],
              110 => Isotope[],
              111 => Isotope[],
              112 => Isotope[],
              113 => Isotope[],
              114 => Isotope[],
              115 => Isotope[],
              116 => Isotope[],
              117 => Isotope[],
              118 => Isotope[] )

    """
    `PeriodicTable.bindingEnergies_Williams2000(Z::Int64)`  ... former  JAC.store_Williams2000(Z::Int64)
        ... to return the `stored values' of Williams et al. (2000), ... for the element with nuclear charge Z.
    """
    function bindingEnergies_Williams2000(Z::Int64)
        # This procedure 'stores' the binding energies of the inner-shell electrons 1s, ..., 3d for all elements from Z = 1, ..., 54;
        # all binding energies are given in eV
        #
        #                        1s_1/2,   2s_1/2,   2p_1/2,   2p_3/2,   3s_1/2,   3p_1/2,   3p_3/2,   3d_3/2,   3d_5/2,   4s_1/2,   4p_1/2,   4p_3/2
        if     Z ==  1   wa = [    13.6,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z ==  2   wa = [    24.6,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z ==  3   wa = [    54.7,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z ==  4   wa = [   111.5,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z ==  5   wa = [   188.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z ==  6   wa = [   284.2,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z ==  7   wa = [   409.9,     37.3,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z ==  8   wa = [   543.1,     41.6,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z ==  9   wa = [   696.7,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 10   wa = [   870.2,     48.5,     21.7,     21.6,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        #
        elseif Z == 11   wa = [  1070.8,     63.5,    30.65,    30.81,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 12   wa = [  1303.0,     88.7,    49.78,    49.50,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 13   wa = [  1559.6,    117.8,    72.95,    72.55,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 14   wa = [  1839.0,    149.7,    99.82,    99.42,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 15   wa = [  2145.5,    189.0,   136.0,    135.0,      -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 16   wa = [  2472.0,    230.9,   163.6,    162.5,      -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 17   wa = [  2822.4,    270.0,   202.0,    200.0,      -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 18   wa = [  3205.9,    326.3,   250.6,    248.4,      29.3,     15.9,     15.7,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 19   wa = [  3608.4,    378.6,   297.3,    294.6,      34.8,     18.3,     18.3,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 20   wa = [  4038.5,    438.4,   349.7,    346.2,      44.3,     25.4,     25.4,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        #
        elseif Z == 21   wa = [  4492.0,    498.0,   403.6,    398.7,      51.1,     28.3,     28.3,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 22   wa = [  4966.0,    560.9,   460.2,    453.8,      58.7,     32.6,     32.6,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 23   wa = [  5465.0,    626.7,   519.8,    512.1,      66.3,     37.2,     37.2,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 24   wa = [  5989.0,    696.0,   583.8,    574.1,      74.1,     42.2,     42.2,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 25   wa = [  6539.0,    769.1,   649.9,    638.7,      82.3,     47.2,     47.2,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 26   wa = [  7112.0,    844.6,   719.9,    706.8,      91.3,     52.7,     52.7,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 27   wa = [  7709.0,    925.1,   793.2,    778.1,     101.0,     58.9,     59.9,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 28   wa = [  8333.0,   1008.6,   870.0,    852.7,     110.8,     68.0,     66.2,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 29   wa = [  8979.0,   1096.7,   952.3,    932.7,     122.5,     77.3,     75.1,     -1.0,     -1.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 30   wa = [  9659.0,   1196.2,  1044.9,   1021.8,     139.8,     91.4,     88.6,     10.2,     10.1,     -1.0,     -1.0,     -1.0]
        #
        elseif Z == 31   wa = [ 10367.0,   1299.0,  1143.2,   1116.4,     159.5,    103.5,    100.0,     18.7,     18.7,     -1.0,     -1.0,     -1.0]
        elseif Z == 32   wa = [ 11103.0,   1414.6,  1248.1,   1217.0,     180.1,    124.9,    120.8,     29.8,     29.2,     -1.0,     -1.0,     -1.0]
        elseif Z == 33   wa = [ 11867.0,   1527.0,  1359.1,   1323.6,     204.7,    146.2,    141.2,     41.7,     41.7,     -1.0,     -1.0,     -1.0]
        elseif Z == 34   wa = [ 12658.0,   1652.0,  1474.3,   1433.9,     229.6,    166.5,    160.7,     55.5,     54.6,     -1.0,     -1.0,     -1.0]
        elseif Z == 35   wa = [ 13474.0,   1782.0,  1596.0,   1550.0,     257.0,    189.0,    182.0,     70.0,     69.0,     -1.0,     -1.0,     -1.0]
        elseif Z == 36   wa = [ 14326.0,   1921.0,  1730.9,   1678.4,     292.8,    222.2,    214.4,     95.0,     93.8,     27.5,     14.1,     14.1]
        #
        else  error("Data not available for Z = $Z")
        end
        return( wa )
    end

end # module


text:
using TypeDBClient
using Test

@testset "TypeDBClient" begin
    include("unit/test_common.jl")
    include("unit/test_request_builder.jl")
    include("unit/test_concepts.jl")
end

ext: jl


text:
if pwd()[end - 2:end] != "src"
    cd("src")
end

module Arrhenius
    using YAML
    using NPZ
    using LinearAlgebra
    using SparseArrays

    include("Constants.jl")
    include("DataStructure.jl")
    include("Solution.jl")
    include("Magic.jl")
    include("Thermo.jl")
    include("Kinetics.jl")
    include("Transport.jl")
end

```
Following codes are used during development phase only.
```
using LinearAlgebra
using DifferentialEquations
using Sundials
using Profile
using ForwardDiff
using SparseArrays

R = Arrhenius.R
one_atm = Arrhenius.one_atm

gas = Arrhenius.CreateSolution("../mechanism/gri30.yaml");
const ns = gas.n_species
@show gas.thermo.nasa_high[Arrhenius.species_index(gas, "O2"), :];
Y0 = zeros(ns);
Y0[Arrhenius.species_index(gas, "CH4")] = 0.1;
Y0[Arrhenius.species_index(gas, "O2")] = 0.2;
Y0[Arrhenius.species_index(gas, "N2")] = 0.7;

T0 = 1200.0
P = Arrhenius.one_atm
wdot = Arrhenius.set_states(gas, T0, P, Y0)
@show wdot

T = T0
Y = Y0
mean_MW = 1. / dot(Y, 1 ./ gas.MW)
ρ_mass = P / R / T * mean_MW
X =  Arrhenius.Y2X(gas, Y, mean_MW)
C =  Arrhenius.Y2C(gas, Y, ρ_mass)

# cp_mole, cp_mass =  Arrhenius.get_cp(gas, T, X, mean_MW)
# cv_mole, cv_mass = Arrhenius.get_cv(cp_mole, cp_mass, mean_MW)
#
# h_mole =  Arrhenius.get_H(gas, T, Y, X)
# u_mole = Arrhenius.get_U(h_mole, Y)
#
# S0 =  Arrhenius.get_S(gas, T, P, X)
# wdot = Arrhenius.wdot_func(gas.reaction, T, C, S0, h_mole; get_qdot=false)

# u0 = vcat(Y0, T0);
# function f(u)
#     # Arrhenius.set_states(gas, u[end], P, u[1:ns])[10]
#     T = u[end]
#     Y = @view(u[1:ns])
#     mean_MW = 1. / dot(Y, 1 ./ gas.MW)
#     ρ_mass = P / R / T * mean_MW
#     X =  Arrhenius.Y2X(gas, Y, mean_MW)
#     C =  Arrhenius.Y2C(gas, Y, ρ_mass)
#     cp_mole, cp_mass =  Arrhenius.get_cp(gas, T, X, mean_MW)
#     h_mole =  Arrhenius.get_H(gas, T, Y, X)
#     S0 =  Arrhenius.get_S(gas, T, P, X)
#     wdot = Arrhenius.wdot_func(gas.reaction, T, C, S0, h_mole; get_qdot=false)
#     # wdot = gas.reaction.vk * qdot
#     return wdot[end]
# end
#
# @time f(u0)
# @time grad = ForwardDiff.gradient(f, u0)


text:
using Test

include("isbn-verifier.jl")

@test ISBN <: AbstractString

@testset "valid ISBN numbers" begin
    # ISBN number
    @test isvalid(ISBN, "3-598-21508-8")
    # ISBN number with a check digit of 10
    @test isvalid(ISBN, "3-598-21507-X")
    # ISBN without separating dashes
    @test isvalid(ISBN, "3598215088")
    # ISBN without separating dashes and X as check digit
    @test isvalid(ISBN, "359821507X")
end

@testset "invalid ISBN numbers" begin
    # invalid ISBN check digit
    @test !isvalid(ISBN, "3-598-21508-9")
    # check digit is a character other than X
    @test !isvalid(ISBN, "3-598-21507-A")
    # invalid character in ISBN
    @test !isvalid(ISBN, "3-598-2K507-0")
    # X is only valid as a check isdigit
    @test !isvalid(ISBN, "3-598-2X507-9")
    # ISBN without check digit and dashes
    @test !isvalid(ISBN, "359821507")
    # too long ISBN and no dashes
    @test !isvalid(ISBN, "3598215078X")
    # ISBN without check digit
    @test !isvalid(ISBN, "3-598-21507")
    # too long ISBN
    @test !isvalid(ISBN, "3-598-21507-XX")
    # check digit of X should not be used for 0
    @test !isvalid(ISBN, "3-598-21515-X")
    # empty ISBN
    @test !isvalid(ISBN, "")
    # invalid character in ISBN
    @test !isvalid(ISBN, "3-598-P1581-X")
    # too short ISBN
    @test !isvalid(ISBN, "00")
    # input is 9 characters
    @test !isvalid(ISBN, "134456729")
    # invalid characters are not ignored
    @test !isvalid(ISBN, "3132P34035")
    # input is too long but contains a valid ISBN
    @test !isvalid(ISBN, "98245726788")
end

@testset "constructing valid ISBN numbers" begin
    # ISBN number
    @test string(isbn"3-598-21508-8") == "3598215088"
    # ISBN number with a check digit of 10
    @test string(isbn"3-598-21507-X") == "359821507X"
    # ISBN without separating dashes
    @test string(isbn"3598215088") == "3598215088"
    # ISBN without separating dashes and X as check digit
    @test string(isbn"359821507X") == "359821507X"
end

@testset "constructing invalid ISBN numbers" begin
    # invalid ISBN check digit
    @test_throws DomainError ISBN("3-598-21508-9")
    # check digit is a character other than X
    @test_throws DomainError ISBN("3-598-21507-A")
    # invalid character in ISBN
    @test_throws DomainError ISBN("3-598-2K507-0")
    # X is only valid as a check isdigit
    @test_throws DomainError ISBN("3-598-2X507-9")
    # ISBN without check digit and dashes
    @test_throws DomainError ISBN("359821507")
    # too long ISBN and no dashes
    @test_throws DomainError ISBN("3598215078X")
    # ISBN without check digit
    @test_throws DomainError ISBN("3-598-21507")
    # too long ISBN
    @test_throws DomainError ISBN("3-598-21507-XX")
    # check digit of X should not be used for 0
    @test_throws DomainError ISBN("3-598-21515-X")
    # empty ISBN
    @test_throws DomainError ISBN("")
    # invalid character in ISBN
    @test_throws DomainError ISBN("3-598-P1581-X")
    # too short ISBN
    @test_throws DomainError ISBN("00")
    # input is 9 characters
    @test_throws DomainError ISBN("134456729")
    # invalid characters are not ignored
    @test_throws DomainError ISBN("3132P34035")
    # input is too long but contains a valid ISBN
    @test_throws DomainError ISBN("98245726788")
end


text:
facts("compare_optimizers") do
  context("comparing optimizers on more than one problem") do
    #methods = BlackBoxOptim.MethodNames
    methods = [:adaptive_de_rand_1_bin_radiuslimited, :generating_set_search, :separable_nes, :xnes
      #:random_search,
      #:simultaneous_perturbation_stochastic_approximation
    ]
    problems = BlackBoxOptim.as_fixed_dim_problem_set(BlackBoxOptim.example_problems, [50, 500])
    ranks, fitnesses = BlackBoxOptim.compare_optimizers(problems; max_time = 60.0, methods = methods)
    @fact size(ranks, 2) --> length(problems)
  end
end


text:
"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
-                   MEASUREMENTS / STATE PREPARATION SETTINGS                  -
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""
    fullbases(N::Int; local_basis = "Pauli")

Generate the full set of bases for a choice of local single-qubit basis set.
Predefined option: "Pauli", # bases = 3^N
"""
function fullbases(N::Int; local_basis = "Pauli")
  local_basis == "Pauli" && (local_basis = ["X","Y","Z"])
  if N >15
    print("The $(N)-qubit set of Pauli bases contains $(3^N) bases.\n This may take a while...\n\n")
  end
  !(local_basis isa AbstractArray) && error("Basis not recognized")
  A = Iterators.product(ntuple(i->local_basis, N)...) |> collect
  B = reverse.(reshape(A,length(A),1))
  return  reduce(hcat, getindex.(B,i) for i in 1:N)
end

"""
    fullpreparations(N::Int; local_input_states="Pauli")

Generate the full set of input states to a channel.
Predefined option: "Pauli", 6^N
                   "Tetra", 4^N
"""
function fullpreparations(N::Int; local_input_states="Pauli")
  if N > 5
    print("The $(N)-qubit set of Pauli eigenstates contains $(6^N) bases.\n This may take a while...\n\n")
  end
  local_input_states == "Pauli" && (local_input_states = ["X+", "X-", "Y+", "Y-", "Z+", "Z-"])
  local_input_states == "Tetra" && (local_input_states = ["T1","T2","T3","T4"])
  !(local_input_states isa AbstractArray) && error("States not recognized")
  A = Iterators.product(ntuple(i->local_input_states, N)...) |> collect
  B = reverse.(reshape(A,length(A),1))
  return  reduce(hcat, getindex.(B,i) for i in 1:N)
end

"""
    randombases(N::Int, nbases::Int; local_basis = "Pauli")

Generate `nbases` measurement bases. By default, each
local basis is randomly selected between `["X","Y","Z"]`, with
`"Z"` being the default basis where the quantum state is written.
"""
function randombases(N::Int, nbases::Int; local_basis = "Pauli")
  local_basis == "Pauli" && (local_basis = ["X","Y","Z"])
  return rand(local_basis, nbases, N)
end

"""
    measurementgates(basis::Vector)

Given as input a measurement basis, returns the corresponding
gate data structure. If the basis is `"Z"`, no action is required.
If not, a quantum gate corresponding to the given basis rotation
is added to the list.

Example:
  basis = `["X","Z","Z","Y"]`

  -> `gate_list = [("basisX", 1),("basisY", 4)]`
"""
function measurementgates(basis::Vector)
  gate_list = Tuple[]
  for j in 1:length(basis)
    if basis[j] ≠ "Z"
      push!(gate_list, ("basis$(basis[j])", j, (dag=true,)))
    end
  end
  return gate_list
end

"""
    randompreparations(N::Int, nshots::Int;
                       local_input_state = ["X+","X-","Y+","Y-","Z+","Z-"],
                       ndistinctstates = nothing)

Generate `nshots` input states to a quantum circuit. By default, each
single-qubit state is randomly selected between the 6 eigenstates of
the Pauli matrices, `["X+","X-","Y+","Y-","Z+","Z-"]`.
If `ndistinctstates` is provided, the output consist of `numprep`
different input states, each being repeated `nshots÷ndistinctstates`
times.
"""
function randompreparations(
  N::Int,
  npreps::Int;
  local_input_states= "Pauli",
)

  local_input_states == "Pauli" && (local_input_states = ["X+", "X-", "Y+", "Y-", "Z+", "Z-"])
  local_input_states == "Tetra" && (local_input_states = ["T1","T2","T3","T4"])
  # One shot per basis
  return rand(local_input_states, npreps, N)
end

"""
    readouterror!(measurement::Union{Vector, Matrix}, p1given0, p0given1)

Add readout error to a single projective measurement.

# Arguments:
  - `measurement`: bit string of projective measurement outcome
  - `p1given0`: readout error probability 0 -> 1
  - `p0given1`: readout error probability 1 -> 0
"""
function readouterror!(
  measurement::Union{Vector,Matrix}, p1given0::Float64, p0given1::Float64
)
  for j in 1:size(measurement)[1]
    if measurement[j] == 0
      measurement[j] = StatsBase.sample([0, 1], Weights([1 - p1given0, p1given0]))
    else
      measurement[j] = StatsBase.sample([0, 1], Weights([p0given1, 1 - p0given1]))
    end
  end
  return measurement
end

"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
-                            PROJECTIVE MEASUREMENTS                           -
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""
    getsamples(T::ITensor, nshots::Int)

Generate `nshots` projective measurements for an input quantum state `T`,
which can either be a wavefunction or density matrix (dense).
"""
function getsamples(T::ITensor, nshots::Int;
    readout_errors = (p1given0 = nothing, p0given1 = nothing))

  p1given0 = readout_errors[:p1given0]
  p0given1 = readout_errors[:p0given1]

  # Get the number of qubits
  N = nsites(T)

  # Get a dense array which can be
  # - Vector with dim 2^N for a wavefunction |ψ⟩
  # - Matrix with dim (2^N,2^N) for a density matrix ρ
  A = array(T)

  # Compute the full probability distribution
  # P(σ) = |⟨σ|ψ⟩|² (Tr[ρ|σ⟩⟨σ|]
  probs = (A isa Vector ? abs2.(A) : real(diag(A)))
  @assert sum(probs) ≈ 1

  # Sample the distribution exactly
  index = StatsBase.sample(0:1<<N-1, StatsBase.Weights(probs), nshots)

  # Map integer to binary vectors and massage the structure
  M = hcat(digits.(index,base=2,pad=N)...)'
  measurements = reverse(M,dims=2)
  if !isnothing(p1given0) || !isnothing(p0given1)
    p1given0 = (isnothing(p1given0) ? 0.0 : p1given0)
    p0given1 = (isnothing(p0given1) ? 0.0 : p0given1)
    for n in 1:nshots
      readouterror!(measurements[n,:], p1given0, p0given1)
    end
  end
  return measurements
end

getsamples(T::ITensor; kwargs...) =
  vcat(getsamples!(T,1; kwargs...)...)

"""
    getsamples!(M::Union{MPS,MPO};
                readout_errors = (p1given0 = nothing,
                                  p0given1 = nothing))

Generate a single projective measurement in the MPS/MPO reference basis.
If `readout_errors` is non-trivial, readout errors with given probabilities
are applied to the measurement outcome.
"""
function getsamples!(M::Union{MPS,MPO}; readout_errors=(p1given0=nothing, p0given1=nothing))
  p1given0 = readout_errors[:p1given0]
  p0given1 = readout_errors[:p0given1]
  orthogonalize!(M, 1)
  measurement = sample(M)
  measurement .-= 1
  if !isnothing(p1given0) || !isnothing(p0given1)
    p1given0 = (isnothing(p1given0) ? 0.0 : p1given0)
    p0given1 = (isnothing(p0given1) ? 0.0 : p0given1)
    readouterror!(measurement, p1given0, p0given1)
  end
  return measurement
end

"""
    PastaQ.getsamples!(M::Union{MPS,MPO}, nshots::Int; kwargs...)

Perform `nshots` projective measurements of a wavefunction
`|ψ⟩` or density operator `ρ` in the MPS/MPO reference basis.
Each measurement consists of a binary vector `σ = (σ₁,σ₂,…)`,
drawn from the probabilty distribution:
- `P(σ) = |⟨σ|ψ⟩|²`  :  if `M = ψ is MPS`
- `P(σ) = ⟨σ|ρ|σ⟩`   :  if `M = ρ is MPO`
"""
function getsamples(M0::Union{MPS,MPO}, nshots::Int; kwargs...)
  nthreads = Threads.nthreads()
  data = [Vector{Vector{Int64}}(undef, 0) for _ in 1:nthreads]
  M = orthogonalize!(copy(M0), 1)

  Threads.@threads for j in 1:nshots
    nthread = Threads.threadid()
    sample_ = getsamples!(M; kwargs...)
    push!(data[nthread], sample_)
  end
  return permutedims(hcat(vcat(data...)...))
end

"""
    getsamples(M::Union{MPS,MPO,ITensor}, bases::Matrix; kwargs...)

Generate a dataset of measurements acccording to a set
of input `bases`, performing `nshots` measurements per basis.
For a single measurement, `Û` is the depth-1
local circuit rotating each qubit, the  data-point `σ = (σ₁,σ₂,…)`
is drawn from the probability distribution:
- `P(σ) = |⟨σ|Û|ψ⟩|²`    :  if `M = ψ is MPS`
- `P(σ) = <σ|Û ρ Û†|σ⟩`  :  if `M = ρ is MPO`
"""
function getsamples(M0::Union{MPS,MPO,ITensor}, bases::Matrix{<:String}, nshots::Int; kwargs...)
  N = nsites(M0)
  @assert N == size(bases)[2]
  nthreads = Threads.nthreads()
  data = [Vector{Vector{Pair{String,Int}}}(undef, 0) for _ in 1:nthreads]
  M = copy(M0)
  !(M isa ITensor) && orthogonalize!(M, 1)

  Threads.@threads for n in 1:size(bases, 1)
    nthread = Threads.threadid()
    meas_gates = measurementgates(bases[n, :])
    M_meas = runcircuit(copy(M), meas_gates)
    measurements = getsamples(M_meas, nshots; kwargs...)
    basisdata = [[bases[n,j] => measurements[k,j] for j in 1:N] for k in 1:nshots]
    append!(data[nthread], basisdata)
  end
  return permutedims(hcat(vcat(data...)...))
end

getsamples(M::Union{MPS,MPO,ITensor}, bases::Vector{<:Vector}, nshots::Int; kwargs...) =
  getsamples(M, permutedims(hcat(bases...)), nshots; kwargs...)

getsamples(M::Union{MPS,MPO,ITensor}, bases::Union{Matrix,Vector{<:Vector}}; kwargs...) =
  getsamples(M, bases, 1; kwargs...)

"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
                            PROCESS TOMOGRAPHY DATA                                -
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""
    projectchannel(U0::MPO,prep::Array)
    projectchannel(U::ITensor, prep::Array)

Project the unitary circuit (MPO) into a state `prep`
made out of single-qubit Pauli eigenstates (e.g. `|ϕ⟩ =|+⟩⊗|0⟩⊗|r⟩⊗…).
The resulting MPS describes the quantum state obtained by applying
the quantum circuit to `|ϕ⟩`. Same for a dense ITensors.
"""
function _projectchannel(U0::Union{MPO,ITensor}, prep::Array)
  U = copy(U0)
  s = U isa MPO ? firstsiteinds(U) : inds(U; plev = 0)
  for j in 1:length(s)
    if U isa MPO
      U[j] = U[j] * state(prep[j], s[j])
    else
      U = U * state(prep[j], s[j])
    end
  end
  U isa ITensor && return noprime!(U)
  return convert(MPS,noprime!(U))
end

function _projectchannel(Λ0::Choi, prep::Array)
  Λ = copy(Λ0.X)
  s = Λ isa MPO ? firstsiteinds(Λ; tags="Input") : inds(Λ; tags = "Input", plev = 0)
  for j in 1:length(s)
    # No conjugate on the gate (transpose input!)
    if Λ isa MPO
      Λ[j] = Λ[j] * dag(state(prep[j], s[j]))
      Λ[j] = Λ[j] * prime(state(prep[j], s[j]))
    else
      Λ = Λ * dag(state(prep[j], s[j]))
      Λ = Λ * prime(state(prep[j], s[j]))
    end
  end
  return Λ
end

projectchannel(M::Union{MPO,ITensor}, prep::AbstractArray) =
  ischoi(M) ? _projectchannel(Choi(M), prep) : _projectchannel(M, prep)

function getsamples(
  M0::Union{LPDO,MPO,ITensor},
  preps_and_bases::Matrix{<:Pair},
  nshots::Int,
  kwargs...
)
  N = nsites(M0)
  @assert N == size(preps_and_bases, 2)

  nthreads = Threads.nthreads()
  data = [Matrix{Pair{String,Pair{String,Int}}}(undef, 0, N) for _ in 1:nthreads]
  M = copy(M0)

  Threads.@threads for k in 1:size(preps_and_bases, 1)
    nthread = Threads.threadid()
    # input state to the channel
    prep  = first.(preps_and_bases[k,:])
    # measurement basis
    basis = last.(preps_and_bases[k,:])

    # project the channel into the input state
    Mproj = projectchannel(M, prep)

    # perform measurement basis rotation
    meas_gates = measurementgates(basis)
    M_meas = runcircuit(Mproj, meas_gates)

    # perform projective measurement
    measurements = getsamples(M_meas, nshots; kwargs...)
    basisdata = [[prep[j] => (basis[j] => measurements[k,j]) for j in 1:N] for k in 1:nshots]
    data[nthread] = vcat(data[nthread], permutedims(hcat(basisdata...)))
  end
  return vcat(data...)
end

function getsamples(
  M0::Union{LPDO,MPO,ITensor},
  preps::Matrix,
  bases::Matrix,
  nshots::Int;
  kwargs...
)

  N = nsites(M0)
  npreps = size(preps, 1)
  nbases = size(bases, 1)
  ntotal = npreps * nbases

  preps_and_bases = Matrix{Pair{String,String}}(undef, 0, N)
  for p in 1:npreps
    x = [preps[p,:] .=> bases[b,:] for b in 1:nbases]
    preps_and_bases = vcat(preps_and_bases, permutedims(hcat(x...)))
  end
  return getsamples(M0, preps_and_bases, nshots; kwargs...)
end

getsamples(M::Union{LPDO,MPO, ITensor}, preps::Vector{<:Vector}, bases::Vector{<:Vector}, nshots::Int; kwargs...) =
  getsamples(M, permutedims(hcat(preps...)), permutedims(hcat(bases...)), nshots; kwargs...)

getsamples(M::Union{LPDO,MPO, ITensor}, preps::Union{Matrix,Vector{<:Vector}}, bases::Union{Matrix,Vector{<:Vector}}; kwargs...) =
  getsamples(M, preps, bases, 1; kwargs...)


text:
using DSGE, ModelConstructors, HDF5, Random, JLD2, FileIO, SMC, Test

path = dirname(@__FILE__)
writing_output = false
if VERSION < v"1.5"
    ver = "111"
elseif VERSION < v"1.6"
    ver = "150"
else
    ver = "160"
end

m = AnSchorfheide()

save = normpath(joinpath(dirname(@__FILE__), "save"))
m <= Setting(:saveroot, save)

data = h5read(joinpath(path, "reference/smc.h5"), "data")

m <= Setting(:n_particles, 400)
m <= Setting(:n_Φ, 100)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, true)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:adaptive_tempering_target_smc, false)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:smc_iteration, 0)
m <= Setting(:use_chand_recursion, true)

@everywhere Random.seed!(42)

println("Estimating AnSchorfheide Model... (approx. 2 minutes)")
DSGE.smc2(m, data, verbose = :none, run_csminwel = false) # Uncomment once new DSGE version comes out: run_csminwel = false, verbose = :none) # us.txt gives equiv to periods 95:174 in our current dataset
println("Estimation done!")

test_file = load(rawpath(m, "estimate", "smc_cloud.jld2"))
test_cloud  = test_file["cloud"]
test_w      = test_file["w"]
test_W      = test_file["W"]

if writing_output
    jldopen(joinpath(path, "reference/smc_cloud_fix=true_version=" * ver * ".jld2"),
            true, true, true, IOStream) do file
        write(file, "cloud", test_cloud)
        write(file, "w", test_w)
        write(file, "W", test_W)
    end
end

saved_file  = load(joinpath(path, "reference/smc_cloud_fix=true_version=" * ver * ".jld2"))
saved_cloud = saved_file["cloud"]
saved_w     = saved_file["w"]
saved_W     = saved_file["W"]

####################################################################
cloud_fields = fieldnames(typeof(test_cloud))
@testset "ParticleCloud Fields: AnSchorf" begin
    @test @test_matrix_approx_eq SMC.get_vals(test_cloud) SMC.get_vals(saved_cloud)
    @test @test_matrix_approx_eq SMC.get_loglh(test_cloud) SMC.get_loglh(saved_cloud)
    @test length(test_cloud.particles) == length(saved_cloud.particles)
    @test test_cloud.tempering_schedule == saved_cloud.tempering_schedule
    @test test_cloud.ESS ≈ saved_cloud.ESS
    @test test_cloud.stage_index == saved_cloud.stage_index
    @test test_cloud.n_Φ == saved_cloud.n_Φ
    @test test_cloud.resamples == saved_cloud.resamples
    @test test_cloud.c == saved_cloud.c
    @test test_cloud.accept == saved_cloud.accept
end

test_particle  = test_cloud.particles[1,:]
saved_particle = saved_cloud.particles[1,:]
N = length(test_particle)
@testset "Individual Particle Fields Post-SMC: AnSchorf" begin
    @test test_particle[1:SMC.ind_para_end(N)] ≈ saved_particle[1:SMC.ind_para_end(N)]
    @test test_particle[SMC.ind_loglh(N)]      ≈ saved_particle[SMC.ind_loglh(N)]
    @test test_particle[SMC.ind_logprior(N)]   ≈ saved_particle[SMC.ind_logprior(N)]
    @test test_particle[SMC.ind_old_loglh(N)] == saved_particle[SMC.ind_old_loglh(N)]
    @test test_particle[SMC.ind_accept(N)]    == saved_particle[SMC.ind_accept(N)]
    @test test_particle[SMC.ind_weight(N)]     ≈ saved_particle[SMC.ind_weight(N)]
end

@testset "Weight Matrices: AnSchorf" begin
    @test @test_matrix_approx_eq test_w saved_w
    @test @test_matrix_approx_eq test_W saved_W
end

####################################################################
# Bridging Test
####################################################################
#=
# This bridging test doesn't work, but the one in SMC.jl works
m = AnSchorfheide()

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

data = h5read(joinpath(path, "reference/smc.h5"), "data")

m <= Setting(:n_particles, 400)
m <= Setting(:n_Φ, 100)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, true)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:adaptive_tempering_target_smc, false)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:smc_iteration, 0)
m <= Setting(:use_chand_recursion, true)

@everywhere Random.seed!(42)

# Estimate with 1st half of sample
m_old = deepcopy(m)
m_old <= Setting(:n_particles, 1000, true, "npart", "")
m_old <= Setting(:data_vintage, "000000")
#DSGE.smc2(m_old, data[:,1:Int(floor(end/2))], run_csminwel = false, verbose = :low)

m_new = deepcopy(m)

# Estimate with 2nd half of sample
m_new <= Setting(:data_vintage, "200218")
m_new <= Setting(:tempered_update_prior_weight, 0.5)
#m_new <= Setting(:tempered_update, true)
old_vint = "000000"
m_new <= Setting(:previous_data_vintage, old_vint)
loadpath = rawpath(m_old, "estimate", "smc_cloud.jld2")
loadpath = replace(loadpath, r"vint=[0-9]{6}" => "vint=" * old_vint)
old_cloud = ParticleCloud(load(loadpath, "cloud"), map(x -> x.key, m.parameters))
m_new <= Setting(:n_particles, 2000, true, "npart", "")
DSGE.smc2(m_new, data,
          old_data = data[:,1:Int(floor(end/2))], old_cloud = old_cloud,
          run_csminwel = false)

loadpath = rawpath(m_new, "estimate", "smc_cloud.jld2")
new_cloud = ParticleCloud(load(loadpath, "cloud"), map(x -> x.key, m.parameters))

# Clean output files up
rm(rawpath(m_new, "estimate", "smc_cloud.jld2"))
rm(rawpath(m_new, "estimate", "smcsave.h5"))
=#


text:
refVals = []
refPrecs = []

#! format: off
# SC ========== Test number 1 reference values and precision match template. =======
# SC ========== /home/jmc/cliMa/cliMa_update/test/Ocean/SplitExplicit/simple_box_ivd.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [  "oce Q_3D",   "u[1]", -1.50907071392284014566e-01,  1.73781231809542469069e-01,  4.68096331745367338212e-03,  1.95920021000241043052e-02 ],
 [  "oce Q_3D",   "u[2]", -1.41693937396426050679e-01,  1.46243917821118285527e-01, -2.51530766341314954843e-03,  1.75572837799599749953e-02 ],
 [  "oce Q_3D",       :η, -7.42455782775250150429e-01,  3.69730515308098417471e-01, -2.12053967392545028277e-04,  2.72181302813857772804e-01 ],
 [  "oce Q_3D",       :θ,  5.28450920905325794569e-03,  9.92875001597955630928e+00,  2.49640763974860568908e+00,  2.18014403700940651021e+00 ],
 [   "oce aux",       :w, -1.93641508217356502321e-04,  1.64389722003552838422e-04,  4.30819768270874037293e-07,  1.44751504460569400376e-05 ],
 [   "oce aux",    :pkin, -8.85137385724656411412e+00,  0.00000000000000000000e+00, -3.26232394846976436753e+00,  2.50203581939964081471e+00 ],
 [   "oce aux",     :wz0, -2.13126994266848044583e-05,  3.03761966887832467763e-05, -1.29724204981140911883e-10,  7.79125622163024900916e-06 ],
 [   "oce aux", "u_d[1]", -1.47643597045305691173e-01,  1.16172928281001314188e-01, -3.40056303283133183873e-05,  1.18325769360269061892e-02 ],
 [   "oce aux", "u_d[2]", -1.41299443151196246760e-01,  1.40392329116052066995e-01, -8.68704744435173440418e-06,  1.16473030652601838852e-02 ],
 [   "oce aux", "ΔGu[1]", -1.82438482660877217949e-06,  1.76131450220044681262e-06, -3.58493343107863405131e-08,  2.63467783463959988175e-07 ],
 [   "oce aux", "ΔGu[2]", -1.21868069222541804602e-06,  2.22130916543395098495e-06,  1.31060034108175441821e-06,  6.57246261117217741340e-07 ],
 [   "oce aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15573163901915703900e+06 ],
 [ "baro Q_2D",   "U[1]", -1.28638134009455562534e+01,  6.10305774663233222554e+01,  4.70737872487603414839e+00,  1.49617848161145641228e+01 ],
 [ "baro Q_2D",   "U[2]", -2.82511172627316504702e+01,  6.47636399564029829889e+01, -2.47937561262140748752e+00,  1.34334006549456468349e+01 ],
 [ "baro Q_2D",       :η, -7.42661610672267546995e-01,  3.69774652167159434413e-01, -2.12086144227697549498e-04,  2.72331664923152527713e-01 ],
 [  "baro aux",  "Gᵁ[1]", -1.76131450220044674486e-03,  1.82438482660877217441e-03,  3.58493343107863313016e-05,  2.63480826100707878241e-04 ],
 [  "baro aux",  "Gᵁ[2]", -2.22130916543395113064e-03,  1.21868069222541802570e-03, -1.31060034108175450969e-03,  6.57278797255507452567e-04 ],
 [  "baro aux", "U_c[1]", -1.28418434709131652482e+01,  6.09600328717265398382e+01,  4.71373811044994983632e+00,  1.49477163957520815529e+01 ],
 [  "baro aux", "U_c[2]", -2.82747850743473101431e+01,  6.47179561285744426868e+01, -2.50652647306075504474e+00,  1.34314270881023727355e+01 ],
 [  "baro aux",     :η_c, -7.42455782775250150429e-01,  3.69730515308098417471e-01, -2.12053967392540691468e-04,  2.72194776802272997429e-01 ],
 [  "baro aux", "U_s[1]", -1.28420859836829617251e+01,  6.09610745104530664662e+01,  4.71375870095391835690e+00,  1.49479160464055684798e+01 ],
 [  "baro aux", "U_s[2]", -2.82756417530404320360e+01,  6.47178077635098105702e+01, -2.50659230751473804943e+00,  1.34314567895758454341e+01 ],
 [  "baro aux",     :η_s, -7.42463591584554216674e-01,  3.69734792872739193026e-01, -2.12048722047518373819e-04,  2.72195582254090906460e-01 ],
 [  "baro aux",  "Δu[1]", -4.52576313913315424631e-04,  4.01397574102695402134e-04, -1.17811088890548056881e-05,  1.01680622362620656394e-04 ],
 [  "baro aux",  "Δu[2]", -2.43950808154872009637e-04,  3.86395073848479621119e-04,  5.11424556860824746545e-05,  7.52605326838232742588e-05 ],
 [  "baro aux",  :η_diag, -7.39883203465319660985e-01,  3.68215473126410952620e-01, -2.11645242010587240206e-04,  2.71980481766462722781e-01 ],
 [  "baro aux",      :Δη, -3.57633535239532118766e-03,  4.28419930737145016053e-03, -4.08725381931473987174e-07,  1.12670824245143363224e-03 ],
 [  "baro aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15578885204060329124e+06 ],
]
parr = [
 [  "oce Q_3D",   "u[1]",    12,    12,    12,    12 ],
 [  "oce Q_3D",   "u[2]",    12,    12,    12,    12 ],
 [  "oce Q_3D",       :η,    12,    12,     8,    12 ],
 [  "oce Q_3D",       :θ,    12,    12,    12,    12 ],
 [   "oce aux",       :w,    12,    12,     8,    12 ],
 [   "oce aux",    :pkin,    12,    12,    12,    12 ],
 [   "oce aux",     :wz0,    12,    12,     8,    12 ],
 [   "oce aux", "u_d[1]",    12,    12,    12,    12 ],
 [   "oce aux", "u_d[2]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[1]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[2]",    12,    12,    12,    12 ],
 [   "oce aux",       :y,    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[1]",    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[2]",    12,    12,    12,    12 ],
 [ "baro Q_2D",       :η,    12,    12,     8,    12 ],
 [  "baro aux",  "Gᵁ[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Gᵁ[2]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[2]",    12,    12,    12,    12 ],
 [  "baro aux",     :η_c,    12,    12,     8,    12 ],
 [  "baro aux", "U_s[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_s[2]",    12,    12,    12,    12 ],
 [  "baro aux",     :η_s,    12,    12,     8,    12 ],
 [  "baro aux",  "Δu[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Δu[2]",    12,    12,    12,    12 ],
 [  "baro aux",  :η_diag,    12,    12,     8,    12 ],
 [  "baro aux",      :Δη,    12,    12,     8,    12 ],
 [  "baro aux",       :y,    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

#! format: on

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
println("Launching")
using MocosSimLauncher

launch()

max_stars_repo_name: Novirmor/MocosSimLauncher

max_issues_repo_name: Novirmor/MocosSimLauncher

max_forks_repo_name: Novirmor/MocosSimLauncher


text:
# Containers for SCIP structures
import Base: pointer

type SCIP_t
    array_ptr_scip::Array{Ptr{SCIP}}
end
type SCIP_BDCHGIDX_t
    array_ptr_scip_bdchgidx::Array{Ptr{SCIP_BDCHGIDX}}
end
type SCIP_BDCHGINFO_t
    array_ptr_scip_bdchginfo::Array{Ptr{SCIP_BDCHGINFO}}
end
type SCIP_BILINTERM_t
    array_ptr_scip_bilinterm::Array{Ptr{SCIP_BILINTERM}}
end
type SCIP_Bool_t
    array_ptr_scip_bool::Array{Ptr{SCIP_Bool}}
end
type SCIP_BOOLARRAY_t
    array_ptr_scip_boolarray::Array{Ptr{SCIP_BOOLARRAY}}
end
type SCIP_BOUNDCHG_t
    array_ptr_scip_boundchg::Array{Ptr{SCIP_BOUNDCHG}}
end
type SCIP_BOUNDTYPE_t
    array_ptr_scip_boundtype::Array{Ptr{SCIP_BOUNDTYPE}}
end
type SCIP_BRANCHRULE_t
    array_ptr_scip_branchrule::Array{Ptr{SCIP_BRANCHRULE}}
end
type SCIP_BRANCHRULEDATA_t
    array_ptr_scip_branchruledata::Array{Ptr{SCIP_BRANCHRULEDATA}}
end
type SCIP_BT_t
    array_ptr_scip_bt::Array{Ptr{SCIP_BT}}
end
type SCIP_BTNODE_t
    array_ptr_scip_btnode::Array{Ptr{SCIP_BTNODE}}
end
type SCIP_CLIQUE_t
    array_ptr_scip_clique::Array{Ptr{SCIP_CLIQUE}}
end
type SCIP_CLOCK_t
    array_ptr_scip_clock::Array{Ptr{SCIP_CLOCK}}
end
type SCIP_COL_t
    array_ptr_scip_col::Array{Ptr{SCIP_COL}}
end
type SCIP_CONFLICTHDLR_t
    array_ptr_scip_conflicthdlr::Array{Ptr{SCIP_CONFLICTHDLR}}
end
type SCIP_CONFLICTHDLRDATA_t
    array_ptr_scip_conflicthdlrdata::Array{Ptr{SCIP_CONFLICTHDLRDATA}}
end
type SCIP_CONS_t
    array_ptr_scip_cons::Array{Ptr{SCIP_CONS}}
end
type SCIP_CONSDATA_t
    array_ptr_scip_consdata::Array{Ptr{SCIP_CONSDATA}}
end
type SCIP_CONSHDLR_t
    array_ptr_scip_conshdlr::Array{Ptr{SCIP_CONSHDLR}}
end
type SCIP_CONSHDLRDATA_t
    array_ptr_scip_conshdlrdata::Array{Ptr{SCIP_CONSHDLRDATA}}
end
type SCIP_CUT_t
    array_ptr_scip_cut::Array{Ptr{SCIP_CUT}}
end
type SCIP_CUTPOOL_t
    array_ptr_scip_cutpool::Array{Ptr{SCIP_CUTPOOL}}
end
type SCIP_DIALOG_t
    array_ptr_scip_dialog::Array{Ptr{SCIP_DIALOG}}
end
type SCIP_DIALOGDATA_t
    array_ptr_scip_dialogdata::Array{Ptr{SCIP_DIALOGDATA}}
end
type SCIP_DIALOGHDLR_t
    array_ptr_scip_dialoghdlr::Array{Ptr{SCIP_DIALOGHDLR}}
end
type SCIP_DIGRAPH_t
    array_ptr_scip_digraph::Array{Ptr{SCIP_DIGRAPH}}
end
type SCIP_DISP_t
    array_ptr_scip_disp::Array{Ptr{SCIP_DISP}}
end
type SCIP_DISPDATA_t
    array_ptr_scip_dispdata::Array{Ptr{SCIP_DISPDATA}}
end
type SCIP_DOMCHG_t
    array_ptr_scip_domchg::Array{Ptr{SCIP_DOMCHG}}
end
type SCIP_EVENT_t
    array_ptr_scip_event::Array{Ptr{SCIP_EVENT}}
end
type SCIP_EVENTDATA_t
    array_ptr_scip_eventdata::Array{Ptr{SCIP_EVENTDATA}}
end
type SCIP_EVENTHDLR_t
    array_ptr_scip_eventhdlr::Array{Ptr{SCIP_EVENTHDLR}}
end
type SCIP_EVENTHDLRDATA_t
    array_ptr_scip_eventhdlrdata::Array{Ptr{SCIP_EVENTHDLRDATA}}
end
type SCIP_EXPR_t
    array_ptr_scip_expr::Array{Ptr{SCIP_EXPR}}
end
type SCIP_EXPRCURV_t
    array_ptr_scip_exprcurv::Array{Ptr{SCIP_EXPRCURV}}
end
type SCIP_EXPRDATA_MONOMIAL_t
    array_ptr_scip_exprdata_monomial::Array{Ptr{SCIP_EXPRDATA_MONOMIAL}}
end
type SCIP_EXPRGRAPH_t
    array_ptr_scip_exprgraph::Array{Ptr{SCIP_EXPRGRAPH}}
end
type SCIP_EXPRGRAPHNODE_t
    array_ptr_scip_exprgraphnode::Array{Ptr{SCIP_EXPRGRAPHNODE}}
end
type SCIP_EXPRINTDATA_t
    array_ptr_scip_exprintdata::Array{Ptr{SCIP_EXPRINTDATA}}
end
type SCIP_EXPRTREE_t
    array_ptr_scip_exprtree::Array{Ptr{SCIP_EXPRTREE}}
end
type SCIP_HASHMAP_t
    array_ptr_scip_hashmap::Array{Ptr{SCIP_HASHMAP}}
end
type SCIP_HASHMAPLIST_t
    array_ptr_scip_hashmaplist::Array{Ptr{SCIP_HASHMAPLIST}}
end
type SCIP_HASHTABLE_t
    array_ptr_scip_hashtable::Array{Ptr{SCIP_HASHTABLE}}
end
type SCIP_HEUR_t
    array_ptr_scip_heur::Array{Ptr{SCIP_HEUR}}
end
type SCIP_HEURDATA_t
    array_ptr_scip_heurdata::Array{Ptr{SCIP_HEURDATA}}
end
type SCIP_HISTORY_t
    array_ptr_scip_history::Array{Ptr{SCIP_HISTORY}}
end
type SCIP_HOLELIST_t
    array_ptr_scip_holelist::Array{Ptr{SCIP_HOLELIST}}
end
type SCIP_INTARRAY_t
    array_ptr_scip_intarray::Array{Ptr{SCIP_INTARRAY}}
end
type SCIP_INTERVAL_t
    array_ptr_scip_interval::Array{Ptr{SCIP_INTERVAL}}
end
type SCIP_LPI_t
    array_ptr_scip_lpi::Array{Ptr{SCIP_LPI}}
end
type SCIP_MESSAGEHDLR_t
    array_ptr_scip_messagehdlr::Array{Ptr{SCIP_MESSAGEHDLR}}
end
type SCIP_MESSAGEHDLRDATA_t
    array_ptr_scip_messagehdlrdata::Array{Ptr{SCIP_MESSAGEHDLRDATA}}
end
type SCIP_NLPI_t
    array_ptr_scip_nlpi::Array{Ptr{SCIP_NLPI}}
end
type SCIP_NLPIPROBLEM_t
    array_ptr_scip_nlpiproblem::Array{Ptr{SCIP_NLPIPROBLEM}}
end
type SCIP_NLPSTATISTICS_t
    array_ptr_scip_nlpstatistics::Array{Ptr{SCIP_NLPSTATISTICS}}
end
type SCIP_NLROW_t
    array_ptr_scip_nlrow::Array{Ptr{SCIP_NLROW}}
end
type SCIP_NODE_t
    array_ptr_scip_node::Array{Ptr{SCIP_NODE}}
end
type SCIP_NODESEL_t
    array_ptr_scip_nodesel::Array{Ptr{SCIP_NODESEL}}
end
type SCIP_NODESELDATA_t
    array_ptr_scip_nodeseldata::Array{Ptr{SCIP_NODESELDATA}}
end
type SCIP_PARAM_t
    array_ptr_scip_param::Array{Ptr{SCIP_PARAM}}
end
type SCIP_PQUEUE_t
    array_ptr_scip_pqueue::Array{Ptr{SCIP_PQUEUE}}
end
type SCIP_PRESOL_t
    array_ptr_scip_presol::Array{Ptr{SCIP_PRESOL}}
end
type SCIP_PRESOLDATA_t
    array_ptr_scip_presoldata::Array{Ptr{SCIP_PRESOLDATA}}
end
type SCIP_PRICER_t
    array_ptr_scip_pricer::Array{Ptr{SCIP_PRICER}}
end
type SCIP_PRICERDATA_t
    array_ptr_scip_pricerdata::Array{Ptr{SCIP_PRICERDATA}}
end
type SCIP_PROBDATA_t
    array_ptr_scip_probdata::Array{Ptr{SCIP_PROBDATA}}
end
type SCIP_PROFILE_t
    array_ptr_scip_profile::Array{Ptr{SCIP_PROFILE}}
end
type SCIP_PROP_t
    array_ptr_scip_prop::Array{Ptr{SCIP_PROP}}
end
type SCIP_PROPDATA_t
    array_ptr_scip_propdata::Array{Ptr{SCIP_PROPDATA}}
end
type SCIP_PTRARRAY_t
    array_ptr_scip_ptrarray::Array{Ptr{SCIP_PTRARRAY}}
end
type SCIP_QUADELEM_t
    array_ptr_scip_quadelem::Array{Ptr{SCIP_QUADELEM}}
end
type SCIP_QUADVARTERM_t
    array_ptr_scip_quadvarterm::Array{Ptr{SCIP_QUADVARTERM}}
end
type SCIP_QUEUE_t
    array_ptr_scip_queue::Array{Ptr{SCIP_QUEUE}}
end
type SCIP_READER_t
    array_ptr_scip_reader::Array{Ptr{SCIP_READER}}
end
type SCIP_READERDATA_t
    array_ptr_scip_readerdata::Array{Ptr{SCIP_READERDATA}}
end
type SCIP_Real_t
    array_ptr_scip_real::Array{Ptr{SCIP_Real}}
end
type SCIP_REALARRAY_t
    array_ptr_scip_realarray::Array{Ptr{SCIP_REALARRAY}}
end
type SCIP_RELAX_t
    array_ptr_scip_relax::Array{Ptr{SCIP_RELAX}}
end
type SCIP_RELAXDATA_t
    array_ptr_scip_relaxdata::Array{Ptr{SCIP_RELAXDATA}}
end
type SCIP_RESOURCEACTIVITY_t
    array_ptr_scip_resourceactivity::Array{Ptr{SCIP_RESOURCEACTIVITY}}
end
type SCIP_RESULT_t
    array_ptr_scip_result::Array{Ptr{SCIP_RESULT}}
end
type SCIP_ROW_t
    array_ptr_scip_row::Array{Ptr{SCIP_ROW}}
end
type SCIP_SEPA_t
    array_ptr_scip_sepa::Array{Ptr{SCIP_SEPA}}
end
type SCIP_SEPADATA_t
    array_ptr_scip_sepadata::Array{Ptr{SCIP_SEPADATA}}
end
type SCIP_SOL_t
    array_ptr_scip_sol::Array{Ptr{SCIP_SOL}}
end
type SCIP_SPARSESOL_t
    array_ptr_scip_sparsesol::Array{Ptr{SCIP_SPARSESOL}}
end
type SCIP_VALUEHISTORY_t
    array_ptr_scip_valuehistory::Array{Ptr{SCIP_VALUEHISTORY}}
end
type SCIP_VAR_t
    array_ptr_scip_var::Array{Ptr{SCIP_VAR}}
end
type SCIP_VARDATA_t
    array_ptr_scip_vardata::Array{Ptr{SCIP_VARDATA}}
end

# Pointer and array access functions
array(scip::SCIP_t) = scip.array_ptr_scip
pointer(scip::SCIP_t) = array(scip)[1]
array(scip_bdchgidx::SCIP_BDCHGIDX_t) = scip_bdchgidx.array_ptr_scip_bdchgidx
pointer(scip_bdchgidx::SCIP_BDCHGIDX_t) = array(scip_bdchgidx)[1]
array(scip_bdchginfo::SCIP_BDCHGINFO_t) = scip_bdchginfo.array_ptr_scip_bdchginfo
pointer(scip_bdchginfo::SCIP_BDCHGINFO_t) = array(scip_bdchginfo)[1]
array(scip_bilinterm::SCIP_BILINTERM_t) = scip_bilinterm.array_ptr_scip_bilinterm
pointer(scip_bilinterm::SCIP_BILINTERM_t) = array(scip_bilinterm)[1]
array(scip_bool::SCIP_Bool_t) = scip_bool.array_ptr_scip_bool
pointer(scip_bool::SCIP_Bool_t) = array(scip_bool)[1]
array(scip_boolarray::SCIP_BOOLARRAY_t) = scip_boolarray.array_ptr_scip_boolarray
pointer(scip_boolarray::SCIP_BOOLARRAY_t) = array(scip_boolarray)[1]
array(scip_boundchg::SCIP_BOUNDCHG_t) = scip_boundchg.array_ptr_scip_boundchg
pointer(scip_boundchg::SCIP_BOUNDCHG_t) = array(scip_boundchg)[1]
array(scip_boundtype::SCIP_BOUNDTYPE_t) = scip_boundtype.array_ptr_scip_boundtype
pointer(scip_boundtype::SCIP_BOUNDTYPE_t) = array(scip_boundtype)[1]
array(scip_branchrule::SCIP_BRANCHRULE_t) = scip_branchrule.array_ptr_scip_branchrule
pointer(scip_branchrule::SCIP_BRANCHRULE_t) = array(scip_branchrule)[1]
array(scip_branchruledata::SCIP_BRANCHRULEDATA_t) = scip_branchruledata.array_ptr_scip_branchruledata
pointer(scip_branchruledata::SCIP_BRANCHRULEDATA_t) = array(scip_branchruledata)[1]
array(scip_bt::SCIP_BT_t) = scip_bt.array_ptr_scip_bt
pointer(scip_bt::SCIP_BT_t) = array(scip_bt)[1]
array(scip_btnode::SCIP_BTNODE_t) = scip_btnode.array_ptr_scip_btnode
pointer(scip_btnode::SCIP_BTNODE_t) = array(scip_btnode)[1]
array(scip_clique::SCIP_CLIQUE_t) = scip_clique.array_ptr_scip_clique
pointer(scip_clique::SCIP_CLIQUE_t) = array(scip_clique)[1]
array(scip_clock::SCIP_CLOCK_t) = scip_clock.array_ptr_scip_clock
pointer(scip_clock::SCIP_CLOCK_t) = array(scip_clock)[1]
array(scip_col::SCIP_COL_t) = scip_col.array_ptr_scip_col
pointer(scip_col::SCIP_COL_t) = array(scip_col)[1]
array(scip_conflicthdlr::SCIP_CONFLICTHDLR_t) = scip_conflicthdlr.array_ptr_scip_conflicthdlr
pointer(scip_conflicthdlr::SCIP_CONFLICTHDLR_t) = array(scip_conflicthdlr)[1]
array(scip_conflicthdlrdata::SCIP_CONFLICTHDLRDATA_t) = scip_conflicthdlrdata.array_ptr_scip_conflicthdlrdata
pointer(scip_conflicthdlrdata::SCIP_CONFLICTHDLRDATA_t) = array(scip_conflicthdlrdata)[1]
array(scip_cons::SCIP_CONS_t) = scip_cons.array_ptr_scip_cons
pointer(scip_cons::SCIP_CONS_t) = array(scip_cons)[1]
array(scip_consdata::SCIP_CONSDATA_t) = scip_consdata.array_ptr_scip_consdata
pointer(scip_consdata::SCIP_CONSDATA_t) = array(scip_consdata)[1]
array(scip_conshdlr::SCIP_CONSHDLR_t) = scip_conshdlr.array_ptr_scip_conshdlr
pointer(scip_conshdlr::SCIP_CONSHDLR_t) = array(scip_conshdlr)[1]
array(scip_conshdlrdata::SCIP_CONSHDLRDATA_t) = scip_conshdlrdata.array_ptr_scip_conshdlrdata
pointer(scip_conshdlrdata::SCIP_CONSHDLRDATA_t) = array(scip_conshdlrdata)[1]
array(scip_cut::SCIP_CUT_t) = scip_cut.array_ptr_scip_cut
pointer(scip_cut::SCIP_CUT_t) = array(scip_cut)[1]
array(scip_cutpool::SCIP_CUTPOOL_t) = scip_cutpool.array_ptr_scip_cutpool
pointer(scip_cutpool::SCIP_CUTPOOL_t) = array(scip_cutpool)[1]
array(scip_dialog::SCIP_DIALOG_t) = scip_dialog.array_ptr_scip_dialog
pointer(scip_dialog::SCIP_DIALOG_t) = array(scip_dialog)[1]
array(scip_dialogdata::SCIP_DIALOGDATA_t) = scip_dialogdata.array_ptr_scip_dialogdata
pointer(scip_dialogdata::SCIP_DIALOGDATA_t) = array(scip_dialogdata)[1]
array(scip_dialoghdlr::SCIP_DIALOGHDLR_t) = scip_dialoghdlr.array_ptr_scip_dialoghdlr
pointer(scip_dialoghdlr::SCIP_DIALOGHDLR_t) = array(scip_dialoghdlr)[1]
array(scip_digraph::SCIP_DIGRAPH_t) = scip_digraph.array_ptr_scip_digraph
pointer(scip_digraph::SCIP_DIGRAPH_t) = array(scip_digraph)[1]
array(scip_disp::SCIP_DISP_t) = scip_disp.array_ptr_scip_disp
pointer(scip_disp::SCIP_DISP_t) = array(scip_disp)[1]
array(scip_dispdata::SCIP_DISPDATA_t) = scip_dispdata.array_ptr_scip_dispdata
pointer(scip_dispdata::SCIP_DISPDATA_t) = array(scip_dispdata)[1]
array(scip_domchg::SCIP_DOMCHG_t) = scip_domchg.array_ptr_scip_domchg
pointer(scip_domchg::SCIP_DOMCHG_t) = array(scip_domchg)[1]
array(scip_event::SCIP_EVENT_t) = scip_event.array_ptr_scip_event
pointer(scip_event::SCIP_EVENT_t) = array(scip_event)[1]
array(scip_eventdata::SCIP_EVENTDATA_t) = scip_eventdata.array_ptr_scip_eventdata
pointer(scip_eventdata::SCIP_EVENTDATA_t) = array(scip_eventdata)[1]
array(scip_eventhdlr::SCIP_EVENTHDLR_t) = scip_eventhdlr.array_ptr_scip_eventhdlr
pointer(scip_eventhdlr::SCIP_EVENTHDLR_t) = array(scip_eventhdlr)[1]
array(scip_eventhdlrdata::SCIP_EVENTHDLRDATA_t) = scip_eventhdlrdata.array_ptr_scip_eventhdlrdata
pointer(scip_eventhdlrdata::SCIP_EVENTHDLRDATA_t) = array(scip_eventhdlrdata)[1]
array(scip_expr::SCIP_EXPR_t) = scip_expr.array_ptr_scip_expr
pointer(scip_expr::SCIP_EXPR_t) = array(scip_expr)[1]
array(scip_exprcurv::SCIP_EXPRCURV_t) = scip_exprcurv.array_ptr_scip_exprcurv
pointer(scip_exprcurv::SCIP_EXPRCURV_t) = array(scip_exprcurv)[1]
array(scip_exprdata_monomial::SCIP_EXPRDATA_MONOMIAL_t) = scip_exprdata_monomial.array_ptr_scip_exprdata_monomial
pointer(scip_exprdata_monomial::SCIP_EXPRDATA_MONOMIAL_t) = array(scip_exprdata_monomial)[1]
array(scip_exprgraph::SCIP_EXPRGRAPH_t) = scip_exprgraph.array_ptr_scip_exprgraph
pointer(scip_exprgraph::SCIP_EXPRGRAPH_t) = array(scip_exprgraph)[1]
array(scip_exprgraphnode::SCIP_EXPRGRAPHNODE_t) = scip_exprgraphnode.array_ptr_scip_exprgraphnode
pointer(scip_exprgraphnode::SCIP_EXPRGRAPHNODE_t) = array(scip_exprgraphnode)[1]
array(scip_exprintdata::SCIP_EXPRINTDATA_t) = scip_exprintdata.array_ptr_scip_exprintdata
pointer(scip_exprintdata::SCIP_EXPRINTDATA_t) = array(scip_exprintdata)[1]
array(scip_exprtree::SCIP_EXPRTREE_t) = scip_exprtree.array_ptr_scip_exprtree
pointer(scip_exprtree::SCIP_EXPRTREE_t) = array(scip_exprtree)[1]
array(scip_hashmap::SCIP_HASHMAP_t) = scip_hashmap.array_ptr_scip_hashmap
pointer(scip_hashmap::SCIP_HASHMAP_t) = array(scip_hashmap)[1]
array(scip_hashmaplist::SCIP_HASHMAPLIST_t) = scip_hashmaplist.array_ptr_scip_hashmaplist
pointer(scip_hashmaplist::SCIP_HASHMAPLIST_t) = array(scip_hashmaplist)[1]
array(scip_hashtable::SCIP_HASHTABLE_t) = scip_hashtable.array_ptr_scip_hashtable
pointer(scip_hashtable::SCIP_HASHTABLE_t) = array(scip_hashtable)[1]
array(scip_heur::SCIP_HEUR_t) = scip_heur.array_ptr_scip_heur
pointer(scip_heur::SCIP_HEUR_t) = array(scip_heur)[1]
array(scip_heurdata::SCIP_HEURDATA_t) = scip_heurdata.array_ptr_scip_heurdata
pointer(scip_heurdata::SCIP_HEURDATA_t) = array(scip_heurdata)[1]
array(scip_history::SCIP_HISTORY_t) = scip_history.array_ptr_scip_history
pointer(scip_history::SCIP_HISTORY_t) = array(scip_history)[1]
array(scip_holelist::SCIP_HOLELIST_t) = scip_holelist.array_ptr_scip_holelist
pointer(scip_holelist::SCIP_HOLELIST_t) = array(scip_holelist)[1]
array(scip_intarray::SCIP_INTARRAY_t) = scip_intarray.array_ptr_scip_intarray
pointer(scip_intarray::SCIP_INTARRAY_t) = array(scip_intarray)[1]
array(scip_interval::SCIP_INTERVAL_t) = scip_interval.array_ptr_scip_interval
pointer(scip_interval::SCIP_INTERVAL_t) = array(scip_interval)[1]
array(scip_lpi::SCIP_LPI_t) = scip_lpi.array_ptr_scip_lpi
pointer(scip_lpi::SCIP_LPI_t) = array(scip_lpi)[1]
array(scip_messagehdlr::SCIP_MESSAGEHDLR_t) = scip_messagehdlr.array_ptr_scip_messagehdlr
pointer(scip_messagehdlr::SCIP_MESSAGEHDLR_t) = array(scip_messagehdlr)[1]
array(scip_messagehdlrdata::SCIP_MESSAGEHDLRDATA_t) = scip_messagehdlrdata.array_ptr_scip_messagehdlrdata
pointer(scip_messagehdlrdata::SCIP_MESSAGEHDLRDATA_t) = array(scip_messagehdlrdata)[1]
array(scip_nlpi::SCIP_NLPI_t) = scip_nlpi.array_ptr_scip_nlpi
pointer(scip_nlpi::SCIP_NLPI_t) = array(scip_nlpi)[1]
array(scip_nlpiproblem::SCIP_NLPIPROBLEM_t) = scip_nlpiproblem.array_ptr_scip_nlpiproblem
pointer(scip_nlpiproblem::SCIP_NLPIPROBLEM_t) = array(scip_nlpiproblem)[1]
array(scip_nlpstatistics::SCIP_NLPSTATISTICS_t) = scip_nlpstatistics.array_ptr_scip_nlpstatistics
pointer(scip_nlpstatistics::SCIP_NLPSTATISTICS_t) = array(scip_nlpstatistics)[1]
array(scip_nlrow::SCIP_NLROW_t) = scip_nlrow.array_ptr_scip_nlrow
pointer(scip_nlrow::SCIP_NLROW_t) = array(scip_nlrow)[1]
array(scip_node::SCIP_NODE_t) = scip_node.array_ptr_scip_node
pointer(scip_node::SCIP_NODE_t) = array(scip_node)[1]
array(scip_nodesel::SCIP_NODESEL_t) = scip_nodesel.array_ptr_scip_nodesel
pointer(scip_nodesel::SCIP_NODESEL_t) = array(scip_nodesel)[1]
array(scip_nodeseldata::SCIP_NODESELDATA_t) = scip_nodeseldata.array_ptr_scip_nodeseldata
pointer(scip_nodeseldata::SCIP_NODESELDATA_t) = array(scip_nodeseldata)[1]
array(scip_param::SCIP_PARAM_t) = scip_param.array_ptr_scip_param
pointer(scip_param::SCIP_PARAM_t) = array(scip_param)[1]
array(scip_pqueue::SCIP_PQUEUE_t) = scip_pqueue.array_ptr_scip_pqueue
pointer(scip_pqueue::SCIP_PQUEUE_t) = array(scip_pqueue)[1]
array(scip_presol::SCIP_PRESOL_t) = scip_presol.array_ptr_scip_presol
pointer(scip_presol::SCIP_PRESOL_t) = array(scip_presol)[1]
array(scip_presoldata::SCIP_PRESOLDATA_t) = scip_presoldata.array_ptr_scip_presoldata
pointer(scip_presoldata::SCIP_PRESOLDATA_t) = array(scip_presoldata)[1]
array(scip_pricer::SCIP_PRICER_t) = scip_pricer.array_ptr_scip_pricer
pointer(scip_pricer::SCIP_PRICER_t) = array(scip_pricer)[1]
array(scip_pricerdata::SCIP_PRICERDATA_t) = scip_pricerdata.array_ptr_scip_pricerdata
pointer(scip_pricerdata::SCIP_PRICERDATA_t) = array(scip_pricerdata)[1]
array(scip_probdata::SCIP_PROBDATA_t) = scip_probdata.array_ptr_scip_probdata
pointer(scip_probdata::SCIP_PROBDATA_t) = array(scip_probdata)[1]
array(scip_profile::SCIP_PROFILE_t) = scip_profile.array_ptr_scip_profile
pointer(scip_profile::SCIP_PROFILE_t) = array(scip_profile)[1]
array(scip_prop::SCIP_PROP_t) = scip_prop.array_ptr_scip_prop
pointer(scip_prop::SCIP_PROP_t) = array(scip_prop)[1]
array(scip_propdata::SCIP_PROPDATA_t) = scip_propdata.array_ptr_scip_propdata
pointer(scip_propdata::SCIP_PROPDATA_t) = array(scip_propdata)[1]
array(scip_ptrarray::SCIP_PTRARRAY_t) = scip_ptrarray.array_ptr_scip_ptrarray
pointer(scip_ptrarray::SCIP_PTRARRAY_t) = array(scip_ptrarray)[1]
array(scip_quadelem::SCIP_QUADELEM_t) = scip_quadelem.array_ptr_scip_quadelem
pointer(scip_quadelem::SCIP_QUADELEM_t) = array(scip_quadelem)[1]
array(scip_quadvarterm::SCIP_QUADVARTERM_t) = scip_quadvarterm.array_ptr_scip_quadvarterm
pointer(scip_quadvarterm::SCIP_QUADVARTERM_t) = array(scip_quadvarterm)[1]
array(scip_queue::SCIP_QUEUE_t) = scip_queue.array_ptr_scip_queue
pointer(scip_queue::SCIP_QUEUE_t) = array(scip_queue)[1]
array(scip_reader::SCIP_READER_t) = scip_reader.array_ptr_scip_reader
pointer(scip_reader::SCIP_READER_t) = array(scip_reader)[1]
array(scip_readerdata::SCIP_READERDATA_t) = scip_readerdata.array_ptr_scip_readerdata
pointer(scip_readerdata::SCIP_READERDATA_t) = array(scip_readerdata)[1]
array(scip_real::SCIP_Real_t) = scip_real.array_ptr_scip_real
pointer(scip_real::SCIP_Real_t) = array(scip_real)[1]
array(scip_realarray::SCIP_REALARRAY_t) = scip_realarray.array_ptr_scip_realarray
pointer(scip_realarray::SCIP_REALARRAY_t) = array(scip_realarray)[1]
array(scip_relax::SCIP_RELAX_t) = scip_relax.array_ptr_scip_relax
pointer(scip_relax::SCIP_RELAX_t) = array(scip_relax)[1]
array(scip_relaxdata::SCIP_RELAXDATA_t) = scip_relaxdata.array_ptr_scip_relaxdata
pointer(scip_relaxdata::SCIP_RELAXDATA_t) = array(scip_relaxdata)[1]
array(scip_resourceactivity::SCIP_RESOURCEACTIVITY_t) = scip_resourceactivity.array_ptr_scip_resourceactivity
pointer(scip_resourceactivity::SCIP_RESOURCEACTIVITY_t) = array(scip_resourceactivity)[1]
array(scip_result::SCIP_RESULT_t) = scip_result.array_ptr_scip_result
pointer(scip_result::SCIP_RESULT_t) = array(scip_result)[1]
array(scip_row::SCIP_ROW_t) = scip_row.array_ptr_scip_row
pointer(scip_row::SCIP_ROW_t) = array(scip_row)[1]
array(scip_sepa::SCIP_SEPA_t) = scip_sepa.array_ptr_scip_sepa
pointer(scip_sepa::SCIP_SEPA_t) = array(scip_sepa)[1]
array(scip_sepadata::SCIP_SEPADATA_t) = scip_sepadata.array_ptr_scip_sepadata
pointer(scip_sepadata::SCIP_SEPADATA_t) = array(scip_sepadata)[1]
array(scip_sol::SCIP_SOL_t) = scip_sol.array_ptr_scip_sol
pointer(scip_sol::SCIP_SOL_t) = array(scip_sol)[1]
array(scip_sparsesol::SCIP_SPARSESOL_t) = scip_sparsesol.array_ptr_scip_sparsesol
pointer(scip_sparsesol::SCIP_SPARSESOL_t) = array(scip_sparsesol)[1]
array(scip_valuehistory::SCIP_VALUEHISTORY_t) = scip_valuehistory.array_ptr_scip_valuehistory
pointer(scip_valuehistory::SCIP_VALUEHISTORY_t) = array(scip_valuehistory)[1]
array(scip_var::SCIP_VAR_t) = scip_var.array_ptr_scip_var
pointer(scip_var::SCIP_VAR_t) = array(scip_var)[1]
array(scip_vardata::SCIP_VARDATA_t) = scip_vardata.array_ptr_scip_vardata
pointer(scip_vardata::SCIP_VARDATA_t) = array(scip_vardata)[1]


text:
include("../intcode.jl")

tape = Tape(parse_data(readline(open("input.txt"))))

const NORTH = 1
const SOUTH = 2
const WEST = 3
const EAST = 4

const dirs = Dict(
    NORTH => [0, +1],
    SOUTH => [0, -1],
    WEST => [-1, 0],
    EAST => [+1, 0],
)

const UNEXPLORED = -1
const WALL = 0
const EMPTY = 1
const OXYGEN = 2

const chars = Dict(
    UNEXPLORED => ' ',
    WALL => '#',
    EMPTY => '.',
    OXYGEN => 'x',
)

function squeeze_map(areamap, pos)
    is_explored = areamap .!= UNEXPLORED
    exploredX = (1:size(areamap, 1))[maximum(is_explored, dims=2)[:]]
    exploredY = (1:size(areamap, 2))[maximum(is_explored, dims=1)[:]]
    lowX, hiX = minimum(exploredX), maximum(exploredX)
    lowY, hiY = minimum(exploredY), maximum(exploredY)
    shifted_pos = pos - [lowX-1, lowY-1]
    areamap[lowX:hiX, lowY:hiY], shifted_pos
end

function plot_map(areamap, pos)
    for y=size(areamap, 2):-1:1
        for x=1:size(areamap, 1)
            if [x, y] == pos
                ch = (areamap[x, y] == OXYGEN) ? '@' : 'O'
            else
                ch = chars[areamap[x, y]]
            end
            print(ch)
        end
        println()
    end
end

function read_dir()
    valid_input = map(string, collect(keys(dirs)))
    while !((input = readline()) in valid_input) end
    parse(Int, input)
end

function steer!(tape)
    areamap = fill(-1, 1001, 1001)
    pos = [501, 501]
    areamap[pos...] = 1
    while true
        out = take!(tape.outch)
        if isa(out, AwaitingInput)
            print("Direction (N=1, S=2, W=3, E=4): ")
            direction = read_dir()
            put!(tape.inch, direction)
            output = take!(tape.outch)
            dir_vec = dirs[direction]
            newpos = pos + dir_vec
            areamap[newpos...] = output
            if output != WALL
                pos = newpos
            end

            plot_map(squeeze_map(areamap, pos)...)

        elseif isa(out, CodeFinished)
            println("FINISHED")
            return
        else
            println("OUTPUT: ", out)
        end
    end
end

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
# TODO, NOTE
# We should be doing each f1xfm **AFTER** each trial as been cycle averaged
using LinearAlgebra

export f1xfm, cycle_mean
# ============================================================================ #
@inline calcf1(cs::Vector{<:Complex}, d::Matrix{<:Real}) = vec(d' * cs)
# ---------------------------------------------------------------------------- #
@inline calcf1(cs::Vector{<:Complex}, d::Vector{<:Real}) = dot(d, cs)
# ---------------------------------------------------------------------------- #
@inline f1basis(npt::Int) = exp.(-im .* 2.0 .* pi .* range(0.0, stop=1.0, length=npt))
# ============================================================================ #
@inline f1xfm(d::Vector{<:Real}) = calcf1(f1basis(length(d)), d) * (2.0 / length(d))
# ============================================================================ #
function f1xfm(d::Array{<:Real}, f1::AbstractFloat, dur::AbstractFloat)
    bpc = floor(Int, size(d, 1) / (dur * f1))
    tmp = cycle_mean(d, bpc)
    return f1xfm(tmp)
end
# ============================================================================ #
function f1xfm(d::Vector{Matrix{<:Real}}, f1::Real, dur::Real)

    siz = size(d[1], 1)
    check_size(x::Matrix{<:Real}) = size(x, 1) == siz

    out = [Vector{Float64}(undef, size(d[k], 2)) for k in 1:length(d)]

    if !all(check_size, d[2:end])
        # psth matricies have different number of bins so the f1basis must be
        # caclulated for each matrix
        @inbounds for k = 1:length(d)
            out[k] = f1xfm(d[k], f1, dur)
        end
    else
        # our sinusodial basis can be reused so pre-allocate
        bpc = floor(Int, size(d[1], 1) / (dur * f1))
        cs = f1basis(bpc)

        @inbounds for k = 1:length(d)
            out[k] = calcf1(cs, cycle_mean(d[k], bpc)) * (4.0 / bpc)
        end
    end

    return out
end
# ============================================================================ #
function f1xfm(d::Vector{Matrix{<:Real}}, f1::Vector{<:Real}, dur::Real)

    length(f1) != length(d) && error("Number of temporal frequencies and number of trial groups *MUST* match")

    out = Vector{Vector{Float64}}(length(d))

    @inbounds for k = 1:length(d)
        out[k] = f1xfm(d[k], f1[k], dur)
    end

    return out
end
# ============================================================================ #
function cycle_pad(npt::Int, bpc::Integer)
    ncycle = floor(Int64, npt / bpc)
    if npt % bpc > 0
        ncycle += 1
        npad = (bpc * ncycle) - npt
    else
        npad = 0
    end
    return npad, ncycle
end
# ---------------------------------------------------------------------------- #
function cycle_mean(d::Vector{T}, bpc::Integer) where T<:Real
    pad, ncycle = cycle_pad(length(d), bpc)
    den = [fill(ncycle, bpc-pad); fill(ncycle-1, pad)]
    return vec(sum(reshape(cat(d, zeros(T, pad), dims=1), bpc, ncycle), dims=2)) ./ den
end
# ---------------------------------------------------------------------------- #
function cycle_mean(d::Matrix{T}, bpc::Integer) where T<:Real
    pad, ncycle = cycle_pad(size(d, 1), bpc)
    ntrial = size(d, 2)

    den = [fill(ncycle, bpc-pad); fill(ncycle-1, pad)]
    tmp = cat(d, zeros(pad, ntrial), dims=1)

    out = Matrix{T}(undef, bpc, ntrial)
    @inbounds @fastmath for k in 1:ntrial
        out[:,k] = sum(reshape(tmp[:,k], bpc, ncycle), dims=2) ./ den
    end

    return out
end
# ============================================================================ #


text:
using JuMP, JuMPChance
using Ipopt
using Gurobi
using Distributions

function Solve_DD_DR_CCLinDist(buses, lines, generators, wc_variances)
info("Starting Explicitly Implemented CC LinDist Model")

bus_set = collect(keys(buses))
line_set = collect(keys(lines))
gen_set = collect(keys(generators))
gen_bus = []
for g in keys(generators)
    push!(gen_bus, generators[g].bus_idx)
end
root_bus = 0
for b in keys(buses)
    if buses[b].is_root
        root_bus = b
    end
end

lines_to = Dict()
for l in keys(lines)
    lines_to[lines[l].to_node] = lines[l]
end

v_root = 1

var_sum_P = sum(wc_variances)
var_sum_Q = sum(wc_variances[b] * buses[b].tanphi^2 for b in bus_set)
var_P = wc_variances
var_Q = [wc_variances[b] * buses[b].tanphi^2 for b in bus_set]

m = Model(solver=IpoptSolver(print_level=3))
# m = Model(solver=GurobiSolver())

@variable(m, v[bus_set] >= 0) #variable for voltage square
@variable(m, fp[bus_set]) #variable for active power flow
@variable(m, fq[bus_set]) # variable for reactive power flow
@variable(m, gp[bus_set]) #variable for active power generation
@variable(m, gq[bus_set]) #variable for reactive power generation
@variable(m, α[bus_set] >= 0) # variable to distribute load deviations (affine frequency control)

@objective(m, Min, sum((gp[g]^2 + α[g]^2*var_sum_P) * buses[g].generator.cost for g in gen_bus))

# Control variable
@constraint(m, sum(α) == 1)

#Voltage and flow constraints for root node
@constraint(m, v[root_bus] == v_root)
@constraint(m, fp[root_bus] == 0)
@constraint(m, fq[root_bus] == 0)

for b in setdiff(bus_set, gen_bus)
# All buses without generation
    @constraint(m, α[b] == 0)
    @constraint(m, gp[b] == 0)
    @constraint(m, gq[b] == 0)
end

for b in bus_set
# All buses
    @constraint(m, buses[b].d_P - gp[b] + sum(fp[k] for k in buses[b].children) == fp[b])
    @constraint(m, buses[b].d_Q - gq[b] + sum(fq[k] for k in buses[b].children) == fq[b])
end

for b in setdiff(bus_set, [root_bus])
# All buses without root
    b_ancestor = buses[b].ancestor[1]
    @constraint(m, v[b] == v[b_ancestor] - 2*(lines_to[b].r * fp[b] + lines_to[b].x * fq[b]))
    @constraint(m, fp[b]^2 + fq[b]^2 <= lines_to[b].s_max^2)
end

# CHANCE CONSTRAINTS:
for b in gen_bus
# All buses with generation
    @constraint(m, gp[b] + z_g * α[b]*sqrt(var_sum_P) <= buses[b].generator.g_P_max)
    @constraint(m, -gp[b] + z_g * α[b]*sqrt(var_sum_P) <= 0)

    @constraint(m, gq[b] + z_g * α[b]*sqrt(var_sum_Q) <= buses[b].generator.g_Q_max)
    @constraint(m, -gq[b] + z_g * α[b]*sqrt(var_sum_Q) <= buses[b].generator.g_Q_max)
end

for b in setdiff(bus_set, [root_bus])
# All buses without root
    l = lines_to[b]
    line_arr = [lines[l] for l in line_set]
    @NLconstraint(m, v[b] + z_v * 2 * sqrt(sum(rPTDF[l.index, b] * l.r * sum(rPTDF[l.index, j] * (var_P[j] + α[j]^2*var_sum_P) for j in bus_set) for l in line_arr) +
                                           sum(rPTDF[l.index, b] * l.x * sum(rPTDF[l.index, j] * (var_Q[j] + α[j]^2*var_sum_Q) for j in bus_set) for l in line_arr))
                      <= buses[b].v_max)
    @NLconstraint(m, -v[b] + z_v * 2 * sqrt(sum(rPTDF[l.index, b] * l.r * sum(rPTDF[l.index, j] * (var_P[j] + α[j]^2*var_sum_P) for j in bus_set) for l in line_arr) +
                                            sum(rPTDF[l.index, b] * l.x * sum(rPTDF[l.index, j] * (var_Q[j] + α[j]^2*var_sum_Q) for j in bus_set) for l in line_arr))
                      <= -buses[b].v_min)
end

solve(m)
objective = getobjectivevalue(m)
# solvetime = getsolvetime(m) # Probably not available for Ipopt
solvetime = 0

info("Objective = $(objective))")

# Prepare Results
bus_results = DataFrame(bus = Any[], d_P = Any[], d_Q = Any[], cosphi=Any[], tanphi = [],  g_P = Any[], g_Q = Any[], g_cost=Any[], v_squared = Any[], α=Any[])
line_results = DataFrame(line = Any[], from = Any[], to = Any[], f_P = Any[], f_Q = Any[], a = Any[])
for b in bus_set
    cost = 0
    if in(b, gen_bus)
        cost = buses[b].generator.cost
    end

    # Numerical smoothing
    α_res = getvalue(α[b])
    α_res < 1e-10 ? α_res = 0 : nothing

    res = [b, buses[b].d_P, buses[b].d_Q, buses[b].cosphi, buses[b].tanphi,  getvalue(gp[b]), getvalue(gq[b]), cost, getvalue(v[b]), α_res]
    push!(bus_results, res)

    # Calc square current on line
    a_res = 0
    fp_res = getvalue(fp[b])
    fq_res = getvalue(fq[b])
    if b != root_bus
        v_res = getvalue(v[buses[b].ancestor[1]])
        a_res = (fp_res^2 + fq_res^2)/v_res^2
        lres = [lines_to[b].index, buses[b].ancestor[1], b, fp_res, fq_res, a_res]
        push!(line_results, lres)
    end
end

sort!(bus_results, cols=:bus)
sort!(line_results, cols=:to)

return objective, bus_results, line_results, solvetime
end


text:
function ueberpruefeberechneteSchnittpunkte(Sx1,Sy1,Sx2,Sy2,x0,y0,N, JJ,nx,ny,x1,y1,theta,X,Y)

vx = cosd(theta);
vy = sind(theta);
Dim = nx;
Winkel = 1;
Nullzeile = 0;
theta = theta*pi/180;
laengenmatrix, SchnittVektorsort,Nullzeile = Schnittpunkte(Dim,Winkel,X,Y,x0,y0,x1,y1,N,JJ,nx,ny,vx,vy,theta,Nullzeile);

if almostequal(Sx1,round(SchnittVektorsort[2,1];digits=2)) && almostequal(Sy1,round(SchnittVektorsort[2,2];digits=2))
    println("Ihr habt den ersten Schnittpunkt korrekt berechnet!")
else
    println("Euer erster Schnittpunkt stimmt noch nicht.")
end
if almostequal(Sx2,round(SchnittVektorsort[3,1];digits=2)) && almostequal(Sy2,round(SchnittVektorsort[3,2];digits=2))
    println("Ihr habt den zweiten Schnittpunkt korrekt berechnet!")
else
    println("Euer zweiter Schnittpunkt stimmt noch nicht.")
end

end