text:
lookup(pix::Gray{<:Normed}, table::Vector{<:Gray{<:Normed}}) = table[pix.val.i+1]

const CINE_LUT = (Gray{N0f16}.([2, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 17, 18,
              19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 33,
              34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 48,
              49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 63,
              64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79,
              79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94,
              94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109,
              110, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124,
              125, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 137, 138,
              139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154,
              156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 167, 168, 169, 170, 171, 172,
              173, 175, 176, 177, 178, 179, 181, 182, 183, 184, 186, 187, 188, 189, 191, 192,
              193, 194, 196, 197, 198, 200, 201, 202, 204, 205, 206, 208, 209, 210, 212, 213,
              215, 216, 217, 219, 220, 222, 223, 225, 226, 227, 229, 230, 232, 233, 235, 236,
              238, 239, 241, 242, 244, 245, 247, 249, 250, 252, 253, 255, 257, 258, 260, 261,
              263, 265, 266, 268, 270, 271, 273, 275, 276, 278, 280, 281, 283, 285, 287, 288,
              290, 292, 294, 295, 297, 299, 301, 302, 304, 306, 308, 310, 312, 313, 315, 317,
              319, 321, 323, 325, 327, 328, 330, 332, 334, 336, 338, 340, 342, 344, 346, 348,
              350, 352, 354, 356, 358, 360, 362, 364, 366, 368, 370, 372, 374, 377, 379, 381,
              383, 385, 387, 389, 391, 394, 396, 398, 400, 402, 404, 407, 409, 411, 413, 416,
              418, 420, 422, 425, 427, 429, 431, 434, 436, 438, 441, 443, 445, 448, 450, 452,
              455, 457, 459, 462, 464, 467, 469, 472, 474, 476, 479, 481, 484, 486, 489, 491,
              494, 496, 499, 501, 504, 506, 509, 511, 514, 517, 519, 522, 524, 527, 529, 532,
              535, 537, 540, 543, 545, 548, 551, 553, 556, 559, 561, 564, 567, 570, 572, 575,
              578, 581, 583, 586, 589, 592, 594, 597, 600, 603, 606, 609, 611, 614, 617, 620,
              623, 626, 629, 632, 635, 637, 640, 643, 646, 649, 652, 655, 658, 661, 664, 667,
              670, 673, 676, 679, 682, 685, 688, 691, 694, 698, 701, 704, 707, 710, 713, 716,
              719, 722, 726, 729, 732, 735, 738, 742, 745, 748, 751, 754, 758, 761, 764, 767,
              771, 774, 777, 781, 784, 787, 790, 794, 797, 800, 804, 807, 811, 814, 817, 821,
              824, 828, 831, 834, 838, 841, 845, 848, 852, 855, 859, 862, 866, 869, 873, 876,
              880, 883, 887, 890, 894, 898, 901, 905, 908, 912, 916, 919, 923, 927, 930, 934,
              938, 941, 945, 949, 952, 956, 960, 964, 967, 971, 975, 979, 982, 986, 990, 994,
              998, 1001, 1005, 1009, 1013, 1017, 1021, 1025, 1028, 1032, 1036, 1040, 1044, 1048, 1052, 1056,
              1060, 1064, 1068, 1072, 1076, 1080, 1084, 1088, 1092, 1096, 1100, 1104, 1108, 1112, 1116, 1120,
              1124, 1128, 1132, 1137, 1141, 1145, 1149, 1153, 1157, 1162, 1166, 1170, 1174, 1178, 1183, 1187,
              1191, 1195, 1200, 1204, 1208, 1212, 1217, 1221, 1225, 1230, 1234, 1238, 1243, 1247, 1251, 1256,
              1260, 1264, 1269, 1273, 1278, 1282, 1287, 1291, 1295, 1300, 1304, 1309, 1313, 1318, 1322, 1327,
              1331, 1336, 1340, 1345, 1350, 1354, 1359, 1363, 1368, 1372, 1377, 1382, 1386, 1391, 1396, 1400,
              1405, 1410, 1414, 1419, 1424, 1428, 1433, 1438, 1443, 1447, 1452, 1457, 1462, 1466, 1471, 1476,
              1481, 1486, 1490, 1495, 1500, 1505, 1510, 1515, 1520, 1524, 1529, 1534, 1539, 1544, 1549, 1554,
              1559, 1564, 1569, 1574, 1579, 1584, 1589, 1594, 1599, 1604, 1609, 1614, 1619, 1624, 1629, 1634,
              1639, 1644, 1649, 1655, 1660, 1665, 1670, 1675, 1680, 1685, 1691, 1696, 1701, 1706, 1711, 1717,
              1722, 1727, 1732, 1738, 1743, 1748, 1753, 1759, 1764, 1769, 1775, 1780, 1785, 1791, 1796, 1801,
              1807, 1812, 1818, 1823, 1828, 1834, 1839, 1845, 1850, 1856, 1861, 1867, 1872, 1878, 1883, 1889,
              1894, 1900, 1905, 1911, 1916, 1922, 1927, 1933, 1939, 1944, 1950, 1956, 1961, 1967, 1972, 1978,
              1984, 1989, 1995, 2001, 2007, 2012, 2018, 2024, 2030, 2035, 2041, 2047, 2053, 2058, 2064, 2070,
              2076, 2082, 2087, 2093, 2099, 2105, 2111, 2117, 2123, 2129, 2135, 2140, 2146, 2152, 2158, 2164,
              2170, 2176, 2182, 2188, 2194, 2200, 2206, 2212, 2218, 2224, 2231, 2237, 2243, 2249, 2255, 2261,
              2267, 2273, 2279, 2286, 2292, 2298, 2304, 2310, 2317, 2323, 2329, 2335, 2341, 2348, 2354, 2360,
              2366, 2373, 2379, 2385, 2392, 2398, 2404, 2411, 2417, 2423, 2430, 2436, 2443, 2449, 2455, 2462,
              2468, 2475, 2481, 2488, 2494, 2501, 2507, 2514, 2520, 2527, 2533, 2540, 2546, 2553, 2559, 2566,
              2572, 2579, 2586, 2592, 2599, 2605, 2612, 2619, 2625, 2632, 2639, 2645, 2652, 2659, 2666, 2672,
              2679, 2686, 2693, 2699, 2706, 2713, 2720, 2726, 2733, 2740, 2747, 2754, 2761, 2767, 2774, 2781,
              2788, 2795, 2802, 2809, 2816, 2823, 2830, 2837, 2844, 2850, 2857, 2864, 2871, 2878, 2885, 2893,
              2900, 2907, 2914, 2921, 2928, 2935, 2942, 2949, 2956, 2963, 2970, 2978, 2985, 2992, 2999, 3006,
              3013, 3021, 3028, 3035, 3042, 3049, 3057, 3064, 3071, 3078, 3086, 3093, 3100, 3108, 3115, 3122,
              3130, 3137, 3144, 3152, 3159, 3166, 3174, 3181, 3189, 3196, 3204, 3211, 3218, 3226, 3233, 3241,
              3248, 3256, 3263, 3271, 3278, 3286, 3294, 3301, 3309, 3316, 3324, 3331, 3339, 3347, 3354, 3362,
              3370, 3377, 3385, 3393, 3400, 3408, 3416, 3423, 3431, 3439, 3447, 3454, 3462, 3470, 3478, 3486,
              3493, 3501, 3509, 3517, 3525, 3533, 3540, 3548, 3556, 3564, 3572, 3580, 3588, 3596, 3604, 3612,
              3620, 3628, 3636, 3644, 3652, 3660, 3668, 3676, 3684, 3692, 3700, 3708, 3716, 3724, 3732, 3740,
              3749, 3757, 3765, 3773, 3781, 3789, 3798, 3806, 3814, 3822, 3830, 3839, 3847, 3855, 3863, 3872,
              3880, 3888, 3897, 3905, 3913, 3922, 3930, 3938, 3947, 3955, 3963, 3972, 3980, 3989, 3997, 4006,
              4014, 4022, 4031, 4039, 4048, 4056, 4064, 4095, 4095, 4095, 4095, 4095, 4095, 4095, 4095, 4095]./(2^12)))


text:
mutable struct PPHMCTune <: SamplerTune
    logfgrad::Union{Function,Missing}
    epsilon::Float64
    delta::Float64
    nleap::Int64
    m::Int64
    randomization::Bool
    adapter::Float64
    moves::Vector{Int}

    PPHMCTune() = new()

    function PPHMCTune(x::Vector{T}, logfgrad::Union{Function,Missing}, epsilon::Float64, nleap::Int64, delta::Float64=0.003) where T <: GeneralNode
        new(logfgrad, epsilon, delta, nleap, 0, true, 0.4,Int[])
    end
end

#PPHMCTune(x::Vector{T}, logfgrad::Function, epsilon::Float64, nleap::Int64, delta::Float64=0.003) where T <: GeneralNode =
#    PPHMCTune(x, logfgrad, epsilon, nleap, delta)

const PPHMCVariate = Sampler{PPHMCTune, T} where T<:GeneralNode

#################### Sampler Constructor ####################

"""
    PPHMC(params::ElementOrVector{Symbol}; args...)

Construct a `Sampler` object for PPHMC sampling. The Parameter is assumed to be
a tree.

Returns a `Sampler{PPHMCTune}` type object.

* params: stochastic node to be updated with the sampler.

* args...: additional keyword arguments to be passed to the PNUTSVariate constructor.
"""
function PPHMC(params::ElementOrVector{Symbol}, epsilon::Float64, nleap::Int64, delta::Float64)
    samplerfx = function (model::Model, block::Integer)
        block = SamplingBlock(model, block, true)

        f = let block = block
            (x, sz, ll, gr) -> mlogpdfgrad!(block, x, sz, ll, gr)
        end
        v = Sampler(block, f, epsilon, nleap, delta)

        sample!(v::PPHMCVariate, f, model.iter <= model.burnin, model.burnin)

        relist(block, v)
    end
    Sampler(params, samplerfx, PPHMCTune())
end

sample!(v::PPHMCVariate; args...) = sample!(v, v.tune.logfgrad, bi)

function sample!(v::PPHMCVariate, logfgrad::Function, adapt::Bool, bi::Int64)
    mt = v.value[1]
    nl = size(mt)[1] - 1
    delta = v.tune.delta
    v.tune.m += 1
    r = randn(nl)
    blv = get_branchlength_vector(mt)
    set_branchlength_vector!(mt, molifier.(blv, delta))
    lf, grad = logfgrad(mt, nl, true, true)
    s = Tree_HMC_State(deepcopy(mt), r, grad, lf)

    currH = hamiltonian(s)
    ovnni = 0

    fac = adapt ? (1-v.tune.adapter)^(1-v.tune.m*1/bi) : 1
    epsilon = v.tune.epsilon * fac
    nl = v.tune.randomization ? rand(1:v.tune.nleap) : v.tune.nleap

    for i in 1:nl
        nni = refraction!(s, epsilon, logfgrad, delta, nl)
        @show hamiltonian(s)
        ovnni += nni
    end

    push!(v.tune.moves, ovnni)
    propH = hamiltonian(s)
    @show propH, currH
    if log(rand()) < propH - currH
        v.value[1] = s.x
    end
    v
end

ext: jl

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
    # Functions and methods for cascade computation

    """
    `Cascade.computeSteps(scheme::Cascade.PhotonExcitationScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})`
        ... computes in turn all the requested transition amplitudes as well as PhotoExcitation.Line's, etc. for all pre-specified
            excitation steps of the cascade. When compared with standard excitation computations, however, the amount of output is
            largely reduced and often just printed into the summary file. A set of  data::Cascade.ExcitationData  is returned.
    """
    function computeSteps(scheme::Cascade.PhotonExcitationScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})
        linesE = PhotoExcitation.Line[]
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        nt = 0;   st = 0;   previousMeanEn = 0.
        for  step  in  stepList
            st = st + 1
            nc = length(step.initialMultiplet.levels) * length(step.finalMultiplet.levels)
            println("\n  $st) Perform $(string(step.process)) amplitude computations for up to $nc excitation lines (without selection rules): ")
            if  printSummary   println(iostream, "\n* $st) Perform $(string(step.process)) amplitude computations for " *
                                                 "up to $nc excitation lines (without selection rules): ")   end

            if  step.process == Basics.PhotoExc()
                newLines = PhotoExcitation.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.grid,
                                                               step.settings, output=true, printout=false)
                append!(linesE, newLines);    nt = length(linesE)
            else   error("Unsupported atomic process for excitation computations.")
            end
            println("     Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, giving now rise " *
                    "to a total of $nt $(string(step.process)) decay lines." )
            if  printSummary   println(iostream, "\n*    Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, " *
                                                 "giving now rise to a total of $nt $(string(step.process)) decay lines." )   end
        end
        #
        data = Cascade.ExcitationData(linesE)
    end

    """
    `Cascade.determineSteps(scheme::Cascade.PhotonExcitationScheme, comp::Cascade.Computation,
                            initialList::Array{Cascade.Block,1}, excitedList::Array{Cascade.Block,1})`
        ... determines all step::Cascade.Step's that need to be computed for this decay cascade. It cycles through all processes of the given
            decay scheme and selects all pairs of blocks due to the selected cascade approach. It checks that at least on pair of levels
            supports a `photo-excitation' within the step. A stepList::Array{Cascade.Step,1} is returned, and for which subsequently all
            required transition amplitudes and rates/cross sections are computed.
    """
    function determineSteps(scheme::Cascade.PhotonExcitationScheme, comp::Cascade.Computation,
                            initialList::Array{Cascade.Block,1}, excitedList::Array{Cascade.Block,1})
        stepList = Cascade.Step[]
        if  comp.approach  in  [Cascade.AverageSCA(), Cascade.SCA()]
            for  initialBlock in initialList
                for  excitedBlock in excitedList
                    for  process  in  comp.scheme.processes
                        if      process == Basics.PhotoExc()
                            if  initialBlock.NoElectrons == excitedBlock.NoElectrons
                                settings = PhotoExcitation.Settings([E1], [UseCoulomb, UseBabushkin], false, false, false, false, LineSelection(),
                                                                    0., 0., 1.0e6, Basics.ExpStokes() )
                                push!( stepList, Cascade.Step(process, settings, initialBlock.confs, excitedBlock.confs,
                                                                                initialBlock.multiplet, excitedBlock.multiplet) )
                            end
                        else    error("stop a")
                        end
                    end
                end
            end
            #
        else  error("Unsupported cascade approach.")
        end
        return( stepList )
    end

    """
    `Cascade.generateBlocks(scheme::Cascade.PhotonExcitationScheme, comp::Cascade.Computation, confs::Array{Configuration,1}; printout::Bool=true)`
        ... generate all block::Cascade.Block's, that need to be computed for this excitation cascade, and compute also the corresponding multiplets.
            The different cascade approches realized different strategies how these block are selected and computed. A blockList::Array{Cascade.Block,1}
            is returned.
    """
    function generateBlocks(scheme::Cascade.PhotonExcitationScheme, comp::Cascade.Computation, confs::Array{Configuration,1}; printout::Bool=true)
        blockList = Cascade.Block[]
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        if    comp.approach == AverageSCA()
            if  printout
            println("\n* Generate blocks for excitation computations:")
            println("\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println("    + orbitals are generated independently for each block for a Dirac-Fock-Slater potential; ")
            println("    + all blocks (multiplets) are generated from single-CSF levels and without any configuration mixing even in the SC; ")
            println("    + only E1 excitations are considered. \n")
            if  printSummary
            println(iostream, "\n* Generate blocks for excitation computations:")
            println(iostream, "\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println(iostream, "    + orbitals are generated independently for each block for a Dirac-Fock-Slater potential; ")
            println(iostream, "    + all blocks (multiplets) are generated from single-CSF levels and without any configuration mixing even in the SC; ")
            println(iostream, "    + only E1 excitations are considered. \n")
            end
            end     # printout
            #
            for  confa  in confs
                print("  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... ")
                if  printSummary   println(iostream, "\n*  Multiplet computations for $(string(confa)[1:end])   with $(confa.NoElectrons) electrons ... ")   end
                    basis     = Basics.performSCF([confa], comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                    multiplet = Basics.perform("computation: mutiplet from orbitals, no CI, CSF diagonal", [confa],  basis.orbitals,
                                               comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
            end
        elseif    comp.approach == SCA()  error("Not yet implemented.")
        else  error("Unsupported cascade approach.")
        end

        return( blockList )
    end

    """
    `Cascade.generateConfigurationsForPhotoexcitation(multiplets::Array{Multiplet,1},  scheme::PhotonExcitationScheme, nm::Nuclear.Model)`
        ... generates all possible (photo-excited) configurations with upto  scheme.NoExcitations of displacements of electron from
            scheme.excitationFromShells  to  scheme.excitationToShells  with regard to the given mutiplets.
            A Tuple(initialConfList::Array{Configuration,1}, confList::Array{Configuration,1}) is returned.
    """
    function generateConfigurationsForPhotoexcitation(multiplets::Array{Multiplet,1},  scheme::PhotonExcitationScheme, nm::Nuclear.Model)
        # Determine all (reference) configurations from multiplets and generate the 'excited' configurations due to the specificed excitations
        initialConfList = Configuration[]
        for mp  in  multiplets
            confList = Basics.extractNonrelativisticConfigurations(mp.levels[1].basis)
            for  conf in confList   if  conf in initialConfList   nothing   else   push!(initialConfList, conf)      end      end
        end
        blockConfList = Basics.generateConfigurations(initialConfList, scheme.excitationFromShells, scheme.excitationToShells, scheme.NoExcitations)
        # Exclude configurations with too low or too high mean energies as well as those that are parity forbidden for the given multipoles
        hasPlus = false;   hasMinus = false
        for  conf  in  initialConfList
            if  Basics.determineParity(conf) == Basics.plus   hasPlus = true    else   hasMinus = true      end
        end
        en     = Float64[];   for conf in initialConfList    push!(en, -Semiempirical.estimate("binding energy", round(Int64, nm.Z), conf))   end
        maxen  = maximum(en);    minen  = minimum(en);
        println(">>> initial configuration(s) have energies from $minen  to  $maxen  [a.u.].")
        #
        newBlockConfList = Configuration[]
        for  conf  in  blockConfList    meanEnergy = -Semiempirical.estimate("binding energy", round(Int64, nm.Z), conf)
            if  minen + scheme.minPhotonEnergy  <= meanEnergy <= maxen + scheme.maxPhotonEnergy
                if      hasPlus  &&  hasMinus                                                                        push!(newBlockConfList, conf)
                elseif  M1 in scheme.multipoles  ||  E2 in scheme.multipoles  ||  M2 in scheme.multipoles            push!(newBlockConfList, conf)
                elseif  hasPlus  &&  E1 in scheme.multipoles  &&   Basics.determineParity(conf) == Basics.minus      push!(newBlockConfList, conf)
                elseif  hasMinus &&  E1 in scheme.multipoles  &&   Basics.determineParity(conf) == Basics.plus       push!(newBlockConfList, conf)
                else    println(">>> exclude $conf because of parity reasons.")
                end
            else        println(">>> exclude $conf with energy $meanEnergy [a.u.] because of energy reasons.")
            end
        end

        return( (initialConfList, newBlockConfList)  )
    end

    """
    `Cascade.perform(scheme::PhotonExcitationScheme, comp::Cascade.Computation)`
        ... to set-up and perform a photo-excitation computation that starts from a given set of initial configurations xor initial multiplets
            and comprises (various) photoexcitation processes into configurations with up-to NoExcitations single-electron excitations with
            regard to the initial multiplets. The results of these excitation are comprised into (output) data::PhotoExcData, while these
            data are only printed during the generation and nothing is returned.

    `Cascade.perform(scheme::PhotonExcitationScheme, comp::Cascade.Computation; output=true, outputToFile::Bool=true)`
        ... to perform the same but to return the complete output in a dictionary that is written to disk and can be used in subsequent
            cascade simulation. The particular output depends on the specifications of the cascade.
    """
    function perform(scheme::PhotonExcitationScheme, comp::Cascade.Computation; output::Bool=false, outputToFile::Bool=true)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        # Perform the SCF and CI computation for the intial-state multiplets if initial configurations are given
        if  comp.initialConfigs != Configuration[]
            basis      = Basics.performSCF(comp.initialConfigs, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplet  = Basics.performCI(basis, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplets = [Multiplet("initial states", multiplet.levels)]
        else
            multiplets = comp.initialMultiplets
        end
        # Print out initial configurations and levels
        Cascade.displayLevels(stdout, multiplets, sa="initial ")
        if  printSummary   Cascade.displayLevels(iostream, multiplets, sa="initial ")                            end
        #
        # Generate subsequent cascade configurations as well as display and group them together
        wa  = Cascade.generateConfigurationsForPhotoexcitation(multiplets, comp.scheme, comp.nuclearModel)
        wb1 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[1], sa="(initial part of the) photo-excited ")
        wb2 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[2], sa="(generated part of the) photo-excited ")
        #
        # Determine first all configuration 'blocks' and from them the individual steps of the cascade
        wc1 = Cascade.generateBlocks(scheme, comp::Cascade.Computation, wb1)
        wc2 = Cascade.generateBlocks(scheme, comp::Cascade.Computation, wb2, printout=false)
        Cascade.displayBlocks(stdout, wc1, sa="for the (initial part of the) excited cascade ");
        Cascade.displayBlocks(stdout, wc2, sa="for the (generated part of the) excited cascade ")
        if  printSummary   Cascade.displayBlocks(iostream, wc1, sa="for the (initial part of the) excited cascade ")
                           Cascade.displayBlocks(iostream, wc2, sa="for the (generated part of the) excited cascade ")    end
        # Determine, modify and compute the transition data for all steps, ie. the PhotoIonization.Line's, etc.
        gMultiplets = Multiplet[];     for block in wc2  push!(gMultiplets, block.multiplet)    end
        we = Cascade.determineSteps(scheme, comp, wc1, wc2)
        Cascade.displaySteps(stdout, we, sa="excited ")
        if  printSummary   Cascade.displaySteps(iostream, we, sa="ionizing ")    end
        wf   = Cascade.modifySteps(we)
        data = Cascade.computeSteps(scheme, comp, wf)
        if output
            results = Base.merge( results, Dict("name"                       => comp.name) )
            results = Base.merge( results, Dict("cascade scheme"             => comp.scheme) )
            results = Base.merge( results, Dict("initial multiplets:"        => multiplets) )
            results = Base.merge( results, Dict("generated multiplets:"      => gMultiplets) )
            results = Base.merge( results, Dict("photo-excited line data:"   => data) )
            #
            #  Write out the result to file to later continue with simulations on the cascade data
            filename = "zzz-cascade-excitation-computations-" * string(Dates.now())[1:13] * ".jld"
            println("\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                    "\n   results = JLD.load(''$filename'')    ... to load the results back from file.")
            if  printSummary   println(iostream, "\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                                                 "\n   results = JLD.load(''$filename'')    ... to load the results back from file." )      end
            JLD2.@save filename results
        end
        ## return( results )
        return( results )
    end


text:
using Oceananigans.Architectures

using Oceananigans.Grids: topology, validate_tupled_argument

struct MultiCPU{G, R, I, ρ, C, γ} <: AbstractCPUArchitecture
    distributed_grid :: G
          local_rank :: R
         local_index :: I
               ranks :: ρ
        connectivity :: C
        communicator :: γ
end

struct MultiGPU{G, R, I, ρ, C, γ} <: AbstractGPUArchitecture
    distributed_grid :: G
          local_rank :: R
         local_index :: I
               ranks :: ρ
        connectivity :: C
        communicator :: γ
end

const AbstractMultiArchitecture = Union{MultiCPU, MultiGPU}

child_architecture(::MultiCPU) = CPU()
child_architecture(::CPU) = CPU()

child_architecture(::MultiGPU) = GPU()
child_architecture(::GPU) = GPU()

#####
##### Converting between index and MPI rank taking k as the fast index
#####

index2rank(i, j, k, Rx, Ry, Rz) = (i-1)*Ry*Rz + (j-1)*Rz + (k-1)

function rank2index(r, Rx, Ry, Rz)
    i = div(r, Ry*Rz)
    r -= i*Ry*Rz
    j = div(r, Rz)
    k = mod(r, Rz)
    return i+1, j+1, k+1  # 1-based Julia
end

#####
##### Rank connectivity graph
#####

struct RankConnectivity{E, W, N, S, T, B}
      east :: E
      west :: W
     north :: N
     south :: S
       top :: T
    bottom :: B
end

RankConnectivity(; east, west, north, south, top, bottom) =
    RankConnectivity(east, west, north, south, top, bottom)

function increment_index(i, R, topo)
    R == 1 && return nothing
    if i+1 > R
        if topo == Periodic
            return 1
        else
            return nothing
        end
    else
        return i+1
    end
end

function decrement_index(i, R, topo)
    R == 1 && return nothing
    if i-1 < 1
        if topo == Periodic
            return R
        else
            return nothing
        end
    else
        return i-1
    end
end

function RankConnectivity(model_index, ranks, topology)
    i, j, k = model_index
    Rx, Ry, Rz = ranks
    TX, TY, TZ = topology

    i_east  = increment_index(i, Rx, TX)
    i_west  = decrement_index(i, Rx, TX)
    j_north = increment_index(j, Ry, TY)
    j_south = decrement_index(j, Ry, TY)
    k_top   = increment_index(k, Rz, TZ)
    k_bot   = decrement_index(k, Rz, TZ)

    r_east  = isnothing(i_east)  ? nothing : index2rank(i_east, j, k, Rx, Ry, Rz)
    r_west  = isnothing(i_west)  ? nothing : index2rank(i_west, j, k, Rx, Ry, Rz)
    r_north = isnothing(j_north) ? nothing : index2rank(i, j_north, k, Rx, Ry, Rz)
    r_south = isnothing(j_south) ? nothing : index2rank(i, j_south, k, Rx, Ry, Rz)
    r_top   = isnothing(k_top)   ? nothing : index2rank(i, j, k_top, Rx, Ry, Rz)
    r_bot   = isnothing(k_bot)   ? nothing : index2rank(i, j, k_bot, Rx, Ry, Rz)

    return RankConnectivity(east=r_east, west=r_west, north=r_north,
                            south=r_south, top=r_top, bottom=r_bot)
end

#####
##### Constructors
#####

function MultiArchitecture(MultiArch, grid, ranks, communicator)
    MPI.Initialized() || error("Must call MPI.Init() before constructing a MultiCPU.")

    validate_tupled_argument(ranks, Int, "ranks")

    Rx, Ry, Rz = ranks
    total_ranks = Rx*Ry*Rz

    mpi_ranks = MPI.Comm_size(communicator)
    local_rank   = MPI.Comm_rank(communicator)

    i, j, k = local_index = rank2index(local_rank, Rx, Ry, Rz)

    if total_ranks != mpi_ranks
        throw(ArgumentError("ranks=($Rx, $Ry, $Rz) [$total_ranks total] inconsistent " *
                            "with number of MPI ranks: $mpi_ranks."))
    end

    local_connectivity = RankConnectivity(local_index, ranks, topology(grid))

    return MultiArch(grid, local_rank, local_index, ranks, local_connectivity, communicator)
end

MultiCPU(; grid, ranks, communicator=MPI.COMM_WORLD) = MultiArchitecture(MultiCPU, grid, ranks, communicator)
MultiGPU(; grid, ranks, communicator=MPI.COMM_WORLD) = MultiArchitecture(MultiGPU, grid, ranks, communicator)

#####
##### Pretty printing
#####

function Base.show(io::IO, arch::MultiCPU)
    c = arch.connectivity
    print(io, "MultiCPU architecture (rank $(arch.local_rank)/$(prod(arch.ranks)-1)) [index $(arch.local_index) / $(arch.ranks)]\n",
              "└── connectivity:",
              isnothing(c.east) ? "" : " east=$(c.east)",
              isnothing(c.west) ? "" : " west=$(c.west)",
              isnothing(c.north) ? "" : " north=$(c.north)",
              isnothing(c.south) ? "" : " south=$(c.south)",
              isnothing(c.top) ? "" : " top=$(c.top)",
              isnothing(c.bottom) ? "" : " bottom=$(c.bottom)")
end


text:
# Risk measures
abstract type AbstractRiskMeasure end
# Expectation
struct Expectation <: AbstractRiskMeasure end
# Worst case
struct WorstCase <: AbstractRiskMeasure end
# Conditional value at risk
struct CVaR <: AbstractRiskMeasure
    β::Float64

    function CVaR(β::Float64)
        @assert 0. <= β <= 1. "β must be in [0,1]"
        return new(β)
    end
end

conditional_value_at_risk(support::Array{Float64,1}, probabilities::Array{Float64,1}, risk::WorstCase) = conditional_value_at_risk(support, probabilities, 0.)
conditional_value_at_risk(support::Array{Float64,1}, probabilities::Array{Float64,1}, risk::Expectation) = conditional_value_at_risk(support, probabilities, 1.)
conditional_value_at_risk(support::Array{Float64,1}, probabilities::Array{Float64,1}, risk::CVaR) = conditional_value_at_risk(support, probabilities, 1. - risk.β)

# From https://github.com/jaantollander/ConditionalValueAtRisk
function conditional_value_at_risk(support::Array{Float64,1}, probabilities::Array{Float64,1}, α::Float64)
    # Value at risk
    var = value_at_risk(support, probabilities, α)
    # Tail values
    if α == 0.
        return var
    else
        tail = support .< var
        return (sum(probabilities[tail] .* support[tail]) - (sum(probabilities[tail]) - α) * var) / α
    end
end

function value_at_risk(support::Array{Float64,1}, probabilities::Array{Float64,1}, α::Float64)
    i = findfirst(cumsum(probabilities[sortperm(support)]) .>= α)
    if i == nothing
        return sort(support)[end]
    else
        return sort(support)[i]
    end
end

# MILP functions
# Decisions
function add_investment_decisions!(m::Model, generations::Vector{AbstractGeneration})
    if !isempty(generations)
        na = length(generations)
        @variable(m, r_g[1:na])
    end
end
function add_investment_decisions!(m::Model, storages::Vector{AbstractStorage})
    if !isempty(storages)
        na = length(storages)
        @variable(m, r_sto[1:na])
    end
end
function add_investment_decisions!(m::Model, converters::Vector{AbstractConverter})
    if !isempty(converters)
        na = length(converters)
        @variable(m, r_c[1:na])
    end
end
function fix_investment_decisions!(m::Model, generations::Vector{Float64}, storages::Vector{Float64}, converters::Vector{Float64})
    # Generation
    if !isempty(generations)
        fix.(m[:r_g], generations)
    end
    # Storages
    if !isempty(storages)
        fix.(m[:r_sto], storages)
    end
    # Converters
    if !isempty(converters)
        fix.(m[:r_c], converters)
    end
end
function add_operation_decisions!(m::Model, demands::Vector{AbstractDemand}, nh::Int64, ns::Int64)
    if !isempty(demands)
        na = length(demands)
        @variables(m, begin
        p_d[1:nh, 1:ns, 1:na]
        end)
    end
end
function add_operation_decisions!(m::Model, generations::Vector{AbstractGeneration}, nh::Int64, ns::Int64)
    if !isempty(generations)
        na = length(generations)
        @variables(m, begin
        p_g[1:nh, 1:ns, 1:na]
        end)
    end
end
function add_operation_decisions!(m::Model, storages::Vector{AbstractStorage}, nh::Int64, ns::Int64)
    if !isempty(storages)
        na = length(storages)
        @variables(m, begin
        p_ch[1:nh, 1:ns, 1:na]   >= 0.
        p_dch[1:nh, 1:ns, 1:na]  >= 0.
        soc[1:nh+1, 1:ns, 1:na]
        end)
    end
end
function add_operation_decisions!(m::Model, converters::Vector{AbstractConverter}, nh::Int64, ns::Int64)
    if !isempty(converters)
        na = length(converters)
        @variable(m, p_c[1:nh, 1:ns, 1:na] >= 0.)
    end
end
function add_operation_decisions!(m::Model, grids::Vector{AbstractGrid}, nh::Int64, ns::Int64)
    if !isempty(grids)
        na = length(grids)
        @variables(m, begin
        p_in[1:nh, 1:ns, 1:na]   >= 0.
        p_out[1:nh, 1:ns, 1:na]  >= 0.
        end)
    end
end
# Investment bounds
function add_investment_constraints!(m::Model, generations::Vector{AbstractGeneration})
    if !isempty(generations)
        na = length(generations)
        @constraints(m, begin
        [a in 1:na], m[:r_g][a] >= generations[a].bounds.lb
        [a in 1:na], m[:r_g][a] <= generations[a].bounds.ub
        end)
    end
end
function add_investment_constraints!(m::Model, storages::Vector{AbstractStorage})
    if !isempty(storages)
        na = length(storages)
        @constraints(m, begin
        [a in 1:na], m[:r_sto][a] >= storages[a].bounds.lb
        [a in 1:na], m[:r_sto][a] <= storages[a].bounds.ub
        end)
    end
end
function add_investment_constraints!(m::Model, converters::Vector{AbstractConverter})
    if !isempty(converters)
        na = length(converters)
        @constraints(m, begin
        [a in 1:na], m[:r_c][a] >= converters[a].bounds.lb
        [a in 1:na], m[:r_c][a] <= converters[a].bounds.ub
        end)
    end
end
# Technical constraint
function add_technical_constraints!(m::Model, storages::Vector{AbstractStorage}, Δh::Int64, nh::Int64, ns::Int64)
    if !isempty(storages)
        na = length(storages)
        @constraints(m, begin
        # Power bounds
        [h in 1:nh, s in 1:ns, a in 1:na], m[:p_dch][h,s,a] <= storages[a].α_p_dch * m[:r_sto][a]
        [h in 1:nh, s in 1:ns, a in 1:na], m[:p_ch][h,s,a]  <= storages[a].α_p_ch * m[:r_sto][a]
        # SoC bounds
        [h in 1:nh+1, s in 1:ns, a in 1:na], m[:soc][h,s,a] <= storages[a].α_soc_max * m[:r_sto][a]
        [h in 1:nh+1, s in 1:ns, a in 1:na], m[:soc][h,s,a] >= storages[a].α_soc_min * m[:r_sto][a]
        # State dynamics
        [h in 1:nh, s in 1:ns, a in 1:na], m[:soc][h+1,s,a] == m[:soc][h,s,a] * (1. - storages[a].η_self * Δh) - (m[:p_dch][h,s,a] / storages[a].η_dch - m[:p_ch][h,s,a] * storages[a].η_ch) * Δh
        # Initial and final states
        soc_ini[s in 1:ns, a in 1:na], m[:soc][1,s,a] == storages[a].soc_ini * m[:r_sto][a]
        end)
    end
end
function add_technical_constraints!(m::Model, converters::Vector{AbstractConverter}, nh::Int64, ns::Int64)
    if !isempty(converters)
        na = length(converters)
        @constraints(m, begin
        # Power bounds
        [h in 1:nh, s in 1:ns, a in 1:na], m[:p_c][h,s,a]  <= m[:r_c][a]
        end)
    end
end
function add_technical_constraints!(m::Model, grids::Vector{AbstractGrid}, nh::Int64, ns::Int64)
    if !isempty(grids)
        na = length(grids)
        @constraints(m, begin
        # Power bounds
        [h in 1:nh, s in 1:ns, a in 1:na], m[:p_in][h,s,a]  <= grids[a].powerMax
        [h in 1:nh, s in 1:ns, a in 1:na], m[:p_out][h,s,a] <= grids[a].powerMax
        end)
    end
end
# Periodicity constraint
function add_periodicity_constraints!(m::Model, storages::Vector{AbstractStorage}, ns::Int64)
    # Storages
    if !isempty(storages)
        na = length(storages)
        @constraints(m, begin
        # Final states
        [s in 1:ns, a in 1:na], m[:soc][end,s,a]  >= m[:soc][1,s,a]
        end)
    end
end
# Power balance
function add_power_balance!(m::Model, mg::Microgrid, ω::Scenarios, type::DataType, nh::Int64, ns::Int64; ispnet::Bool=false)
    # !!! All the decision variables are defined positive !!!
    balance = AffExpr.(zeros(nh,ns))
    # Demands and generation
    if !ispnet
        for (k,a) in enumerate(mg.demands)
            if a.carrier isa type
                add_to_expression!.(balance, ω.demands[k].power[:,1,:])
            end
        end
        # Generation
        for (k,a) in enumerate(mg.generations)
            if a.carrier isa type
                add_to_expression!.(balance, .- m[:r_g][k] .* ω.generations[k].power[:,1,:])
            end
        end
    else
        for (k,a) in enumerate(mg.demands)
            if a.carrier isa type
                add_to_expression!.(balance, m[:p_d][:,:,k])
            end
        end
        # Generation
        for (k,a) in enumerate(mg.generations)
            if a.carrier isa type
                add_to_expression!.(balance, .- m[:p_g][:,:,k])
            end
        end
    end
    # Storages
    for (k,a) in enumerate(mg.storages)
        if a.carrier isa type
            add_to_expression!.(balance, m[:p_ch][:,:,k] .- m[:p_dch][:,:,k])
        end
    end
    # Converters
    for (k,a) in enumerate(mg.converters)
        if type == Electricity
            if a isa Heater
                add_to_expression!.(balance, m[:p_c][:,:,k])
            elseif a isa Electrolyzer
                add_to_expression!.(balance, m[:p_c][:,:,k])
            elseif a isa FuelCell
                add_to_expression!.(balance, .- m[:p_c][:,:,k])
            end
        elseif type == Heat
            if a isa Heater
                add_to_expression!.(balance, .- m[:p_c][:,:,k] * a.η_E_H)
            elseif a isa Electrolyzer
                add_to_expression!.(balance, .- m[:p_c][:,:,k] * a.η_E_H)
            elseif a isa FuelCell
                add_to_expression!.(balance, .- m[:p_c][:,:,k] * a.η_H2_H / a.η_H2_E)
            end
        elseif type == Hydrogen
            if a isa Electrolyzer
                add_to_expression!.(balance, .- m[:p_c][:,:,k] * a.η_E_H2)
            elseif a isa FuelCell
                add_to_expression!.(balance, m[:p_c][:,:,k] / a.η_H2_E)
            end
        end
    end
    # Grids
    for (k,a) in enumerate(mg.grids)
        if a.carrier isa type
            add_to_expression!.(balance, .- m[:p_in][:,:,k] + m[:p_out][:,:,k])
        end
    end
    # Energy balance constraint
    if type == Electricity
        @constraint(m, electricity, balance .<= 0.)
    elseif type == Heat
        @constraint(m, heat, balance .<= 0.)
    elseif type == Hydrogen
        @constraint(m, hydrogen, balance .== 0.)
    end
end
# Renewable share
function add_renewable_share!(m::Model, mg::Microgrid, ω::Scenarios, probabilities::Vector{Float64}, risk::AbstractRiskMeasure, nh::Int64, ns::Int64)
    total = zeros(ns)
    for (k,a) in enumerate(mg.demands)
        if a.carrier isa Electricity
            total .= total .+ sum(ω.demands[k].power[h,1,:] for h in 1:nh)
        elseif a.carrier isa Heat
            total .= total .+ sum(ω.demands[k].power[h,1,:] for h in 1:nh) ./ mg.converters[isin(mg.converters, Heater)[2]].η_E_H
        end
    end
    for (k,a) in enumerate(mg.grids)
        if a.carrier isa Electricity
            @expression(m, share[s in 1:ns], sum(m[:p_in][h,s,k] for h in 1:nh) - (1. - mg.parameters.renewable_share) * total[s])
        end
    end
    # Constraint according to CVaR
    @variables(m, begin
    ζ_s
    α_s[1:ns] >= 0.
    end)
    @constraints(m, begin
    [s in 1:ns], α_s[s] >= m[:share][s] - ζ_s
    ζ_s + 1 / (1 - beta(risk)) * sum(probabilities[s] * α_s[s] for s in 1:ns) <= 0.
    end)
end
# Objective
function add_design_objective!(m::Model, mg::Microgrid, ω::Scenarios, probabilities::Vector{Float64}, risk::AbstractRiskMeasure, nh::Int64, ns::Int64)
    # CAPEX
    capex = compute_capex(m, mg, ω)
    # OPEX
    opex = compute_opex(m, mg, ω, nh, ns)
    # Objective according to the CVaR
    @variables(m, begin
    ζ_o
    α_o[1:ns] >= 0.
    end)
    @constraint(m, [s in 1:ns], α_o[s] >= capex + opex[s] - ζ_o)
    @objective(m, Min, ζ_o + 1 / (1 - beta(risk)) * sum(probabilities[s] * α_o[s] for s in 1:ns))
end
# Capex
function compute_capex(m::Model, mg::Microgrid, ω::Scenarios)
    cost = AffExpr(0.)
    # Generations
    for (k,a) in enumerate(mg.generations)
        add_to_expression!(cost, Γ(mg.parameters.τ, a.lifetime) * ω.generations[k].cost[1] * m[:r_g][k])
    end
    # Storages
    for (k,a) in enumerate(mg.storages)
        add_to_expression!(cost, Γ(mg.parameters.τ, a.lifetime) * ω.storages[k].cost[1] * m[:r_sto][k])
    end
    # Converters
    for (k,a) in enumerate(mg.converters)
        add_to_expression!(cost, Γ(mg.parameters.τ, a.lifetime) * ω.converters[k].cost[1] * m[:r_c][k])
    end
    return cost
end
# Grids
function compute_opex(m::Model, mg::Microgrid, ω::Scenarios, nh::Int64, ns::Int64)
    cost = AffExpr.(zeros(ns))
    for (k,a) in enumerate(mg.grids)
        add_to_expression!.(cost, sum((m[:p_in][h,:,k] .* ω.grids[k].cost_in[h,1,:] .- m[:p_out][h,:] .* ω.grids[k].cost_out[h,1,:]) .* mg.parameters.Δh  for h in 1:nh))
    end
    return cost
end


text:
type PlotWindow
	renderwindow::RenderWindow
	graphs::Array{Graph}
	view::View
	event::Event
	last_mousepos::Vector2f
	xaxis::Axis
	yaxis::Axis
	ppu::Real
	followgraph::AnimatedGraph
	task
end

function PlotWindow(window::RenderWindow)
	width = get_size(window).x; height = get_size(window).y
	graphs = Graph[]

	view = View(Vector2f(0, 0), Vector2f(width, -height))
	event = Event()

	xaxis = Axis(view, true)
	yaxis = Axis(view, false)

	PlotWindow(window, graphs, view, event, Vector2f(0, 0), xaxis, yaxis, 20, AnimatedGraph(x -> 0), nothing)
end

function PlotWindow(plotname::AbstractString, width::Integer, height::Integer)
	settings = ContextSettings()
	settings.antialiasing_level = 4
	window = RenderWindow(VideoMode(width, height), plotname, settings, window_resize)
	# set_framerate_limit(window, 60)

	PlotWindow(window)
end

function add_graph(window::PlotWindow, graph::Graph)
	push!(window.graphs, graph)
end

function redraw(window::PlotWindow, fullredraw=false)
	center = get_center(window.view)
	left = Int(round(center.x - get_size(window.view).x/2))
	right = Int(round(center.x + get_size(window.view).x/2))

	if window.followgraph in window.graphs
		set_center(window.view, Vector2f(window.followgraph.xval - get_size(window.view).x/2 + 20, 0))
	end

	update(window.xaxis, window)
	update(window.yaxis, window)

	for j = 1:length(window.graphs)
		i = left
		if typeof(window.graphs[j]) == AnimatedGraph
			if window.graphs[j].xval > right
				advancex(window.graphs[j])
			else
				advance(window.graphs[j], window.ppu)
			end
		end
		while i <= right
			if typeof(window.graphs[j]) == AnimatedGraph
				if i > window.graphs[j].xval
					break
				end
			end
			# if typeof(window.graphs[j]) == StaticGraph
				if i % window.graphs[j].accuracy == 0
					if !fullredraw
						if haskey(window.graphs[j].points, i)
							add_point(window.graphs[j], i, window.graphs[j].points[i])
						else
							add_point(window.graphs[j], i, window.ppu)
						end
					else
						add_point(window.graphs[j], i, window.ppu)
					end
				end
				i += 1
			# else
			# 	break
			# end
		end
	end
	garbagecollect(window)
end

function check_input(window::PlotWindow)
	while pollevent(window.renderwindow, window.event)
		if get_type(window.event) == EventType.CLOSED
			SFML.close(window.renderwindow)
		end
		if get_type(window.event) == EventType.MOUSE_WHEEL_MOVED
			oldcenter = Vector2f(get_center(window.view).x / window.ppu, get_center(window.view).y / window.ppu)
			delta = get_mousewheel(window.event).delta
			window.ppu += delta / 10
			if window.ppu < 0
				window.ppu = 0.01
			end
			set_center(window.view, Vector2f(oldcenter.x*window.ppu, oldcenter.y*window.ppu))
			redraw(window, true)
		end
		if get_type(window.event) == EventType.MOUSE_BUTTON_PRESSED
			mouse_event = get_mousebutton(window.event)
			window.last_mousepos = Vector2f(mouse_event.x, mouse_event.y)
		end
		if get_type(window.event) == EventType.MOUSE_BUTTON_RELEASED
			garbagecollect(window)
		end
		if get_type(window.event) == EventType.MOUSE_MOVED
			if is_mouse_pressed(0)
				mouse_event = get_mousemove(window.event)
				mousepos = Vector2f(mouse_event.x, mouse_event.y)
				move(window.view, Vector2f(window.last_mousepos.x - mousepos.x, mousepos.y - window.last_mousepos.y))
				window.last_mousepos = mousepos
			end
		end
	end
end

function garbagecollect(window::PlotWindow)
	center = get_center(window.view)
	width = Vector2f(center.x - get_size(window.view).x/2, center.x + get_size(window.view).x/2)
	for i = 1:length(window.graphs)
		# if typeof(window.graphs[i]) == StaticGraph
			remove_points_outside(window.graphs[i], width)
		# end
	end
end

function isopen(window::PlotWindow)
	SFML.isopen(window.renderwindow)
end

function follow(window::PlotWindow, graph::Graph)
	window.followgraph = window.graphs[find(window.graphs .== graph)[1]]
	nothing
end

function unfollow(window::PlotWindow)
	window.followgraph = AnimatedGraph(x -> 0)
end

function draw(window::PlotWindow)
	set_view(window.renderwindow, window.view)
	draw(window, window.xaxis)
	draw(window, window.yaxis)
	for i = 1:length(window.graphs)
		draw(window.renderwindow, window.graphs[i])
	end
end

function close(window::PlotWindow)
	SFML.close(window.renderwindow)
end

function waitfor(window::PlotWindow)
	Base.wait(window.task)
end

export PlotWindow, add_graph, redraw, check_input, isopen, draw, close, waitfor, follow


text:
using Test
using DataDeps
using ParameterEstimocean
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity

@testset "SyntheticObservations tests" begin
    @info "  Generating synthetic observations with an Oceananigans.Simulation..."
    # Generate synthetic observations
    Nz = 16
    Lz = 128
    Qᵇ = 1e-8
    Qᵘ = -1e-5
    Δt = 20.0
    f₀ = 1e-4
    N² = 1e-5

    stop_time = 6Δt
    save_interval = 2Δt
    N_ensemble = 3

    experiment_name = "convective_adjustment_test"
    data_path = experiment_name * ".jld2"

    # "True" parameters to be estimated by calibration
    convective_κz = 1.0
    convective_νz = 0.9
    background_κz = 1e-4
    background_νz = 1e-5

    grid = RectilinearGrid(size=Nz, z=(-Lz, 0), topology=(Flat, Flat, Bounded))
    closure = ConvectiveAdjustmentVerticalDiffusivity(; convective_κz, background_κz, convective_νz, background_νz)
    coriolis = FPlane(f=f₀)

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))

    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        tracers = :b,
                                        buoyancy = BuoyancyTracer(),
                                        boundary_conditions = (; u=u_bcs, b=b_bcs),
                                        coriolis = coriolis,
                                        closure = closure)

    set!(model, b = (x, y, z) -> N² * z)

    simulation = Simulation(model; Δt, stop_time)

    init_with_parameters(file, model) = file["parameters"] = (; Qᵇ, Qᵘ, Δt, N², tracers=keys(model.tracers))

    simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                          schedule = TimeInterval(save_interval),
                                                          filename = experiment_name,
                                                          with_halos = true,
                                                          overwrite_existing = true,
                                                          init = init_with_parameters)

    run!(simulation)

    #####
    ##### Test synthetic observations construction
    #####

    @testset "SyntheticObservations construction" begin
        @info "    Testing construction of SyntheticObservations from Oceananigans data..."
        data_path = experiment_name * ".jld2"

        # field_names
        b_observations = SyntheticObservations(data_path, field_names=:b)
        ub_observations = SyntheticObservations(data_path, field_names=(:u, :b))
        uvb_observations = SyntheticObservations(data_path, field_names=(:u, :v, :b))

        @test keys(b_observations.field_time_serieses) == tuple(:b)
        @test keys(ub_observations.field_time_serieses) == tuple(:u, :b)
        @test keys(uvb_observations.field_time_serieses) == tuple(:u, :v, :b)

        # transformations and normalizations
        field_names = (:u, :v, :b)
        transformation = ZScore()
        uvb_observations = SyntheticObservations(data_path; field_names, transformation)
        @test all(t.normalization isa ZScore for t in values(uvb_observations.transformation))

        transformation = (u = ZScore(), v = ZScore(), b = ZScore())
        uvb_observations = SyntheticObservations(data_path; field_names, transformation)
        @test all(t.normalization isa ZScore for t in values(uvb_observations.transformation))

        transformation = (u = nothing, v = ZScore(), b = RescaledZScore(0.1))
        uvb_observations = SyntheticObservations(data_path; field_names, transformation)
        @test uvb_observations.transformation[:u].normalization isa Nothing
        @test uvb_observations.transformation[:v].normalization isa ZScore
        @test uvb_observations.transformation[:b].normalization isa RescaledZScore
        @test uvb_observations.transformation[:b].normalization.scale === 0.1

        # Regridding
        coarse_grid = RectilinearGrid(size=Int(Nz/2), z=(-Lz, 0), topology=(Flat, Flat, Bounded))
        fine_grid = RectilinearGrid(size=2Nz, z=(-Lz, 0), topology=(Flat, Flat, Bounded))

        for regrid in [(1, 1, Int(Nz/2)), coarse_grid]
            coarsened_observations = SyntheticObservations(data_path, field_names=(:u, :v, :b); regrid)
            @test size(coarsened_observations.grid) === (1, 1, Int(Nz/2))
        end

        for regrid in [(1, 1, 2Nz), fine_grid]
            refined_observations = SyntheticObservations(data_path, field_names=(:u, :v, :b); regrid)
            @test size(refined_observations.grid) === (1, 1, 2Nz)
        end

        # Test regridding LESbrary observations
        data_path = datadep"two_day_suite_2m/free_convection_instantaneous_statistics.jld2";
        for Nz in (8, 16, 32, 64, 128, 256, 512)
            observations = SyntheticObservations(data_path; field_names=(:u, :v, :b), regrid=(1, 1, Nz))
            @test size(observations.grid) === (1, 1, Nz)

            raw_obs = SyntheticObservations(data_path; field_names=(:u, :v, :b))
            regrid = RectilinearGrid(size=Nz, z=(-raw_obs.grid.Lz, 0), topology=(Flat, Flat, Bounded))
            observations = SyntheticObservations(data_path; field_names=(:u, :v, :b), regrid)
            @test size(observations.grid) === (1, 1, Nz)
        end
    end

    @testset "SyntheticObservations display" begin
        data_path = experiment_name * ".jld2"
        observations = SyntheticObservations(data_path, field_names=(:u, :v, :b))
        @show observations
        @test true
    end
end


text:
can_do_fast_attack(knight_awake) = !knight_awake
can_spy(knight_awake, archer_awake, prisoner_awake) = knight_awake || archer_awake || prisoner_awake
can_signal_prisoner(archer_awake, prisoner_awake) = prisoner_awake && !archer_awake
can_free_prisoner(knight_awake, archer_awake, prisoner_awake, dog_present) = !knight_awake && !archer_awake && prisoner_awake || dog_present && !archer_awake