text:
#
# fix-fitsio.jl --
#
# Some fixes and additional routines which should be part of FITSIO.jl package.
#
# ----------------------------------------------------------------------------

using FITSIO
using FITSIO.Libcfitsio

import FITSIO: libcfitsio, fits_assert_ok, fits_assert_open
import FITSIO.Libcfitsio: fits_get_errstatus

# The exported functions cfitsio_datatype and fits_bitpix deal with conversion
# between CFITSIO type code or BITPIX value and actual Julia data types.
# They can be used as follows (assuming `T` is a Julia data type, while
# `code` and `bitpix` are integers):
#
#     cfitsio_datatype(T) --------> code (e.g., TBYTE, TFLOAT, etc.)
#     cfitsio_datatype(code) -----> T
#
#     fits_bitpix(T) -------------> bitpix (e.g., BYTE_IMG, FLOAT_IMG, etc.)
#     fits_bitpix(bitpix) --------> T
#
export cfitsio_datatype, fits_bitpix

# The following table gives the correspondances between CFITSIO "types",
# the BITPIX keyword and Julia types.
#
#     -------------------------------------------------
#     BITPIX   CFISTIO         Julia     Comments
#     -------------------------------------------------
#              int             Cint
#              long            Clong
#              LONGLONG        Int64     64-bit integer
#     -------------------------------------------------
#        8     BYTE_IMG        UInt8
#       16     SHORT_IMG       Int16
#       32     LONG_IMG        Int32
#       64     LONGLONG_IMG    Int64
#      -32     FLOAT_IMG       Float32
#      -64     DOUBLE_IMG      Float64
#     -------------------------------------------------
#              TBIT
#              TBYTE           Cuchar = UInt8
#              TSBYTE          Cchar = Int8
#              TLOGICAL        Bool
#              TSHORT          Cshort
#              TUSHORT         Cushort
#              TINT            Cint
#              TUINT           Cuint
#              TLONG           Clong
#              TLONGLONG       Int64
#              TULONG          Culong
#              TFLOAT          Cfloat
#              TDOUBLE         Cdouble
#              TCOMPLEX        Complex{Cfloat}
#              TDBLCOMPLEX     Complex{Cdouble}
#     -------------------------------------------------

# Conversion to a C `int`:
cint(x) = convert(Cint, x)
cint(x::Cint) = x

# BITPIX routines and table.
const _BITPIX = Dict{Cint, DataType}()
for (sym, val, T) in ((:BYTE_IMG,        8,       UInt8),
                      (:SHORT_IMG,      16,       Int16),
                      (:LONG_IMG,       32,       Int32),
                      (:LONGLONG_IMG,   64,       Int64),
                      (:FLOAT_IMG,     -32,       Float32),
                      (:DOUBLE_IMG,    -64,       Float64))
    val = cint(val)
    _BITPIX[val] = T
    @eval begin
        fits_bitpix(::Type{$T}) = $val
    end
end
fits_bitpix(code::Integer) = get(_BITPIX, cint(code), Void)

# Data type routines and table.
const _DATATYPE = Dict{Cint, DataType}()
for (sym, val, T) in ((:TBIT       ,   1, Void),
                      (:TBYTE      ,  11, UInt8),
                      (:TSBYTE     ,  12, Int8),
                      (:TLOGICAL   ,  14, Bool),
                      (:TSTRING    ,  16, AbstractString),
                      (:TUSHORT    ,  20, Cushort),          # Uint16
                      (:TSHORT     ,  21, Cshort),           # Int16
                      (:TUINT      ,  30, Cuint),            # Uint32
                      (:TINT       ,  31, Cint),             # Int32
                      (:TULONG     ,  40, Culong),
                      (:TLONG      ,  41, Clong),
                      (:TFLOAT     ,  42, Cfloat),           # Float32
                      (:TLONGLONG  ,  81, Int64),
                      (:TDOUBLE    ,  82, Cdouble),          # Float64
                      (:TCOMPLEX   ,  83, Complex{Cfloat}),  # Complex64
                      (:TDBLCOMPLEX, 163, Complex{Cdouble})) # Complex128
    val = cint(val)
    _DATATYPE[val] = T
    if T == AbstractString
        @eval cfitsio_datatype{S<:AbstractString}(::Type{S}) = $val
    elseif T != Void
        @eval cfitsio_datatype(::Type{$T}) = $val
    end
end
cfitsio_datatype(code::Integer) = get(_DATATYPE, cint(code), Void)


text:
const TEST_FLOAT32 = false
const TEST_DOUBLE64 = false
const TEST_STATIC_SIZED = false

using Manifolds
using ManifoldsBase
using ManifoldsBase: number_of_coordinates

using LinearAlgebra
using Distributions
using DoubleFloats
using ForwardDiff
using Quaternions
using Random
using ReverseDiff
using StaticArrays
using Statistics
using StatsBase
using Test
using LightGraphs
using SimpleWeightedGraphs

using Manifolds.ManifoldTests

ext: jl


text:
const use_return  = haskey(ENV, "USE_RETURN" ) ? parse(Bool, ENV["USE_RETURN"] ) : false
const USE_GPU     = haskey(ENV, "USE_GPU"    ) ? parse(Bool, ENV["USE_GPU"]    ) : false
const do_viz      = haskey(ENV, "DO_VIZ"     ) ? parse(Bool, ENV["DO_VIZ"]     ) : true
const do_save     = haskey(ENV, "DO_SAVE"    ) ? parse(Bool, ENV["DO_SAVE"]    ) : false
const do_save_viz = haskey(ENV, "DO_SAVE_VIZ") ? parse(Bool, ENV["DO_SAVE_VIZ"]) : false
const nx          = haskey(ENV, "NX"         ) ? parse(Int , ENV["NX"]         ) : 512
###
using ParallelStencil
using ParallelStencil.FiniteDifferences1D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 1)
else
    @init_parallel_stencil(Threads, Float64, 1)
end
using Plots, Printf, LinearAlgebra, MAT

@parallel function compute_flux!(qHx, qHx2, H, D, θr_dτ, dx)
    @all(qHx)  = (@all(qHx) * θr_dτ - D * @d(H) / dx) / (1.0 + θr_dτ)
    @all(qHx2) = -D * @d(H) / dx
    return
end

@parallel function compute_update!(H, Hold, qHx, dτ_ρ, dt, dx)
    @inn(H) = (@inn(H) +  dτ_ρ * (@inn(Hold) / dt - @d(qHx) / dx)) / (1.0 + dτ_ρ / dt)
    return
end

@parallel function check_res!(ResH, H, Hold, qHx2, dt, dx)
    @all(ResH) = -(@inn(H) - @inn(Hold)) / dt - @d(qHx2) / dx
    return
end

@views function diffusion_1D_()
    # Physics
    lx     = 10.0       # domain size
    D      = 1          # diffusion coefficient
    ttot   = 1.0        # total simulation time
    dt     = 0.2        # physical time step
    # Numerics
    # nx     = 2*256      # numerical grid resolution
    tol    = 1e-8       # tolerance
    itMax  = 1e5        # max number of iterations
    nout   = 10         # tol check
    CFL    = 0.99       # CFL number
    # Derived numerics
    dx     = lx / nx      # grid size
    Vpdτ   = CFL * dx
    Re     = π + sqrt(π^2 + (lx^2 / D / dt)) # Numerical Reynolds number
    θr_dτ  = lx / Vpdτ / Re
    dτ_ρ   = Vpdτ * lx / D / Re
    xc     = LinRange(dx/2, lx - dx/2, nx)
    # Array allocation
    qHx    = @zeros(nx-1)
    qHx2   = @zeros(nx-1)
    ResH   = @zeros(nx-2)
    dH     = @zeros(nx-2)
    # Initial condition
    H0     = Data.Array(exp.(-(xc .- lx/2).^2 / D))
    Hold   = @ones(nx) .* H0
    H      = @ones(nx) .* H0
    t = 0.0; it = 0; ittot = 0; nt = Int(ceil(ttot/dt))
    # Physical time loop
    while it < nt
        iter = 0; err = 2 * tol
        # Pseudo-transient iteration
        while err > tol && iter < itMax
            @parallel compute_flux!(qHx, qHx2, H, D, θr_dτ, dx)
            @parallel compute_update!(H, Hold, qHx, dτ_ρ, dt, dx)
            iter += 1
            if iter % nout == 0
                @parallel check_res!(ResH, H, Hold, qHx2, dt, dx)
                err = norm(ResH) / sqrt(length(ResH))
            end
        end
        ittot += iter; it += 1; t += dt
        Hold .= H
        if isnan(err) error("NaN") end
    end
    # Analytic solution
    Hana = 1 / sqrt(4 * (ttot + 1 / 4)) * exp.(-(xc .- lx/2).^2 / (4 * D * (ttot + 1 / 4)))
    @printf("Total time = %1.2f, time steps = %d, nx = %d, iterations tot = %d, error vs analytic = %1.2e \n", round(ttot, sigdigits=2), it, nx, ittot, norm(Array(H) - Hana) / sqrt(nx))
    # Visualise
    if do_viz plot(xc, Array(H0), linewidth=3); display(plot!(xc, [Array(H) Array(Hana)], legend=false, framestyle=:box, linewidth=3, xlabel="lx", ylabel="H", title="linear diffusion (nt=$it, iters=$ittot)")) end
    if do_save
        !ispath("../output") && mkdir("../output")
        open("../output/out_diff_1D_lin.txt","a") do io
            println(io, "$(nx) $(ittot) $(nt)")
        end
    end
    if do_save_viz
        !ispath("../out_visu") && mkdir("../out_visu")
        matwrite("../out_visu/diff_1D_lin.mat", Dict("H_1D"=> Array(H), "xc_1D"=> Array(xc)); compress = true)
    end
    return xc, H
end

if use_return
    xc, H = diffusion_1D_();
else
    diffusion_1D = begin diffusion_1D_(); return; end
end


text:
module SolveDSGE

include("types.jl")
include("extra_functions.jl")
include("givens.jl")
include("householder.jl")
include("vec_to_vech.jl")
include("vech_to_vec.jl")
include("derivative.jl")
include("hessian.jl")
include("reorder_complex_schur.jl")
include("check_model.jl")
include("solve_re.jl")

export derivative,
       hessian,
       solve_re,
       check_model_form,
       convert_second_order

export Blanchard_Kahn_Form,
       Blanchard_Kahn_Soln,
       Klein_Form,
       Klein_Soln,
       Structural_Form,
       Structural_Soln,
       Sims_Form,
       Sims_Soln,
       Gomme_Klein_Form,
       Gomme_Klein_Soln,
       Lombardo_Sutherland_Form,
       Lombardo_Sutherland_Soln

end

max_stars_repo_name: dingliumath/SolveDSGE

max_issues_repo_name: dingliumath/SolveDSGE

max_forks_repo_name: dingliumath/SolveDSGE


text:
# Convert from Intervals and numbers to LazySets
to_set(a::AbstractArray) = reduce(×, to_set.(a))
to_set(i::IntervalArithmetic.Interval) = Interval(i)
to_set(s::LazySet) = s
to_set(x) = Singleton([x])


text:
    #
    # Two-photon emission of initially-unpolarized atoms.
    #
    """
    `struct  MultiPhotonDeExcitation.ReducedChannel_2pEmission`
        ... defines a type for a two-photon emission channel for the emission of photons with well-defined
            multipolarities.

        + K              ::AngularJ64             ... Rank K of the channel.
        + omega1         ::Float64                ... omega1.
        + omega2         ::Float64                ... omega2.
        + multipole1     ::EmMultipole            ... Multipole M1.
        + multipole2     ::EmMultipole            ... Multipole M2.
        + gauge          ::EmGauge                ... Gauge for dealing with the (coupled) radiation field.
        + Jsym           ::LevelSymmetry          ... Symmetry of the Green function channel used in the summation.
        + amplitude      ::Complex{Float64}       ... reduced two-photon emission amplitude U^(K, 2gamma, emission) (..)
                                                      associated with the given channel.
    """
    struct  ReducedChannel_2pEmission
        K                ::AngularJ64
        omega1           ::Float64
        omega2           ::Float64
        multipole1       ::EmMultipole
        multipole2       ::EmMultipole
        gauge            ::EmGauge
        Jsym             ::LevelSymmetry
        amplitude        ::Complex{Float64}
    end

    """
    `struct  MultiPhotonDeExcitation.Sharing_2pEmission`
        ... defines a type for a 2pEmission sharing to help characterize energy sharing between the two emitted photons.

        + omega1           ::Float64         ... Energy of the emitted photon 1.
        + omega2           ::Float64         ... Energy of the emitted photon 2.
        + weight           ::Float64         ... Gauss-Lengendre weight of this sharing for energy-integrated quantities.
        + differentialRate ::EmProperty      ... differential rate of this energy sharing.
        + channels         ::Array{MultiPhotonDeExcitation.ReducedChannel_2pEmission,1}
                                             ... List of 2pEmission channels of this sharing.
    """
    struct  Sharing_2pEmission
        omega1             ::Float64
        omega2             ::Float64
        weight             ::Float64
        differentialRate   ::EmProperty
        channels           ::Array{MultiPhotonDeExcitation.ReducedChannel_2pEmission,1}
    end

    """
    `struct  MultiPhotonDeExcitation.Line_2pEmission`
        ... defines a type for a two-photon absorption line by monochromatic light that may include the definition of channels.

        + initialLevel     ::Level                ... initial-(state) level
        + finalLevel       ::Level                ... final-(state) level
        + totalRate        ::EmProperty           ... Total rate for the two-photon transition.
        + sharings         ::Array{MultiPhotonDeExcitation.Sharing_2pEmission,1}
                                                  ... List of MultiPhotonDeExcitation.Sharing_2pEmission's of this line.
    """
    struct  Line_2pEmission
        initialLevel       ::Level
        finalLevel         ::Level
        totalRate          ::EmProperty
        sharings           ::Array{MultiPhotonDeExcitation.Sharing_2pEmission,1}
    end

    # `Base.show(io::IO, line::MultiPhotonDeExcitation.Line_2pEmission)`
    #   ... prepares a proper printout of the variable line::MultiPhotonDeExcitation.Line_2pEmission.
    function Base.show(io::IO, line::MultiPhotonDeExcitation.Line_2pEmission)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "totalRate:         $(line.totalRate)  ")
        println(io, "sharings:          $(line.sharings)  ")
    end

    """
    `MultiPhotonDeExcitation.computeAmplitudesProperties_2pEmission(line::MultiPhotonDeExcitation.Line_2pEmission,
                             grid::Radial.Grid, settings::MultiPhotonDeExcitation.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::MultiPhotonDeExcitation.Line_2pEmission
            is returned for which the amplitudes and properties are now evaluated.
    """
    function  computeAmplitudesProperties_2pEmission(line::MultiPhotonDeExcitation.Line_2pEmission,
                                                     grid::Radial.Grid, settings::MultiPhotonDeExcitation.Settings)
        newSharings = MultiPhotonDeExcitation.Sharing_2pEmission[]
        for sharing in line.sharings
            newChannels = MultiPhotonDeExcitation.ReducedChannel_2pEmission[]
            for ch in sharing.channels
                amplitude = MultiPhotonDeExcitation.computeReducedAmplitudeEmission(ch.K, line.finalLevel, ch.omega2, ch.multipole2,
                                                                ch.Jsym, ch.omega1, ch.multipole1, line.initialLevel, ch.gauge, grid, settings.green)
                push!( newChannels, MultiPhotonDeExcitation.ReducedChannel_2pEmission(ch.K, ch.omega1, ch.omega2, ch.multipole1, ch.multipole2,
                                                                                ch.gauge, ch.Jsym, amplitude) )
            end
            # Calculate the differential rate
            diffRate = EmProperty(-1., -1.)
            push!( newSharings, MultiPhotonDeExcitation.Sharing_2pEmission( sharing.omega1, sharing.omega2, sharing.weight, diffRate, newChannels) )
        end
        # Calculate the totalRate
        totalRate = EmProperty(-1., -1.)
        line = MultiPhotonDeExcitation.Line_2pEmission( line.initialLevel, line.finalLevel, totalRate, newSharings)
        return( line )
    end

    """
    `MultiPhotonDeExcitation.computeLines(process::TwoPhotonEmission, finalMultiplet::Multiplet, initialMultiplet::Multiplet,
                                          grid::Radial.Grid, settings::MultiPhotonDeExcitation.Settings; output=true)`
        ... to compute the multiphoton transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{MultiPhotonDeExcitation.Line_2pEmission,1} is returned.
    """
    function  computeLines(process::TwoPhotonEmission, finalMultiplet::Multiplet, initialMultiplet::Multiplet,
                           grid::Radial.Grid, settings::MultiPhotonDeExcitation.Settings; output=true)
        println("")
        printstyled("MultiPhotonDeExcitation.computeLines(::TwoPhotonEmission): The computation of amplitudes starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        lines = MultiPhotonDeExcitation.determineLines_2pEmission(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    MultiPhotonDeExcitation.displayLines_2pEmission(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = MultiPhotonDeExcitation.Line_2pEmission[]
        for  line in lines
            newLine = MultiPhotonDeExcitation.computeAmplitudesProperties_2pEmission(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        MultiPhotonDeExcitation.displayTotalRates_2pEmission(stdout, lines, settings)
        MultiPhotonDeExcitation.displayDifferentialRates_2pEmission(stdout, lines, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   MultiPhotonDeExcitation.displayTotalRates_2pEmission(iostream, lines, settings)
                           MultiPhotonDeExcitation.displayDifferentialRates_2pEmission(iostream, lines, settings)     end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `MultiPhotonDeExcitation.computeEnergyDiffCs(omega1::Float64, omega2::Float64, line::MultiPhotonDeExcitation.Line_2pEmission,
                                                 settings::MultiPhotonDeExcitation.Settings)`
        ... to compute the energy differential cross section for the given line and the energy sharing omega1. A dcs::Float64 is returned.
    """
    function computeEnergyDiffCs(omega1::Float64, omega2::Float64, line::MultiPhotonDeExcitation.Line_2pEmission, gauge::EmGauge,
                                 settings::MultiPhotonDeExcitation.Settings)
        dcs = 0.
        for  mp1 in settings.multipoles
            for  mp2 in settings.multipoles
                symi        = LevelSymmetry(line.initialLevel.J, line.initialLevel.parity);    symf = LevelSymmetry(line.finalLevel.J, line.finalLevel.parity)
                symmetries  = AngularMomentum.allowedTotalSymmetries(symf, mp2, mp1, symi)
                Klist       = oplus(line.finalLevel.J, line.initialLevel.J)
                for Jsym in symmetries
                    wk  = ComplexF64(0.);   wm = ComplexF64(0.)
                    for  K in Klist
                        wk =  wk + MultiPhotonDeExcitation.getReducedAmplitudeEmission(K, line.finalLevel, omega2, mp2, Jsym, omega1, mp1,
                                                                                line.initialLevel, gauge, line.channels) +
                               AngularMomentum.phaseFactor([K, +1, line.finalLevel.J, +1, line.initialLevel.J]) *
                                   MultiPhotonDeExcitation.getReducedAmplitudeEmission(K, line.finalLevel, omega1, mp1, Jsym, omega2, mp2,
                                                                                line.initialLevel, gauge, line.channels)
                    end
                    dcs = dcs + abs( wk )^2
                end
            end
        end
        dcs = dcs * 2pi * Defaults.getDefaults("alpha")^2 / (Basics.twice(line.initialLevel.J) + 1) * omega1 * omega2

        return( dcs )
    end

    """
    `MultiPhotonDeExcitation.computeReducedAmplitudeEmission(K::AngularJ64, finalLevel::Level, omega2::Float64, multipole2::EmMultipole, Jsym::LevelSymmetry,
                                                                                omega1::Float64, multipole1::EmMultipole, initialLevel::Level,
                                                                               gauge::EmGauge, grid::Radial.Grid, greenChannels::Array{AtomicState.GreenChannel,1})`
        ... to compute the reduced amplitude U^{K, 2gamma emission} (K, Jf, omega2, multipole2, Jsym, omega1, multipole1, Ji) by means of the
            given Green function channels. An amplitude::Complex{Float64} is returned.
    """
    function computeReducedAmplitudeEmission(K::AngularJ64, finalLevel::Level, omega2::Float64, multipole2::EmMultipole, Jsym::LevelSymmetry,
                                                                               omega1::Float64, multipole1::EmMultipole, initialLevel::Level,
                                            gauge::EmGauge, grid::Radial.Grid, greenChannels::Array{AtomicState.GreenChannel,1})
        U = Complex(0.);    found = false
        for channel in greenChannels
            if  Jsym == channel.symmetry
                found = true
                for (i, nuLevel) in enumerate(channel.gMultiplet.levels)
                    U = U + PhotoEmission.amplitude("emission", multipole2, gauge, omega2, finalLevel, nuLevel, grid) *
                            PhotoEmission.amplitude("emission", multipole1, gauge, omega1, nuLevel, initialLevel, grid) /
                            (initialLevel.energy + omega1 - nuLevel.energy)
                end
            end
        end

        if    found
              U = U * AngularMomentum.Wigner_6j(initialLevel.J, finalLevel.J, K, AngularJ64(multipole2.L), AngularJ64(multipole1.L), Jsym.J)
              ##x println(" ")
        else  println("No Green function cannel found for amplitude U^{K, 2gamma emission} (K, Jf, omega2, multipole2, Jsym = $Jsym, omega1, multipole1, Ji) ")
        end

        return( U )
    end

    """
    `MultiPhotonDeExcitation.getReducedAmplitudeEmission(K::AngularJ64, finalLevel::Level, omega2::Float64, multipole2::EmMultipole, Jsym::LevelSymmetry,
                                                                                omega1::Float64, multipole1::EmMultipole, initialLevel::Level,
                                                                           gauge::EmGauge, channels::Array{MultiPhotonDeExcitation.ReducedChannel_2pEmission,1})`
        ... to get/return the reduced amplitude U^{K, 2gamma emission} (K, Jf, omega2, multipole2, Jsym, omega1, multipole1, Ji) from the calculated list
            of channels. An amplitude::Complex{Float64} is returned.
    """
    function getReducedAmplitudeEmission(K::AngularJ64, finalLevel::Level, omega2::Float64, multipole2::EmMultipole, Jsym::LevelSymmetry,
                                                                           omega1::Float64, multipole1::EmMultipole, initialLevel::Level,
                                                                           gauge::EmGauge, channels::Array{MultiPhotonDeExcitation.ReducedChannel_2pEmission,1})
        U = Complex(0.);    found = false
        for channel in channels
            if  K == channel.K  &&  omega1 == channel.omega1  &&  omega2 == channel.omega2  &&  Jsym == channel.Jsym  &&
                multipole1 == channel.multipole1  &&  multipole2 == channel.multipole2      &&
                (gauge == channel.gauge  ||  EmGauge("Magnetic")  == channel.gauge)
                U = channel.amplitude;    found = true
            end
        end

        if    found
              println("U^{K, 2gamma emission} (..) = $U found. ")
        else  println("No U^{K, 2gamma emission} (K, Jf, omega2, multipole2, Jsym = $Jsym, omega1, multipole1, Ji) amplitude found.")
        end

        return( U )
    end

    """
    `MultiPhotonDeExcitation.determineLines_2pEmission(finalMultiplet::Multiplet, initialMultiplet::Multiplet,
                                                       settings::MultiPhotonDeExcitation.Settings)`
        ... to determine a list of MultiPhotonDeExcitation.Line_2pEmission's for transitions between the levels from the given
            initial- and final-state multiplets and by taking into account the particular selections and settings for this computation;
            an Array{MultiPhotonDeExcitation.Line_2pEmission,1} is returned. Apart from the level specification, all physical
            properties are set to zero during this initialization process.
    """
    function  determineLines_2pEmission(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::MultiPhotonDeExcitation.Settings)
        lines = MultiPhotonDeExcitation.Line_2pEmission[]
        for  iLevel  in  initialMultiplet.levels
            for  fLevel  in  finalMultiplet.levels
                if  Basics.selectLevelPair(iLevel, fLevel, settings.lineSelection)
                    energy   = iLevel.energy - fLevel.energy
                    sharings = MultiPhotonDeExcitation.determineSharingsAndChannels(fLevel, iLevel, energy, settings)
                    push!( lines, MultiPhotonDeExcitation.Line_2pEmission(iLevel, fLevel, EmProperty(0., 0.,), sharings) )
                end
            end
        end
        return( lines )
    end

    """
    `MultiPhotonDeExcitation.determineSharingsAndChannels(finalLevel::Level, initialLevel::Level, energy::Float64, settings::MultiPhotonDeExcitation.Settings)`
        ... to determine a list of MultiPhotonDeExcitation 2pEmission Sharing's and Channel's for a transitions from the initial to
            final level and by taking into account the particular settings of for this computation;
            an Array{MultiPhotonDeExcitation.Sharing_2pEmission,1} is returned.
    """
    function determineSharingsAndChannels(finalLevel::Level, initialLevel::Level, energy::Float64, settings::MultiPhotonDeExcitation.Settings)
        sharings = MultiPhotonDeExcitation.Sharing_2pEmission[];    eSharings = Basics.determineEnergySharings(energy, settings.NoEnergySharings)
        for  es in eSharings
            omega1    = es[1];    omega2 = es[2];    weight = es[3]
            channels  = MultiPhotonDeExcitation.ReducedChannel_2pEmission[];
            symi      = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
            for  mp1 in settings.multipoles
                for  mp2 in settings.multipoles
                    symmetries  = AngularMomentum.allowedTotalSymmetries(symf, mp2, mp1, symi)
                    Klist       = oplus(symf.J, symi.J)
                    for  symx in symmetries
                        for  gauge in settings.gauges
                            # Include further restrictions if appropriate
                            if     string(mp1)[1] == 'E' || string(mp2)[1] == 'E'  &&   gauge == Basics.UseCoulomb
                                for K in Klist  push!(channels, ReducedChannel_2pEmission(K, omega1, omega2, mp1, mp2, Basics.Coulomb, symx, 0.) )    end
                            elseif string(mp1)[1] == 'E' || string(mp2)[1] == 'E'  &&   gauge == Basics.UseBabushkin
                                for K in Klist  push!(channels, ReducedChannel_2pEmission(K, omega1, omega2, mp1, mp2, Basics.Babushkin, symx, 0.) )  end
                            elseif string(mp1)[1] == 'M' && string(mp2)[1] == 'M'
                                for K in Klist  push!(channels, ReducedChannel_2pEmission(K, omega1, omega2, mp1, mp2, Basics.Magnetic, symx, 0.) )   end
                            end
                        end
                    end
                end
            end
            push!(sharings, MultiPhotonDeExcitation.Sharing_2pEmission(omega1, omega2, weight, EmProperty(0., 0.), channels) )
        end
        return( sharings )
    end

    """
    `MultiPhotonDeExcitation.displayDifferentialRates_2pEmission(stream::IO, lines::Array{MultiPhotonDeExcitation.Line_2pEmission,1},
                                                                 settings::MultiPhotonDeExcitation.Settings)`
        ... to display all differential rates, etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayDifferentialRates_2pEmission(stream::IO, lines::Array{MultiPhotonDeExcitation.Line_2pEmission,1},
                                                  settings::MultiPhotonDeExcitation.Settings)
        #
        # First, print lines and sharings
        nx = 130
        println(stream, " ")
        println(stream, "  Energy-differential rates of selected two-photon emission lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(38, "Energies (all in " * TableStrings.inUnits("energy") * ")"; na=4);
        sb = sb * TableStrings.flushleft(38, "  i -- f          omega1        omega2"; na=4)
        sa = sa * TableStrings.center(14, "Weight"; na=0);                            sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(34, "Cou -- diff. rate -- Bab"; na=3)
        sb = sb * TableStrings.center(34, TableStrings.inUnits("rate") * "        " *
                                          TableStrings.inUnits("rate"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";      isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy)) * "    "
            #
            for  (is, sharing)  in  enumerate(line.sharings)
                if  is == 1     sb = sa     else    sb = TableStrings.hBlank( length(sa) )    end
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.omega1))                   * "    "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.omega2))                   * "    "
                sb = sb * @sprintf("%.4e",                                              sharing.weight)                    * "       "
                sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", sharing.differentialRate.Coulomb))   * "   "
                sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", sharing.differentialRate.Babushkin)) * "   "
                println(stream,  sb )
            end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `MultiPhotonDeExcitation.displayLines_2pEmission(lines::Array{MultiPhotonDeExcitation.Line_2pEmission,1})`
        ... to display a list of lines, sharings & (reduced) channels that have been selected due to the prior settings.
            A neat table of all selected transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines_2pEmission(lines::Array{MultiPhotonDeExcitation.Line_2pEmission,1})
        nx = 157
        println(" ")
        println("  Selected two-photon emission lines with given photon splitting:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(34, "Energies (all in " * TableStrings.inUnits("energy") *")"; na=5);
        sb = sb * TableStrings.flushleft(34, "  i -- f      omega1      omega2  "; na=5)
        sa = sa * TableStrings.flushleft(77, "List of multipoles, gauges & intermediate level symmetries"; na=4)
        sb = sb * TableStrings.flushleft(77, "(K-rank, multipole_1, Jsym, multipole_2, gauge)           "; na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", energy)) * "  "
            #
            for  sharing  in  line.sharings
                sb =      @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.omega1))   * "  "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.omega2))   * "     "
                mpGaugeList = Tuple{AngularJ64, Basics.EmMultipole, LevelSymmetry, Basics.EmMultipole, Basics.EmGauge}[]
                for  channel in  sharing.channels
                    push!( mpGaugeList, (channel.K, channel.multipole1, channel.Jsym, channel.multipole2, channel.gauge) )
                end
                wa = TableStrings.twoPhotonGaugeTupels(75, mpGaugeList)
                sc = sa * sb * wa[1];    println( sc )
                for  i = 2:length(wa)
                    sc = TableStrings.hBlank( length(sa*sb) ) * wa[i];    println( sc )
                end
            end
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `MultiPhotonDeExcitation.displayTotalRates_2pEmission(stream::IO, lines::Array{MultiPhotonDeExcitation.Line_2pEmission,1},
                                                          settings::MultiPhotonDeExcitation.Settings)`
        ... to display all total rates, etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayTotalRates_2pEmission(stream::IO, lines::Array{MultiPhotonDeExcitation.Line_2pEmission,1}, settings::MultiPhotonDeExcitation.Settings)
        #
        # First, print lines and sharings
        nx = 88
        println(stream, " ")
        println(stream, "  Total rates of selected two-photon emission lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=2);                         sb = sb * TableStrings.hBlank(21)
        sa = sa * TableStrings.center(12, "i--Energy--f"; na=3)
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(30, "Cou --   total rate   -- Bab"; na=4)
        sb = sb * TableStrings.center(30, TableStrings.inUnits("rate") * "          " *
                                          TableStrings.inUnits("rate"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";      isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", energy))                 * "      "
            #
            sb = sa * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", line.totalRate.Coulomb))   * "     "
            sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", line.totalRate.Babushkin)) * "   "
            println(stream,  sb )
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end


text:
using Turing
using Distributions
using Random
using LinearAlgebra
using StatsPlots
using MCMCChains
using DataFrames
using Printf

# geberate moc data from logit model
n = 1_000
num_var = 2
Random.seed!(99)
X = rand(Uniform(-sqrt(3.0), 2.0*sqrt(3.0)), n, num_var) # case-specific variables Z_i Matrix
X = [ones(n) X]
β = Vector([0.0, 0.5, -0.5]) # true coefficients
# ポアソン分布の平均をモデル化
λ = exp.(X*β)
y = rand.(Poisson.(λ)) # observable choice 0 or 1

# set prior, β ~ MvNormal(β0, A0), no σ²
n, k = size(X)
β0 = zeros(k)
A0 = 0.01 * I

@model function poisson_regression_model(X, y)
    # Set the priors only for β.
    β ~ MvNormal(β0, inv(A0))

    # Write the likelihood, can't vectorize?
    n = size(X, 1)
    for i = 1:n
        y[i] ~ Poisson(exp(X[i, :]' * β))
    end
end

# Draw random number from posterior.
model = poisson_regression_model(X, y)

n_draws = 5_000
# n_chains = 4
n_chains = 1
n_tune = 1_000
chn = sample(model, NUTS(), MCMCThreads(), n_draws, n_chains; discard_adapt=false)
plot(chn)
# chn = sample(model, NUTS(), MCMCThreads(), n_draws, n_chains)[n_tune:n_draws, :, :]
β_chn = chn[n_tune+1:n_draws, :, :] # delete initial and exclude odd (hamiltonian bluh bluh) chain
β_chn = group(β_chn, :β)

summarize(β_chn)
quantile(β_chn)
plot(β_chn)

# Marginal effect ... interpret β just like as ln(y) = βx model (equivalent with λ = exp(x))


text:
function Integrate(f, x0, tspan, N, M, tol, imax)

    # Get length of state vector
    L = length(x0)

    # Transformation of independant variable
    ω₁ = 0.5*(tspan[2] + tspan[1])
    ω₂ = 0.5*(tspan[2] - tspan[1])

    g(τ,x) =   ω₂*f(ω₁ + ω₂*τ,x)

    # Compute Chebyshev-Gauss-Lobatto Nodes
    τs = SizedVector{M + 1}(zeros(M + 1))
    @inbounds for j in 0:M
        τs[j + 1] = -cos(j*π/M)
    end

    # Compute constant matricies
    Cx  = SizedMatrix{M + 1, N + 1}(zeros(M + 1, N + 1))
    ChebyshevInterpMat!(Cx, τs)
    CI1 = SizedMatrix{N + 1, N}(zeros(N + 1, N))
    initCI1!(CI1)
    Cf  = SizedMatrix{N, M + 1}(zeros(N, M + 1))
    initCf!(Cf, τs)
    C   = SizedMatrix{N + 1, M + 1}(zeros(N + 1, M + 1))
    mul!(C, CI1, Cf)

    # Initialize required matricies
    Xold  = SizedMatrix{M + 1, L}(zeros(M + 1, L))
    XoldT = Xold'
    Xnew  = SizedMatrix{M + 1, L}(zeros(M + 1, L))
    X0    = SizedMatrix{N + 1, L}(zeros(N + 1, L))
    X0[1,:] .= x0
    G     = SizedMatrix{M + 1, L}(zeros(M + 1, L))
    β     = SizedMatrix{N + 1, L}(zeros(N + 1, L))

    # Set Initial Guess for Xold (Could provide as input in future)
    Xold[1,:] .= x0

    # Begin integration loop
    stop = false
    eold = Inf
    enew = Inf
    iter = 0
    while !stop
        iter += 1
        println(iter)

        # Compute forcing function
        @inbounds for i in 1:M+1
            L == 1 ? G[i,1]  = g(τs[i], Xold[i]) :
                     G[i,:] .= g(τs[i], @view(XoldT[:, i]))
        end

        # Update Coefficients
        #β = X0 + C*G
        mul!(β,C,G)
        β .= β + X0

        # Update State
        #Xnew = Cx*β
        mul!(Xnew, Cx, β)

        # Compute error
        enew = 0
        for col in 1:L
            for row in 1:M+1
                temp = Xnew[row,col] - Xold[row,col]
                if temp > enew
                    enew = temp
                end
            end
        end

        if (enew <= tol && eold <= tol) || iter >= imax
            stop = true
        else
            Xold .= Xnew
            eold  = enew
        end
    end

    return (ω₁ .+ ω₂*τs, Xnew)
end

function initCI1!(CI1)
    # Get N
    dims = size(CI1)
    N = dims[2]

    # Initialize R and S
    R = Diagonal(SizedVector{N + 1}(zeros(N + 1)))
    S = SizedMatrix{N + 1, N}(zeros(N + 1, N))

    # Fill Matricies
    R[1,1] = 1
    R[2,2] = 0.5
    S[1,1] = 1
    S[1,2] = -0.25
    S[2,1] = 2
    S[2,3] = -1
    @inbounds for i in 3:N + 1
        R[i,i] = 1/(2*(i - 1))
        S[i,i - 1] = 1
        if i < N
            S[i,i + 1] = -1
        end
        if i < N + 1
            S[1,i] = (1/(2*(i - 2)) - 1/(2*i))*(-1)^i
        end
    end

    # Fill CI1
    #CI1 .= R*S
    lmul!(R,S)
    CI1 .= S
end

function initCf!(Cf, τs)
    # Get N and M
    dims = size(Cf)
    N = dims[1]
    M = dims[2] - 1

    # Set p variable depending on number of nodes and check Cf size for validity
    p = 1
    if M > N
        p = 2
    elseif N > M
        throw(ArgumentError("Improper size of Cf matrix!"))
    end

    # Initialize temporary matricies
    V  = Diagonal(SizedVector{N}(zeros(N)))
    W  = Diagonal(SizedVector{M + 1}(zeros(M + 1)))
    Tm = SizedMatrix{N, M + 1}(zeros(N, M + 1))

    # Fill T
    T!(Tm, τs)

    # Fill V
    @inbounds for i in 1:N
        # Get numerator
        i == 1 ? num = 1 : i == N ? num = p : num = 2

        # Set value
        V[i, i] = num/M
    end

    # Fill W
    @inbounds for i in 1:M + 1
        # Get val
        i == 1 ? val = 0.5 : i == M + 1 ? val = 0.5 : val = 1

        # Set value
        W[i, i] = val
    end

    # Compute Cf
    #Cf .= V*Tm*W
    rmul!(Tm,W)
    lmul!(V,Tm)
    Cf .= Tm
end


text:
module CmdTools

function execute(cmd::Cmd)
    out = Pipe()
    err = Pipe()

    process = run(pipeline(ignorestatus(cmd), stdout=out, stderr=err))
    close(out.in)
    close(err.in)
    stdout = @async String(read(out))
    stderr = @async String(read(err))
    (
      out = fetch(stdout),
      err = fetch(stderr),
      code = process.exitcode
    )
end

end # module CmdTools


text:
module DebugArtifacts

export debug_artifact

# Adapt this script to debug why artifact installation might be failing

using Pkg, Pkg.Artifacts, Pkg.BinaryPlatforms, Pkg.PlatformEngines, InteractiveUtils

function get_meta(artifact_name, artifacts_toml, platform = platform_key_abi())
    artifact_dict = Pkg.Artifacts.load_artifacts_toml(artifacts_toml)
    meta = artifact_meta(artifact_name, artifact_dict, artifacts_toml; platform=platform)

    treehash = meta["git-tree-sha1"]
    tarball_url = first(meta["download"])["url"]
    tarball_hash = first(meta["download"])["sha256"]
    return treehash, tarball_url, tarball_hash
end

function print_download_artifact_cmd(treehash, tarball_url, tarball_hash)
    println()
    println("using Pkg.Artifacts")
    println("""
download_artifact(Base.SHA1($(repr(treehash))),
                  $(repr(tarball_url)),
                  $(repr(tarball_hash));
                  verbose = true)"""
            )
end

function print_download_artifact_cmd(artifact_name::String, platform = platform_key_abi())
    artifacts_toml_url = "https://raw.githubusercontent.com/JuliaBinaryWrappers/$(artifact_name)_jll.jl/master/Artifacts.toml"
    mktempdir() do tmp_dir
        artifacts_toml = joinpath(tmp_dir, "Artifacts.toml")
        Pkg.PlatformEngines.download(artifacts_toml_url, artifacts_toml; verbose = false)
        treehash, tarball_url, tarball_hash = get_meta(artifact_name, artifacts_toml, platform)
        @info "Run the following commands to install the artifact `$(artifact_name)_jll`:"
        print_download_artifact_cmd(treehash, tarball_url, tarball_hash)
    end
end

# Debug code contributed by @staticfloat
function debug_artifact(artifact_name::String, platform = platform_key_abi())
    @info "Platform: $(platform)"
    InteractiveUtils.versioninfo()
    println()

    # Initialize Pkg code
    probe_platform_engines!(; verbose=true)

    # Change these to whatever you need them to be, to debug your artifacts code
    artifacts_toml_url = "https://raw.githubusercontent.com/JuliaBinaryWrappers/$(artifact_name)_jll.jl/master/Artifacts.toml"

    mktempdir() do tmp_dir
        # First, download Artifacts.toml file to temporary directory
        artifacts_toml = joinpath(tmp_dir, "Artifacts.toml")
        @info("Downloading Artifacts.toml to $(artifacts_toml)...")
        Pkg.PlatformEngines.download(artifacts_toml_url, artifacts_toml; verbose=true)

        # Extract artifact metadata for the current platform (you can override this)
        @info("Extracting artifact info for platform $(triplet(platform))...")
        treehash, tarball_url, tarball_hash = get_meta(artifact_name, artifacts_toml, platform)

        @info("Found meta object with git-tree-sha1 $(treehash), attempting download...")
        destdir = joinpath(tmp_dir, "unpacked")
        download_verify_unpack(tarball_url, tarball_hash, destdir; verbose=true)

        if !Sys.iswindows()
            @info("Double-checking git-tree-sha (this is skipped on Windows)")
            calc_treehash = bytes2hex(Pkg.GitTools.tree_hash(destdir))
            if calc_treehash != treehash
                @warn("Calculated treehash ($(calc_treehash)) != promised treehash ($(treehash))!")
            end
        end
    end

end

end # module


text:
using PyPlot

function linescatter()
    x = [1, 2, 3, 4]
    y = [10, 15, 13, 17]
    plot(x,y)
end
linescatter()

function multiple_scatter_traces()
    x=1:4
    y1 = [10, 15, 13, 17]
    y2 = [16, 5, 11, 9]
    y3 = [12, 9, 15, 12]
    y4 = [5, 10, 8, 12]
    plot(y1, label="marker", "o")
    hold(true)
    plot(y2, label="line+marker", "-o")
    plot(y3, label="line", "-")
    plot(y4, label="dash", "--")
end
multiple_scatter_traces()

function data_labels()
    x = 1:5
    y1 = [1, 6, 3, 6, 1]
    y2 = [4, 1, 7, 1, 4]
    y1_labels = ["A-1", "A-2", "A-3", "A-4", "A-5"]
    y2_labels = ["B-a", "B-b", "B-c", "B-d", "B-e"]

    f, ax = subplots(1, 1)
    title("Data Labels on the Plot")
    plot(x, y1, label="Team A")
    plot(x, y2, label="Team B")
#     in python ax.legend(loc="upper left", shadow=true)
    ax[:legend](loc="upper left", shadow=true)

    for idx= 1:length(x)
        ax[:annotate](s=y1_labels[idx], xy = (x[idx], y1[idx]))
        ax[:annotate](s=y2_labels[idx], xy = (x[idx], y2[idx]))
        idx = idx + 1
    end
end
data_labels()

function area1()
    fig, ax = subplots()
    y1=[0, 2, 3, 5] #, fill="tozeroy")
    y2=[3, 5, 1, 7] #, fill="tonexty")
    plot(y1, label="lines")
    plot(y2, label="lines")
    ax[:legend](loc="upper center", shadow=true)
end
area1()

function matrix_subplots()

end

matrix_subplots()

function advanced_layouts()

end
advanced_layouts()

function grouped_bar_example()

end

grouped_bar_example()

function stacked_bar_example()

end
stacked_bar_example()

function two_hists()

end

two_hists()

function box_plot()

end
box_plot()


text:
# const getf_scheduler = Libdl.dlsym(polo_lib, :getf_scheduler)
# const getx_scheduler = Libdl.dlsym(polo_lib, :getx_scheduler)

struct Scheduler <: ParameterServer
    fval::Float64
    x::Vector{Float64}
    psopts::ParameterServerOptions

    function (::Type{Scheduler})(; kw...)
        scheduler =  new(0., Vector{Float64}(), ParameterServerOptions(; kw...))
    end
end

function initialize!(scheduler::Scheduler, x₀::AbstractVector)
    resize!(scheduler.x,length(x₀))
    scheduler.x .= x₀
end
getx(scheduler::Scheduler) = scheduler.x

function (proxgrad::ProxGradient{Scheduler})(x₀::AbstractVector,termination::AbstractTermination,logger::AbstractLogger)
    initialize!(proxgrad.execution, x₀)
    xbegin = pointer(x₀, 1)
    xend = pointer(x₀, length(x₀) + 1)
    termination_c = @cfunction(POLO.termination_wrapper, Cint,
                               (Cint, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Nothing}))
    log_c = @cfunction(POLO.log_wrapper, Nothing,
                       (Cint, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Nothing}))
    return ccall(POLO.run_scheduler, Nothing,
                 (Ptr{Nothing},Ptr{Cdouble}, Ptr{Cdouble},
                  Ptr{Nothing}, Any,
                  Ptr{Nothing}, Any),
                 proxgrad,
                 xbegin, xend,
                 termination_c, termination,
                 log_c, logger)
end
function (proxgrad::ProxGradient{Scheduler})(x₀::AbstractVector)
    termination = POLO.Utility.MaxIteration(100)
    return proxgrad(x₀,termination,POLO.Utility.ProgressLogger.Value(termination))
end

paramserver_handle(::Scheduler) = POLO.proxgradient_scheduler
execution_handle(::Scheduler) = POLO.run_scheduler
delete_handle(::Scheduler) = POLO.delete_proxgradient_scheduler
# getf_handle(::Scheduler) = getf_scheduler
# getx_handle(::Scheduler) = getx_scheduler


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability, Plots

include(@relpath "lorenz_model.jl")
include(@relpath "lorenz_specifications.jl")

S = lorenz_model()
X0, options = lorenz_specification()
options = Options(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# reachability algorithm
algorithm_reach = TMJets(:abs_tol=>1e-27, :orderT=>10, :orderQ=>2,
                         :max_steps=>50_000)

# compute flowpipe
options[:mode] = "reach"
solution = solve(problem, options; op=algorithm_reach)

# project flowpipe to x₁ and x₃
solution.options[:plot_vars] = [1, 3]
solution_proj = project(solution)

# plot projection
plot(solution_proj)
savefig("lorenz")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
"""
Struct for Barrier Down and In Option

		barOption=BarrierOptionDownIn(T::num1,K::num2,barrier::num3,isCall::Bool=true) where {num1 <: Number,num2 <: Number,num3 <: Number}

Where:\n
		T	=	Time to maturity of the Option.
		K	=	Strike Price of the Option.
		barrier	=	Down Barrier of the Option.
		isCall  = true for CALL, false for PUT.
"""
struct BarrierOptionDownIn{num1 <: Number, num2 <: Number, num3 <: Number, numtype <: Number} <: BarrierPayoff{numtype}
    T::num1
    K::num2
    barrier::num3
    isCall::Bool
    function BarrierOptionDownIn(T::num1, K::num2, barrier::num3, isCall::Bool = true) where {num1 <: Number, num2 <: Number, num3 <: Number}
        if T <= 0.0
            error("Time to Maturity must be positive")
        elseif K <= 0.0
            error("Strike Price must be positive")
        elseif barrier <= 0.0
            error("Barrier must be positive")
        else
            zero_typed = zero(num1) + zero(num2) + zero(num3)
            return new{num1, num2, num3, typeof(zero_typed)}(T, K, barrier, isCall)
        end
    end
end

export BarrierOptionDownIn;

function payout(S::abstractArray, barrierPayoff::BarrierOptionDownIn) where {abstractArray <: AbstractArray{num_}} where {num_ <: Number}
    iscall = barrierPayoff.isCall ? 1 : -1
    zero_typed = zero(eltype(S)) * barrierPayoff.K * barrierPayoff.barrier
    return max(iscall * (S[end] - barrierPayoff.K), zero_typed) * (findfirst(x -> x < barrierPayoff.barrier, S) != nothing)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test

@testset begin
    @test false
    return
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using Test: @test, @testset

import Random

using JuMP: JuMP, @NLconstraint, @objective, @variable, @NLexpression
using SparseArrays: spzeros
using PartiallyObservedInverseGames.ForwardOptimalControl: solve_optimal_control
using PartiallyObservedInverseGames.InverseOptimalControl: solve_inverse_optimal_control
using PartiallyObservedInverseGames.TrajectoryVisualization: visualize_trajectory
using PartiallyObservedInverseGames.CostUtils: symbol

unique!(push!(LOAD_PATH, joinpath(@__DIR__, "utils")))
import TestUtils
import TestDynamics

#========================================== Cost Library ===========================================#

function register_shared_forward_cost_expressions!(opt_model, x, u; prox_min_regularization = 0.1)
    T = size(x, 2)
    @NLexpression(
        opt_model,
        regularized_sq_dist[t = 1:T],
        (x[1, t] - obstacle[1])^2 + (x[2, t] - obstacle[2])^2 + prox_min_regularization
    )
end

function add_forward_objective!(opt_model, x, u; weights)
    T = size(x, 2)
    register_shared_forward_cost_expressions!(opt_model, x, u)

    # Avoid a point. Assumes x = [px, py, ...]. Functional form is -log(|(x, y) - p|^2).
    @variable(opt_model, prox_cost[t = 1:T])
    @NLconstraint(opt_model, [t = 1:T], prox_cost[t] == -log(opt_model[:regularized_sq_dist][t]))

    g̃ = (;
        state_goal = sum(x[1:2, T_activate_goalcost:T] .^ 2),
        state_velocity = sum(x[3, :] .^ 2),
        state_proximity = sum(prox_cost),
        control_Δv = sum(u[1, :] .^ 2),
        control_Δθ = sum(u[2, :] .^ 2),
        control = sum(u .^ 2),
    )

    @objective(opt_model, Min, sum(weights[k] * g̃[symbol(k)] for k in keys(weights)))
end

function add_forward_objective_gradients!(opt_model, x, u; weights)
    n_states, T = size(x)
    n_controls = size(u, 1)
    register_shared_forward_cost_expressions!(opt_model, x, u)
    @variable(opt_model, dproxdx[1:T])
    @NLconstraint(
        opt_model,
        [t = 1:T],
        dproxdx[t] == -2 * (x[1, t] - obstacle[1]) / opt_model[:regularized_sq_dist][t]
    )
    @variable(opt_model, dproxdy[1:T])
    @NLconstraint(
        opt_model,
        [t = 1:T],
        dproxdy[t] == -2 * (x[2, t] - obstacle[2]) / opt_model[:regularized_sq_dist][t]
    )

    dJ̃dx = (;
        state_goal = 2 * [
            zeros(2, T_activate_goalcost - 1) x[1:2, T_activate_goalcost:T]
            zeros(n_states - 2, T)
        ],
        state_velocity = 2 * [spzeros(T, 2) x[3, :] spzeros(T)]',
        state_proximity = [dproxdx dproxdy spzeros(T, n_states - 2)]',
        control_Δv = spzeros(n_states, T),
        control_Δθ = spzeros(n_states, T),
        control = spzeros(n_states, T),
    )
    dJdx = sum(weights[k] * dJ̃dx[symbol(k)] for k in keys(weights))

    dJ̃du = (;
        state_goal = spzeros(n_controls, T),
        state_velocity = spzeros(n_controls, T),
        state_proximity = spzeros(n_controls, T),
        control_Δv = 2 * [u[1, :] zeros(T)]',
        control_Δθ = 2 * [spzeros(T) u[2, :]]',
        control = 2 * u,
    )
    dJdu = sum(weights[k] * dJ̃du[symbol(k)] for k in keys(weights))

    (; dx = dJdx, du = dJdu)
end

control_system = TestDynamics.Unicycle(0.1)

x0 = [-1, 1, 1, 0]
T = 100
obstacle = (-0.5, 0.25) # Point to avoid.
T_activate_goalcost = T
cost_model = (
    weights = (;
        state_goal = 100,
        state_velocity = 1,
        state_proximity = 1,
        control_Δv = 10,
        control_Δθ = 1,
    ),
    add_objective! = add_forward_objective!,
    add_objective_gradients! = add_forward_objective_gradients!,
)

#====================================== forward optimal control ====================================#

@testset "Forward Solution" begin
    forward_converged, forward_solution, forward_model =
        solve_optimal_control(control_system, cost_model, x0, T)
    @test forward_converged
    visualize_trajectory(control_system, forward_solution.x)
    global forward_solution
end

#===================================== Inverse Optimal Control =====================================#

observation_model = (; σ = 0.0, expected_observation = identity)

y = TestUtils.noisy_observation(observation_model, forward_solution.x)

inverse_converged, inverse_solution, inverse_model = solve_inverse_optimal_control(
    y;
    init = (; forward_solution.u),
    control_system,
    cost_model,
    observation_model,
)

#============================================== Tests ==============================================#

@testset "Inverse Solution" begin
    @test inverse_converged
    TestUtils.test_inverse_solution(inverse_solution.weights, cost_model.weights)
    TestUtils.test_data_fidelity(inverse_model, observation_model, forward_solution.x, y)
end


text:
# ---
# title: 1559. Detect Cycles in 2D Grid
# id: problem1559
# author: Indigo
# date: 2021-01-22
# difficulty: Hard
# categories: Depth-first Search
# link: <https://leetcode.com/problems/detect-cycles-in-2d-grid/description/>
# hidden: true
# ---
#
# Given a 2D array of characters `grid` of size `m x n`, you need to find if
# there exists any cycle consisting of the **same value**  in `grid`.
#
# A cycle is a path of **length 4  or more** in the grid that starts and ends at
# the same cell. From a given cell, you can move to one of the cells adjacent to
# it - in one of the four directions (up, down, left, or right), if it has the
# **same value** of the current cell.
#
# Also, you cannot move to the cell that you visited in your last move. For
# example, the cycle `(1, 1) -> (1, 2) -> (1, 1)` is invalid because from `(1,
# 2)` we visited `(1, 1)` which was the last visited cell.
#
# Return `true` if any cycle of the same value exists in `grid`, otherwise,
# return `false`.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/07/15/1.png)**
#
#
#
#     Input: grid = [["a","a","a","a"],["a","b","b","a"],["a","b","b","a"],["a","a","a","a"]]
#     Output: true
#     Explanation: There are two valid cycles shown in different colors in the image below:
#     ![](https://assets.leetcode.com/uploads/2020/07/15/11.png)
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/07/15/22.png)**
#
#
#
#     Input: grid = [["c","c","c","a"],["c","d","c","c"],["c","c","e","c"],["f","c","c","c"]]
#     Output: true
#     Explanation: There is only one valid cycle highlighted in the image below:
#     ![](https://assets.leetcode.com/uploads/2020/07/15/2.png)
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2020/07/15/3.png)**
#
#
#
#     Input: grid = [["a","b","b"],["b","z","b"],["b","b","a"]]
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m <= 500`
#   * `1 <= n <= 500`
#   * `grid` consists only of lowercase English letters.
#
#
## @lc code=start
using LeetCode

function contains_cycle(grid::Vector{Vector{String}})
    dirs = ((1, 0, 3), (-1, 0, 0), (0, 1, 1), (0, -1, 2))
    function bfs(grid::Vector{Vector{String}}, r::Int, c::Int)
        sym = grid[r][c]
        q_pos = Queue{Pair{Int, Int}}()
        q_dir = Queue{Int}()
        enqueue!(q_pos, Pair(r, c))
        enqueue!(q_dir, 0)
        visited[r, c] = true
        while !isempty(q_pos)
            pos = dequeue!(q_pos)
            dir = dequeue!(q_dir)
            for d in dirs
                nx, ny = pos.first + d[1], pos.second + d[2]
                if 0 < nx <= nr && 0 < ny <= nc && dir != d[3] && sym == grid[nx][ny]
                    if visited[nx, ny]
                        return true
                    else
                        visited[nx, ny] = true
                        enqueue!(q_pos, Pair(nx, ny))
                        enqueue!(q_dir, 3 - d[3])
                    end
                end
            end
        end
        false
    end
    nr, nc = length(grid), length(grid[1])
    visited = fill(false, nr, nc)
    for r in 1 : nr
        for c in 1 : nc
            if !visited[r, c] && bfs(grid, r, c)
                ## println(grid[r][c])
                return true
            end
        end
    end
    return false
end
## @lc code=end


text:
### CONCRETE TYPE: DIRECT PROX EVALUATION
# prox! is computed using a Cholesky factorization of A'A + I/(lambda*gamma)
# or AA' + I/(lambda*gamma), according to which matrix is smaller.
# The factorization is cached and recomputed whenever gamma changes

using LinearAlgebra
using SparseArrays
using SuiteSparse

mutable struct LeastSquaresDirect{R <: Real, C <: RealOrComplex{R}, M <: AbstractMatrix{C}, V <: AbstractVector{C}, F <: Factorization} <: LeastSquares
    A::M # m-by-n matrix
    b::V
    Atb::V
    lambda::R
    gamma::R
    shape::Symbol
    S::M
    res::Vector{C} # m-sized buffer
    q::Vector{C} # n-sized buffer
    fact::F
    function LeastSquaresDirect{R, C, M, V, F}(A::M, b::V, lambda::R) where {R <: Real, C <: RealOrComplex{R}, M <: AbstractMatrix{C}, V <: AbstractVector{C}, F <: Factorization}
        if size(A, 1) != length(b)
            error("A and b have incompatible dimensions")
        end
        if lambda <= 0
            error("lambda must be positive")
        end
        m, n = size(A)
        if m >= n
            S = A'*A
            shape = :Tall
        else
            S = A*A'
            shape = :Fat
        end
        new(A, b, A'*b, lambda, -1, shape, S, zeros(C, m), zeros(C, n))
    end
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {R <: Real, C <: Union{R, Complex{R}}, M <: DenseMatrix{C}, V <: AbstractVector{C}}
    LeastSquaresDirect{R, C, M, V, Cholesky{C, M}}(A, b, lambda)
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {R <: Real, C <: Union{R, Complex{R}}, I <: Integer, M <: SparseMatrixCSC{C, I}, V <: AbstractVector{C}}
    LeastSquaresDirect{R, C, M, V, SuiteSparse.CHOLMOD.Factor{C}}(A, b, lambda)
end

# Adjoint/Transpose versions
function LeastSquaresDirect(A::M, b::V, lambda::R) where {R <: Real, C <: Union{R, Complex{R}}, M <: TransposeOrAdjoint{<:DenseMatrix{C}}, V <: AbstractVector{C}}
    LeastSquaresDirect(copy(A), b, lambda)
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {R <: Real, C <: Union{R, Complex{R}}, I <: Integer, M <: TransposeOrAdjoint{<:SparseMatrixCSC{C, I}}, V <: AbstractVector{C}}
    LeastSquaresDirect(copy(A), b, lambda)
end

function LeastSquaresDirect(A::M, b::V, lambda::R) where {R <: Real, C <: Union{R, Complex{R}}, M <: AbstractMatrix{C}, V <: AbstractVector{C}}
    @warn "Could not infer type of Factorization for $M in LeastSquaresDirect, this type will be type-unstable"
    LeastSquaresDirect{R, C, M, V, Factorization}(A, b, lambda)
end

function (f::LeastSquaresDirect)(x::AbstractVector)
    mul!(f.res, f.A, x)
    f.res .-= f.b
    return (f.lambda/2)*norm(f.res, 2)^2
end

function prox!(y::AbstractVector{C}, f::LeastSquaresDirect{R, C, M, V, F}, x::AbstractVector{C}, gamma::R=R(1)) where {R, C, M, V, F}
    # if gamma different from f.gamma then call factor_step!
    if gamma != f.gamma
        factor_step!(f, gamma)
    end
    solve_step!(y, f, x, gamma)
    mul!(f.res, f.A, y)
    f.res .-= f.b
    return (f.lambda/2)*norm(f.res, 2)^2
end

function factor_step!(f::LeastSquaresDirect{R, C, M, V, F}, gamma::R) where {R, C, M, V, F}
    lamgam = f.lambda*gamma
    f.fact = cholesky(f.S + I/lamgam)
    f.gamma = gamma
end

function factor_step!(f::LeastSquaresDirect{R, C, M, V, F}, gamma::R) where {R, C, M <: SparseMatrixCSC, V, F}
    lamgam = f.lambda*gamma
    f.fact = ldlt(f.S; shift = R(1)/lamgam)
    f.gamma = gamma
end

function solve_step!(y::AbstractVector{C}, f::LeastSquaresDirect{R, C, M, V, F}, x::AbstractVector{C}, gamma::R) where {R, C, M, V, F <: Cholesky{C, M}}
    lamgam = f.lambda*gamma
    f.q .= f.Atb .+ x./lamgam
    # two cases: (1) tall A, (2) fat A
    if f.shape == :Tall
        # y .= f.fact\f.q
        y .= f.q
        LAPACK.trtrs!('U', 'C', 'N', f.fact.factors, y)
        LAPACK.trtrs!('U', 'N', 'N', f.fact.factors, y)
    else # f.shape == :Fat
        # y .= lamgam*(f.q - (f.A'*(f.fact\(f.A*f.q))))
        mul!(f.res, f.A, f.q)
        LAPACK.trtrs!('U', 'C', 'N', f.fact.factors, f.res)
        LAPACK.trtrs!('U', 'N', 'N', f.fact.factors, f.res)
        mul!(y, adjoint(f.A), f.res)
        y .-= f.q
        y .*= -lamgam
    end
end

function solve_step!(y::AbstractVector{C}, f::LeastSquaresDirect{R, C, M, V, F}, x::AbstractVector{C}, gamma::R) where {R, C, M, V, F}
    lamgam = f.lambda*gamma
    f.q .= f.Atb .+ x./lamgam
    # two cases: (1) tall A, (2) fat A
    if f.shape == :Tall
        y .= f.fact\f.q
    else # f.shape == :Fat
        # y .= lamgam*(f.q - (f.A'*(f.fact\(f.A*f.q))))
        mul!(f.res, f.A, f.q)
        f.res .= f.fact\f.res
        mul!(y, adjoint(f.A), f.res)
        y .-= f.q
        y .*= -lamgam
    end
end

function gradient!(y::AbstractVector{C}, f::LeastSquaresDirect{R, C, M, V, F}, x::AbstractVector{C}) where {R, C, M, V, F}
    mul!(f.res, f.A, x)
    f.res .-= f.b
    mul!(y, adjoint(f.A), f.res)
    y .*= f.lambda
    fy = (f.lambda/2)*dot(f.res, f.res)
end

function prox_naive(f::LeastSquaresDirect{R, C}, x::AbstractVector{C}, gamma::R=R(1)) where {R, C <: RealOrComplex{R}}
    lamgam = f.lambda*gamma
    y = (f.A'*f.A + I/lamgam)\(f.Atb + x/lamgam)
    fy = (f.lambda/2)*norm(f.A*y-f.b)^2
    return y, fy
end


text:
struct FuzzyVector
    numbers::Vector{FuzzyNumber}
end

Base.length(FV::FuzzyVector) = size(FV.numbers, 1)
Base.getindex(FV::FuzzyVector, i::Int64) = FV.numbers[i]
Base.show(io::IO, FV::FuzzyVector) = println(io, "fuzzy vector length $(length(FV))")
Base.iterate(FV::FuzzyVector, state=1) = state > length(FV.numbers) ? nothing : FV.numbers[state]

# FUZZY ARITHMETIC ============================================================
function Base.:+(A⃗::FuzzyVector, B⃗::FuzzyVector)
	p = length(A⃗)
	result_vector = Vector{FuzzyNumber}(undef, p)
	for k = 1:p
		result_vector[k] = A⃗[k] + B⃗[k]
	end
	FuzzyVector(result_vector)
end
function Base.:-(A⃗::FuzzyVector)
	p = length(A⃗)
	result_vector = Vector{FuzzyNumber}(undef, p)
	for k = 1:p
		result_vector[k] = -A⃗[k]
	end
	FuzzyVector(result_vector)
end
Base.:-(A⃗::FuzzyVector, B⃗::FuzzyVector) = A⃗ + (-B⃗)
# =============================================================================

function Base.:(==)(A⃗::FuzzyVector, B⃗::FuzzyVector)
    if length(A⃗) != length(B⃗)
        return false
    end
    for i = 1:length(A⃗)
        if A⃗[i] != B⃗[i]
            return false
        end
    end
    return true
end

function Base.:∪(FV::FuzzyVector, vec::Vector{FuzzyNumber})
   append!(FV.numbers, vec)
   return FV
end

function cut(X::FuzzyVector, α::Real)
	p = length(X)
	cuts = Vector{Interval}(undef, p)
	for j = 1:p
		cuts[j] = cut(X[j], α)
	end
	cuts
end

function cut(X::Vector{FuzzyVector}, p::Int, α::Real)
	N = length(X)
	cuts = Vector{Interval}(undef, N)
	for i = 1:N
		cuts[i] = cut(X[i][p], α)
	end
	cuts
end

function cut(X::Vector{FuzzyVector}, α::Real)
	N = length(X)
	p = length(X[1])
	cuts = Vector{Vector{Interval}}(undef, N)
	for i = 1:N
		cuts[i] = [cut(X[i][j], α) for j = 1:p]
	end
	cuts
end

function draw(ax::Axis3, FV::FuzzyVector; step=0.001, colormap=:jet1)
	A₁ = FV[1]
	A₂ = FV[2]
	X = collect(-4:step:4)
	Y = collect(-3:step:3)
	# X = collect(A₁.grades[1][1]:step:A₁.grades[1][2])
	# Y = collect(A₂.grades[1][1]:step:A₂.grades[1][2])
	Z = [min(A₁(x), A₂(y)) for x in X, y in Y]
	surface!(ax, X, Y, Z, colormap=colormap)#, axis=(type=Axis3,))

    # marking peak
	# max_indices = argmax(Z)
	# max_x = X[max_indices[1]]
	# max_y = Y[max_indices[2]]
	# μ = Z[max_indices]
    # annotate!(fig, [(max_x, max_y, μ, (μ, 8, :black, :center))])
end

function draw2d(FV::FuzzyVector; step::Float64=0.01, fig=nothing, c=:jet1, alpha::Real=nothing, marker=nothing, peak_text::Bool=false, font=Plots.font("Times", 8))
    gr(xtickfont=font, ytickfont=font, legendfont=font)
	A₁ = FV[1]
	A₂ = FV[2]
	# x1 = max(A₁.grades[1].left, -3.5)
	# x2 = min(A₁.grades[1].right, 3.5)
	# y1 = max(A₂.grades[1].left, -3.5)
	# y2 = min(A₂.grades[1].right, 3.5)
	x1 = A₁.grades[1].left
	x2 = A₁.grades[1].right
	y1 = A₂.grades[1].left
	y2 = A₂.grades[1].right
	X = collect(x1:step:x2)
	Y = collect(y1:step:y2)
	f(x, y) = min(A₁(x), A₂(y))
	if isnothing(fig)
		fig = Plots.contour(X, Y, f, c=c, aspect_ratio=1.0, seriesalpha=alpha, fill=true)
	else
		fig = Plots.contour(fig, X, Y, f, c=c, aspect_ratio=1.0, seriesalpha=alpha, fill=true)
		fig = Plots.contour(fig, [0, 0], [0, 0], [0, 1], c=c, aspect_ratio=1.0, seriesalpha=1.0, fill=true, colorbar=true)
	end

	# plot peak
	center_x = peak_center(A₁) # or can be centroid()
	center_y = peak_center(A₂) # or can be centroid()
	if !isnothing(marker)
		fig = Plots.scatter(fig, (center_x, center_y), legend=false, m=marker, color=:grey)
	end
	if peak_text
		Plots.annotate!([(center_x, center_y, ("$(round(center_x, digits=2)), $(round(center_y, digits=2))", 8, :black, :center))])
	end

	fig = Plots.plot(fig, dpi=200)
	fig
end

function draw_contour_2d(FVs::Vector{FuzzyVector}; step::Float64=0.01, fig=nothing, contour_colors=[:jet1, :jet1], peak_colors=["black", "red"], alpha::Vector{Float64}=nothing, marker=nothing, peak_text::Bool=false, offset::Real=0, xlim=nothing, ylim=nothing, font=Plots.font("Times", 8))
    gr(xtickfont=font, ytickfont=font, legendfont=font)

	centers = []
	for (i, FV) in enumerate(FVs)
		A₁ = FV[1]
		A₂ = FV[2]
		x1 = A₁.grades[1].left
		x2 = A₁.grades[1].right
		y1 = A₂.grades[1].left
		y2 = A₂.grades[1].right
		X = collect(x1:step:x2)
		Y = collect(y1:step:y2)
		f(x, y) = min(A₁(x), A₂(y))
		center_x = peak_center(A₁) # or can be centroid()
		center_y = peak_center(A₂) # or can be centroid()
		push!(centers, (center_x, center_y))
		if isnothing(fig)
			fig = Plots.contour(X, Y, f, c=contour_colors[i], aspect_ratio=1.0, seriesalpha=alpha[i], xlim=xlim, ylim=ylim, fill=true)
		else
			fig = Plots.contour(fig, X, Y, f, c=contour_colors[i], aspect_ratio=1.0, seriesalpha=alpha[i], fill=true, xlim=xlim, ylim=ylim)
			fig = Plots.contour(fig, [0, 0], [0, 0], [0, 1], c=contour_colors[i], aspect_ratio=1.0, seriesalpha=1.0, fill=true, colorbar=true, xlim=xlim, ylim=ylim)
		end
	end

	# plot peak
	for (i, (center_x, center_y)) in enumerate(centers)
		if !isnothing(marker)
			fig = Plots.scatter(fig, (center_x, center_y), legend=false, m=marker, c=peak_colors[i])
		end
		if peak_text
			# fig = Plots.annotate(fig, [(center_x, center_y+offset, ("$(round(center_x, digits=2)), $(round(center_y, digits=2))", 8, peak_colors[i], :center, "Times"))])
			Plots.annotate!([(center_x, center_y+offset, ("$(round(center_x, digits=2)), $(round(center_y, digits=2))", 8, peak_colors[i], :center, "Times"))])
		end
	end

	fig = Plots.plot(fig, dpi=200)
	fig
end

function draw_3d(FVs::Vector{FuzzyVector}; step::Float64=0.01, fig=nothing, color=:grey, alpha::Float64=nothing, peak_text::Bool=false, offset::Real=0, xlim=nothing, ylim=nothing, font=Plots.font("Times", 8))
	centers = []
	for (i, FV) in enumerate(FVs)
		A₁ = FV[1]
		A₂ = FV[2]
		x1 = A₁.grades[1].left
		x2 = A₁.grades[1].right
		y1 = A₂.grades[1].left
		y2 = A₂.grades[1].right
		X = collect(x1:step:x2)
		Y = collect(y1:step:y2)
		f(x, y) = min(A₁(x), A₂(y))
		center_x = peak_center(A₁) # or can be centroid()
		center_y = peak_center(A₂) # or can be centroid()
		push!(centers, (center_x, center_y))
		fig = Plots.surface(fig, X, Y, f, aspect_ratio=1.0, color=color, alpha=alpha, fill=true, xlim=xlim, ylim=ylim)
		# fig = Plots.surface(fig, [0, 0], [0, 0], [0, 1], c=color, aspect_ratio=1.0, alpha=1.0, fill=true, colorbar=true, xlim=xlim, ylim=ylim)
	end

	# plot peak
	# for (i, (center_x, center_y)) in enumerate(centers)
	# 	if peak_text
	# 		# fig = Plots.annotate(fig, [(center_x, center_y+offset, ("$(round(center_x, digits=2)), $(round(center_y, digits=2))", 8, peak_colors[i], :center, "Times"))])
	# 		Plots.annotate!([(center_x, center_y+offset, ("$(round(center_x, digits=2)), $(round(center_y, digits=2))", 8, peak_colors[i], :center, "Times"))])
	# 	end
	# end

	fig = Plots.plot(fig, dpi=200)
	fig
end

function draw_center_3d(FVs::Vector{FuzzyVector}; step::Float64=0.01, fig=nothing, contour_colors=[:jet1, :jet1], peak_colors=["black", "red"], alpha::Vector{Float64}=nothing, marker=nothing, peak_text::Bool=false, offset::Real=0, xlim=nothing, ylim=nothing, font=Plots.font("Times", 8))
	centers = []
	for (i, FV) in enumerate(FVs)
		A₁ = FV[1]
		A₂ = FV[2]
		x1 = A₁.grades[1].left
		x2 = A₁.grades[1].right
		y1 = A₂.grades[1].left
		y2 = A₂.grades[1].right
		X = collect(x1:step:x2)
		Y = collect(y1:step:y2)
		f(x, y) = min(A₁(x), A₂(y))
		center_x = peak_center(A₁) # or can be centroid()
		center_y = peak_center(A₂) # or can be centroid()
		push!(centers, (center_x, center_y))
		fig = Plots.surface(fig, X, Y, f, c=contour_colors[i], aspect_ratio=1.0, seriesalpha=alpha[i], fill=true, xlim=xlim, ylim=ylim)
		fig = Plots.surface(fig, [0, 0], [0, 0], [0, 1], c=contour_colors[i], aspect_ratio=1.0, seriesalpha=1.0, fill=true, colorbar=true, xlim=xlim, ylim=ylim)
	end

	# plot peak
	# for (i, (center_x, center_y)) in enumerate(centers)
	# 	if !isnothing(marker)
	# 		fig = Plots.scatter(fig, (center_x, center_y), legend=false, m=marker, c=peak_colors[i])
	# 	end
	# 	if peak_text
	# 		# fig = Plots.annotate(fig, [(center_x, center_y+offset, ("$(round(center_x, digits=2)), $(round(center_y, digits=2))", 8, peak_colors[i], :center, "Times"))])
	# 		Plots.annotate!([(center_x, center_y+offset, ("$(round(center_x, digits=2)), $(round(center_y, digits=2))", 8, peak_colors[i], :center, "Times"))])
	# 	end
	# end

	fig = Plots.plot(fig, dpi=200)
	fig
end

function draw2d_makie(ax::Axis, FV::FuzzyVector; step=0.01, colormap=:jet1, alpha=nothing)
	A₁ = FV[1]
	A₂ = FV[2]
	x1 = max(A₁.grades[1].left, -3.5)
	x2 = min(A₁.grades[1].right, 3.5)
	y1 = max(A₂.grades[1].left, -3.5)
	y2 = min(A₂.grades[1].right, 3.5)
	X = collect(x1:step:x2)
	Y = collect(y1:step:y2)
	Z = [min(A₁(x), A₂(y)) for x in X, y in Y]
	contour!(ax, X, Y, Z, colormap=colormap)#, alpha=alpha)
end

function u_uncertainty(FV::FuzzyVector; agg::String="average")
	p = length(FV)
	Us = Vector{Float64}(undef, p)
	for k = 1:p
		Us[k] = u_uncertainty(FV[k])
	end
	if agg == "average"
		return sum(Us) / p
	elseif agg == "max"
		return maximum(Us)
	elseif agg == "sum"
		return sum(Us)
	end
end


text:
using Test
@isdefined(TestHelpers) || include("../TestHelpers.jl")

TestHelpers.@timed_testset "layers/" begin
    include("layers/conv2d.jl")
    include("layers/flatten.jl")
    include("layers/linear.jl")
    include("layers/masked_relu.jl")
    include("layers/pool.jl")
    include("layers/relu.jl")
end


text:
#=
using DelimiteFiles
=#

# NOTE: These simulation setting are identical to that of test-sims-6-5 EXCEPT
# that W_true is different. The W_true here has low abundance (< 10% each) for
# 2 pairs (i.e. 4) features. Thus, they should be harder to recover, and
# uncertainty in R_i should be greater.

# Simulation Name
simname = basename("$(@__DIR__)")
println(simname); flush(stdout)

# NOTE: write to scratchdir
function outdir_suffix(pthin, batchprop, alpha)
  return "pthin$(pthin)-batchprop$(batchprop)-alpha$(alpha)-N2000" # TODO:N2000
end

# NOTE: modify
settings = let
  Z = Matrix{Bool}(readdlm(joinpath(@__DIR__, "Z.txt"), Int, comments=true))
  W = readdlm(joinpath(@__DIR__, "W.txt"), comments=true)
  N = [2000, 2000]
  [Dict(:simname => simname,
        :repfam_dist_scale => 1.0,
        :N => N,
        :Z => Z,
        :W => W,
        :thin_samps => 2,
        :nburn => 1000,
        :nsamps => 2000,
        :Lmcmc => Dict(0 => 2, 1 => 2),
        :Kmcmc => 10,
        :pthin => pthin,
        :alpha => alpha,
        :batchprop => batchprop,
        :dataseed => 1,
        :mcmcseed => 1,
        :outdir_suffix => outdir_suffix(pthin, batchprop, alpha))
   for pthin in [2]  # pthin{2,10} didn't work
   for batchprop in [.8]  # batchprop{.05} didn't work
   for alpha in [sum(N) * 1000.0]]  # alpha{10} didn't work
end


text:
module TestDataFrameRow

using Test, DataFrames
using DataFrames: columns

df = DataFrame(a=Union{Int, Missing}[1, 2, 3, 1, 2, 2],
               b=[2.0, missing, 1.2, 2.0, missing, missing],
               c=["A", "B", "C", "A", "B", missing],
               d=CategoricalArray([:A, missing, :C, :A, missing, :C]))
df2 = DataFrame(a = [1, 2, 3])
sdf = view(df, [5, 3], [3, 1, 2])

@test names(DataFrameRow(df, 1, :)) == [:a, :b, :c, :d]
@test names(DataFrameRow(df, 3, [3, 2])) == [:c, :b]
@test copy(DataFrameRow(df, 3, [3, 2])) == (c = "C", b = 1.2)
@test copy(DataFrameRow(sdf, 2, [3, 2])) == (b = 1.2, a = 3)
@test copy(DataFrameRow(sdf, 2, :)) == (c = "C", a = 3, b = 1.2)
@test DataFrameRow(df, 3, [3, 2]) == df[3, [3, 2]] == view(df, 3, [3, 2])
@test DataFrameRow(sdf, 2, [3, 2]) == sdf[2, [3, 2]] == view(sdf, 2, [3, 2])
@test DataFrameRow(sdf, 2, :) == sdf[2, :] == view(sdf, 2, :)
@test DataFrameRow(df, 3, 2:3) === df[3, 2:3]
@test view(df, 3, 2:3) === df[3, 2:3]
@test_throws ArgumentError DataFrameRow(df, 1, :a)
@test_throws ArgumentError DataFrameRow(df, 1, 1)
@test_throws BoundsError DataFrameRow(df, 1, 1:10)
@test_throws BoundsError DataFrameRow(df, 1, [1:10;])
@test_throws BoundsError DataFrameRow(df, 100, 1:2)
@test_throws BoundsError DataFrameRow(df, 100, [1:2;])
@test_throws BoundsError DataFrameRow(df, 100, :)
@test_throws MethodError DataFrameRow(df, true, 1:2)
if VERSION ≥ v"1.0.0"
    # this test throws a warning on Julia 0.7
    @test_throws ArgumentError DataFrameRow(sdf, true, 1:2)
end

# getindex
r = DataFrameRow(df, 2, :)
@test r[:] === r
@test view(r, :) === r
@test r[3] == "B"
@test_throws BoundsError r[5]
@test view(r, 3)[] == "B"
view(r, 3)[] = "BB"
@test df.c[2] == "BB"
@test_throws MethodError r[true]
@test_throws MethodError view(r, true)
@test copy(r[[:c,:a]]) == (c = "BB", a = 2)
@test copy(view(r, [:c,:a])) == (c = "BB", a = 2)
@test copy(r[[true, false, true, false]]) == (a = 2, c = "BB")
r.c = "B"
@test df.c[2] == "B"
@test_throws BoundsError copy(r[[true, false, true]])

r = DataFrameRow(sdf, 2, [3, 1])
@test r[:] == r
@test view(r, :) === r
@test r[2] == "C"
@test_throws BoundsError r[4]
@test view(r, 2)[] == "C"
view(r, 2)[] = "CC"
@test df.c[3] == "CC"
@test_throws MethodError r[true]
@test_throws MethodError view(r, true)
@test copy(r[[:c,:b]]) == (c = "CC", b = 1.2)
@test copy(view(r, [:c,:b])) == (c = "CC", b = 1.2)
@test copy(view(r, [:c,:b])) == (c = "CC", b = 1.2)
@test copy(r[[false, true]]) == (c = "CC",)
r.c = "C"
@test df.c[3] == "C"

#
# Equality
#
@test DataFrameRow(df, 1, :) != DataFrameRow(df2, 1, :)
@test DataFrameRow(df, 1, [:a]) == DataFrameRow(df2, 1, :)
@test DataFrameRow(df, 1, [:a]) == DataFrameRow(df2, big(1), :)
@test DataFrameRow(df, 1, :) != DataFrameRow(df, 2, :)
@test DataFrameRow(df, 1, :) != DataFrameRow(df, 3, :)
@test DataFrameRow(df, 1, :) == DataFrameRow(df, 4, :)
@test ismissing(DataFrameRow(df, 2, :) == DataFrameRow(df, 5, :))
@test isequal(DataFrameRow(df, 2, :), DataFrameRow(df, 5, :))
@test ismissing(DataFrameRow(df, 2, :) != DataFrameRow(df, 6, :))
@test !isequal(DataFrameRow(df, 2, :), DataFrameRow(df, 6, :))

# isless()
df4 = DataFrame(a=[1, 1, 2, 2, 2, 2, missing, missing],
                b=Union{Float64, Missing}[2.0, 3.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0],
                c=[:B, missing, :A, :C, :D, :D, :A, :A])
@test isless(DataFrameRow(df4, 1, :), DataFrameRow(df4, 2, :))
@test !isless(DataFrameRow(df4, 2, :), DataFrameRow(df4, 1, :))
@test !isless(DataFrameRow(df4, 1, :), DataFrameRow(df4, 1, :))
@test isless(DataFrameRow(df4, 1, :), DataFrameRow(df4, 3, :))
@test !isless(DataFrameRow(df4, 3, :), DataFrameRow(df4, 1, :))
@test isless(DataFrameRow(df4, 3, :), DataFrameRow(df4, 4, :))
@test !isless(DataFrameRow(df4, 4, :), DataFrameRow(df4, 3, :))
@test isless(DataFrameRow(df4, 4, :), DataFrameRow(df4, 5, :))
@test !isless(DataFrameRow(df4, 5, :), DataFrameRow(df4, 4, :))
@test !isless(DataFrameRow(df4, 6, :), DataFrameRow(df4, 5, :))
@test !isless(DataFrameRow(df4, 5, :), DataFrameRow(df4, 6, :))
@test isless(DataFrameRow(df4, 7, :), DataFrameRow(df4, 8, :))
@test !isless(DataFrameRow(df4, 8, :), DataFrameRow(df4, 7, :))

# hashing
@test hash(DataFrameRow(df, 1, :)) != hash(DataFrameRow(df, 2, :))
@test hash(DataFrameRow(df, 1, :)) != hash(DataFrameRow(df, 3, :))
@test hash(DataFrameRow(df, 1, :)) == hash(DataFrameRow(df, 4, :))
@test hash(DataFrameRow(df, 2, :)) == hash(DataFrameRow(df, 5, :))
@test hash(DataFrameRow(df, 2, :)) != hash(DataFrameRow(df, 6, :))

# check that hashrows() function generates the same hashes as DataFrameRow
df_rowhashes, _ = DataFrames.hashrows(Tuple(columns(df)), false)
@test df_rowhashes == [hash(dr) for dr in eachrow(df)]

# test incompatible frames
@test_throws UndefVarError DataFrameRow(df, 1, :) == DataFrameRow(df3, 1, :)

# test RowGroupDict
N = 20
d1 = rand(map(Int64, 1:2), N)
df5 = DataFrame([d1], [:d1])
df6 = DataFrame(d1 = [2,3])

# test_group("group_rows")
gd = DataFrames.group_rows(df5)
@test length(unique(gd.groups)) == 2

# getting groups for the rows of the other frames
@test length(gd[DataFrameRow(df6, 1, :)]) > 0
@test_throws KeyError gd[DataFrameRow(df6, 2, :)]
@test isempty(DataFrames.findrows(gd, df6, (gd.df[1],), (df6[1],), 2))

# grouping empty frame
gd = DataFrames.group_rows(DataFrame(x=Int[]))
@test length(unique(gd.groups)) == 0

# grouping single row
gd = DataFrames.group_rows(df5[1:1,:])
@test length(unique(gd.groups)) == 1

# getproperty, setproperty! and propertynames
r = DataFrameRow(df, 1, :)
@test Base.propertynames(r) == names(df)
@test r.a === 1
@test r.b === 2.0
@test copy(r[[:a,:b]]) === (a=1, b=2.0)

r.a = 2
@test r.a === 2
r.b = 1
@test r.b === 1.0

# keys, values and iteration, size
@test keys(r) == names(df)
@test values(r) == (df[1, 1], df[1, 2], df[1, 3], df[1, 4])
@test collect(pairs(r)) == [:a=>df[1, 1], :b=>df[1, 2], :c=>df[1, 3], :d=>df[1, 4]]

@test haskey(r, :a)
@test haskey(r, 1)
@test !haskey(r, :zzz)
@test !haskey(r, 0)
@test !haskey(r, 1000)
@test_throws ArgumentError haskey(r, true)

x = DataFrame(ones(5,4))
dfr = view(x, 2, 2:3)
@test names(dfr) == names(x)[2:3]
dfr = view(x, 2, [4,2])
@test names(dfr) == names(x)[[4,2]]

x = DataFrame(ones(10,10))
r = x[3, [8, 5, 1, 3]]
@test length(r) == 4
@test lastindex(r) == 4
@test ndims(r) == 1
@test ndims(typeof(r)) == 1
@test size(r) == (4,)
@test size(r, 1) == 4
@test_throws BoundsError size(r, 2)
@test keys(r) == [:x8, :x5, :x1, :x3]
r[:] = 0.0
r[1:2] = 2.0
@test values(r) == (2.0, 2.0, 0.0, 0.0)
@test collect(pairs(r)) == [:x8 => 2.0, :x5 => 2.0, :x1 => 0.0, :x3 => 0.0]

# convert
df = DataFrame(a=[1, missing], b=[2.0, 3.0])
dfr = DataFrameRow(df, 1, :)
@test convert(Vector, dfr)::Vector{Union{Float64, Missing}} == [1.0, 2.0]
@test convert(Vector{Int}, dfr)::Vector{Int} == [1, 2]
@test Vector(dfr)::Vector{Union{Float64, Missing}} == [1.0, 2.0]
@test Vector{Int}(dfr)::Vector{Int} == [1, 2]

# copy
df = DataFrame(a=Union{Int, Missing}[1, 2, 3, 1, 2, 2],
               b=[2.0, missing, 1.2, 2.0, missing, missing],
               c=["A", "B", "C", "A", "B", missing])
@test copy(DataFrameRow(df, 1, :)) == (a = 1, b = 2.0, c = "A")
@test isequal(copy(DataFrameRow(df, 2, :)), (a = 2, b = missing, c = "B"))

# parent and parentindices
@test parent(df[2, :]) === df
@test parentindices(df[2, :]) == (2, Base.OneTo(3))
@test parent(df[1, 1:3]) === df
@test parentindices(df[1, [3,2]]) == (1, [3, 2])
sdf = view(df, [4,3], [:c, :a])
@test parent(sdf[2, :]) === df
@test parentindices(sdf[2, :]) == (3, [3, 1])
@test parent(sdf[1, 1:2]) === df
@test parentindices(sdf[1, [2, 2]]) == (4, [1, 1])

# iteration and collect
ref = ["a", "b", "c"]
df = DataFrame(permutedims(ref))
dfr = df[1, :]
@test Base.IteratorEltype(dfr) == Base.EltypeUnknown()
@test collect(dfr) == ref
@test eltype(collect(dfr)) === String
for (v1, v2) in zip(ref, dfr)
    @test v1 == v2
end
for (i, v) in enumerate(dfr)
    @test v == ref[i]
end
dfr = DataFrame(a=1, b=true, c=1.0)[1,:]
@test eltype(collect(dfr)) === Real

@testset "duplicate column" begin
    df = DataFrame([11:16 21:26 31:36 41:46])
    sdf = view(df, [3,1,4], [3,1,4])
    dfr1 = df[2, [2,2,2]]
    dfr2 = sdf[2, [2,2,2]]
    @test names(dfr1) == fill(:x2, 3)
    @test names(dfr2) == fill(:x1, 3)
    @test values(dfr1) == (22, 22, 22)
    @test values(dfr2) == (11, 11, 11)
    @test dfr1.x2 == 22
    dfr1.x2 = 100
    @test values(dfr1) == (100, 100, 100)
    @test df[2, 2] == 100
    @test_throws KeyError dfr1.x1
    @test dfr2.x1 == 11
    dfr2.x1 = 200
    @test values(dfr2) == (200, 200, 200)
    @test df[1, 1] == 200
    @test_throws KeyError dfr2.x2
end

@testset "conversion and push!" begin
        df = DataFrame(x=1, y=2)

        @test df == DataFrame(df[1, :])
        @test df[1:1, [2,1]] == DataFrame(df[1, [2,1]])
        @test df[1:1, 1:1] == DataFrame(df[1, 1:1])
        @test_throws ArgumentError DataFrame(df[1, [1,1]])

        @test_throws ArgumentError push!(df, df[1, 1:1])
        @test df == DataFrame(x=1, y=2)

        @test_throws ArgumentError push!(df, df[1, [2,2]])
        @test df == DataFrame(x=1, y=2)

        @test_throws ArgumentError push!(df, df[1, [2,1,2]])
        @test df == DataFrame(x=1, y=2)

        @test push!(df, df[1, :]) == DataFrame(x=[1, 1], y=[2, 2])
        @test push!(df, df[1, [2,1]]) == DataFrame(x=[1, 1, 1], y=[2, 2, 2])
end

@testset "show" begin
    df = DataFrame(a=nothing, b=1)

    @test sprint(show, DataFrameRow(df, 1, :)) == """
        DataFrameRow
        │ Row │ a       │ b     │
        │     │ Nothing │ $(Int) │
        ├─────┼─────────┼───────┤
        │ 1   │         │ 1     │"""

    df = DataFrame(a=1:3, b=["a", "b", "c"], c=Int64[1,0,1])
    dfr = df[2, 2:3]

    @test sprint(show, dfr) == """
        DataFrameRow
        │ Row │ b      │ c     │
        │     │ String │ Int64 │
        ├─────┼────────┼───────┤
        │ 2   │ b      │ 0     │"""

    @test sprint(show, "text/html", dfr) == "
DataFrameRow

" *
                               "

	b	c
" *
                               "
	String	Int64

" *
                               "

1 rows × 2 columns

2	" *
                               "b	0

"

    @test sprint(show, "text/latex", dfr) == """
        \\begin{tabular}{r|cc}
        \t& b & c\\\\
        \t\\hline
        \t& String & Int64\\\\
        \t\\hline
        \t2 & b & 0 \\\\
        \\end{tabular}
        """

    @test sprint(show, "text/csv", dfr) == """
        \"b\",\"c\"
        \"b\",0
        """

    @test sprint(show, "text/tab-separated-values", dfr) == """
        \"b\"\t\"c\"
        \"b\"\t0
        """
end

end # module


text:
include("r3frac.jl")
include("r3mv.jl")
include("LatticeVars.jl")
include("AtomicVars.jl")
include("AtomicSpeciesVars.jl")
include("MuffinTinRadialVars.jl")
include("APWLOVars.jl")
include("readspecies.jl")

function main()
    Nspecies = 2
    atsp_vars = AtomicSpeciesVars(Nspecies)
    mtr_vars = MuffinTinRadialVars(Nspecies)
    #apwlo_vars = APWLOVars(Nspecies, mtr_vars.maxlapw) # XXX use lmaxapw instead of maxlapw?
    apwlo_vars = APWLOVars(Nspecies, mtr_vars.lmaxapw)

    readspecies!(1, "DATA_species/Si.in", atsp_vars, mtr_vars, apwlo_vars)
    readspecies!(2, "DATA_species/Pt.in", atsp_vars, mtr_vars, apwlo_vars)

    println(atsp_vars.spsymb)
    println(atsp_vars.nstsp)
end

main()


text:
VERSION >= v"0.6-" && __precompile__(true)

module Strategems
using Base.Dates
using Temporal
using Indicators

export
    # universe definitions
    Universe, gather!, get_overall_index,
    # parameter sets
    ParameterSet, get_n_runs, get_param_combos, get_run_params, generate_dict,
    # indicators
    Indicator, calculate,
    # signals
    Signal, prep_signal, ↑, ↓, @signal,
    # rules
    Rule, @rule, →,
    # portfolios
    Portfolio,#, update_portfolio!,
    # order
    AbstractOrder, MarketOrder, LimitOrder, StopOrder, liquidate, long, buy, short, sell,
    # strategy results
    Results,
    # summary statistic calculations
    cum_pnl,
    # strategies
    Strategy, generate_trades, generate_trades!, backtest, backtest!, optimize, optimize!

include("universe.jl")
include("paramset.jl")
include("indicator.jl")
include("signal.jl")
include("rule.jl")
include("portfolio.jl")
include("order.jl")
include("results.jl")
include("stats.jl")
include("strategy.jl")

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libtirpc_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libtirpc")
JLLWrappers.@generate_main_file("libtirpc", UUID("86ddda40-0562-5640-9cf4-1268408005c8"))
end  # module libtirpc_jll


text:
@testset "Loopless FBA" begin

    model = load_model(model_paths["e_coli_core.json"])

    sol = flux_balance_analysis_dict(
        model,
        GLPK.Optimizer;
        modifications = [add_loopless_constraints()],
    )

    @test isapprox(
        sol["BIOMASS_Ecoli_core_w_GAM"],
        0.8739215069684292,
        atol = TEST_TOLERANCE,
    )
end


text:
module SpinModels

using LinearAlgebra, Statistics
using NNlib

export SpinModel, Pairwise, SRBM, RBM
export paramnames, getparams, veccopy, zerosum!
include("types.jl")

export energy, ∂x_energy, EnergyBuffer, ΔE₁, ΔEnergyBuffer
include("energy.jl")

export ratiomatch, RatioMatchBuffer
include("losses/ratiomatch.jl")

export SpinSampler, acceptrate, currentstate, batchsize, spinstates, spinsnum
export MetropolisHastings, metropolishastings!
export GibbsWithGradients, gibbswithgradients!
export ZanellaMH, zanella_mh!
export EnergyDescent, energydescent!
include("samplers/SpinSamplers.jl")
using .SpinSamplers

end # module


text:
"""
    swrfstar(data::Array{<:Real,2}, target::Array{<:Integer, 1}, m::Signed=-1,
                  dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
                  f_type::String="continuous")::Array{Float64,1}

Compute feature weights using SWRF* algorithm. The f_type argument specifies whether the features are continuous or discrete
and can either have the value of "continuous" or "discrete".

---
# Reference:
- Matthew E. Stokes and Shyam Visweswaran.
Application of a spatially-weighted Relief algorithm for ranking genetic predictors of disease.
BioData mining, 5(1):20–20, Dec 2012. 23198930[pmid].
"""
function swrfstar(data::Array{<:Real,2}, target::Array{<:Integer, 1}, m::Signed=-1,
                  dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
                  f_type::String="continuous")::Array{Float64,1}

    # Initialize feature weights vector.
    weights = zeros(Float64, 1, size(data, 2))

    # Compute vectors of maximum and minimum feature values.
    max_f_vals = maximum(data, dims=1)
    min_f_vals = minimum(data, dims=1)

    # Sample m examples without replacement.
    sample_idxs = StatsBase.sample(1:size(data, 1), if (m==-1) size(data,1) else m end, replace=false)
    if (m == -1) m = size(data, 1) end

    # Compute probabilities of class values in training set.
    classes_map = countmap(target)
    num_samples = length(target)
    p_classes = Array{Float64}(undef, length(keys(classes_map)), 2)
    @inbounds for (idx, p) = enumerate(classes_map)  # Compute matrix representation of probabilities.
        p_classes[idx, :] = [p[1], p[2]/num_samples]
    end

    # Go over sampled indices.
    @inbounds for idx = sample_idxs

        # Get samples from same class and different classes.
        samples_same_class = data[(target .== target[idx]) .& (1:length(target) .!= idx), :]
        samples_other_class = data[target .!= target[idx], :]

        # Compute distances to examples with same class value.
        distances_same = dist_func(data[idx:idx, :], samples_same_class)

        # Get class values of examples with different class value.
        target_other = target[target .!= target[idx]]

        # Compute t and u parameter values.
        t_same = Statistics.mean(distances_same)
        u_same = Statistics.std(distances_same)

        # Compute distances to examples with different class value.
        distances_other = dist_func(data[idx:idx, :], samples_other_class)

        # Compute t and u parameter values.
        t_other = Statistics.mean(distances_other)
        u_other = Statistics.std(distances_other)

        # Compute weights for examples from same class.
        neigh_weights_same = 2.0 ./ (1 .+ exp.(-(t_same.-distances_same)/(u_same/4.0 + eps(Float64))) .+ eps(Float64)) .- 1

        # Compute weights for examples from different class.
        neigh_weights_other = 2.0 ./ (1 .+ exp.(-(t_other.-distances_other)/(u_other/4.0 + eps(Float64))) .+ eps(Float64)) .- 1

        # Get probabilities of classes not equal to class of sampled example.
        p_classes_other = p_classes[p_classes[:, 1] .!= target[idx], 2]

        # Get other classes and compute diff sum weights for examples from different classes.
        classes_other = p_classes[p_classes[:, 1] .!= target[idx], 1]
        p_weights = p_classes_other./(1 .- p_classes[p_classes[:, 1] .== target[idx], 2] .+ eps(Float64))

        # Map weights vector of classes of samples with different class values to construct weights vector.
        weights_map = hcat(classes_other, p_weights)
        weights_mult = [weights_map[findall(weights_map[:, 1] .== t), 2][1] for t in target_other]

        ### Weights Update ###

        if f_type == "continuous"
            # If features continuous.

            # Penalty term
            penalty = sum(neigh_weights_same.*(abs.(data[idx:idx, :] .- samples_same_class)./(max_f_vals .- min_f_vals .+ eps(Float64))), dims=1)

            # Reward term
            reward = sum(weights_mult.*(neigh_weights_other.*(abs.(data[idx:idx, :] .- samples_other_class)./(max_f_vals .- min_f_vals .+ eps(Float64)))), dims=1)

            # Weights update
            weights = weights .- penalty./(m*size(samples_same_class, 1) + eps(Float64)) .+ reward./(m*size(samples_other_class, 1) + eps(Float64))

        elseif f_type == "discrete"
            # If features discrete.

            # Penalty term
            penalty = sum(neigh_weights_same.*Int64.(data[idx:idx, :] .!= samples_same_class), dims=1)

            # Reward term
            reward = sum(neigh_weights_other.*weights_mult .* Int64.(data[idx:idx, :] .!= samples_other_class), dims=1)

            # Weights update
            weights = weights .- penalty./(m*size(samples_same_class, 1) + eps(Float64)) .+ reward./(m*size(samples_other_class, 1) + eps(Float64))

        else
            throw(DomainError(f_type, "f_type can only be equal to \"continuous\" or \"discrete\"."))
        end

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

    end

    # Return computed feature weights.
    return vec(weights)

end


text:
function vKITTI_Conv_RGB(w,x)
    x = relu(conv4(w[1],x;padding=1).+w[2])
    x = relu(conv4(w[3],x;padding=1).+w[4])
    x = pool(x)
    x = relu(conv4(w[5],x;padding=1).+w[6])
    x = relu(conv4(w[7],x;padding=1).+w[8])
    x = pool(x)
    x = relu(conv4(w[9],x;padding=1).+w[10])
    x = relu(conv4(w[11],x;padding=1).+w[12])
    x = relu(conv4(w[13],x;padding=1).+w[14])
    x = relu(conv4(w[15],x;padding=1).+w[16])
    x = pool(x)
    x = relu(conv4(w[17],x;padding=1).+w[18])
    x = relu(conv4(w[19],x;padding=1).+w[20])
    x = relu(conv4(w[21],x;padding=1).+w[22])
    x = relu(conv4(w[23],x;padding=1).+w[24])

    x = relu(conv4(w[25],x).+w[26])
    x = deconv4(w[27],x;padding=1,stride=6).+w[28]
end

function vKITTI_Conv_RGB_w()
    w = Array[
        #Conv64
        xavier(Float32,3,3,3,64),        zeros(Float32,1,1,64,1),
        xavier(Float32,3,3,64,64),       zeros(Float32,1,1,64,1),
        #Conv128
        xavier(Float32,3,3,64,128),      zeros(Float32,1,1,128,1),
        xavier(Float32,3,3,128,128),     zeros(Float32,1,1,128,1),
        #Conv256
        xavier(Float32,3,3,128,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        xavier(Float32,3,3,256,256),     zeros(Float32,1,1,256,1),
        #Conv512
        xavier(Float32,3,3,256,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        xavier(Float32,3,3,512,512),     zeros(Float32,1,1,512,1),
        #Scores
        xavier(Float32,1,1,512,13),      zeros(Float32,1,1,13,1),
        #Deconv
        bilinear(Float32,32,32,13,13),   zeros(Float32,1,1,13,1),
    ]
    w[1:24] = load(pretrainedPath("VGG19.jld"),"w")[1:24]
    return w
end

function vKITTI_Conv_RGB_loss(w,x,y)
    pixelwiseSoftloss(Array(vKITTI_Conv_RGB(w,x)),y)
end

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
using Plots
using LaTeXStrings

# Probability of v±dv
p(v;m=1.,R=8.3145,T=298.15) = m/(R*T) * v * exp(-m*v^2/(2R*T))

# Probability of v < x, given m and T
pcum(x;m = 1., T = 298.15, R = 8.3145) = -exp(-m*x^2/(2*R*T)) + 1

# Inverse Cumulative Probability of v < x, given m and T
pinv(x;m=1., T=298.15, R=8.3145) = sqrt(-(2*R*T/m)*log(1-x))

# Start plot
plot(layout=(1,2))

# Collect x
x = collect(0:0.1:200)

# Set y with m = 1, T = 298.15K, R = 8.3145 J/Kmol
y = pcum.(x)

# plot
plot!(x,y,label="",subplot=1)
plot!(xlabel=L"x / m s^{-1}",ylabel=L"p~(v<x)",subplot=1)

xl = 70
yl = pcum(xl)
plot!([-10,xl],[yl,yl],label="",color="black",linestyle=:dash,subplot=1)
plot!([xl,xl],[-10,yl],label="",color="black",linestyle=:dash,subplot=1)
scatter!([xl],[yl],label="",color="black",markerstyle=:circle,subplot=1)

annotate!(xl,-0.07,text(L"v~(p_x)",8),subplot=1)
annotate!(-12,yl,text(L"p_x",8),subplot=1)

plot!(xlim=[0,200],ylim=[0,1],subplot=1)

# Sample points from inverse of cumulative distribution
N = 20_000
x = rand(N)
histogram!(pinv.(x),alpha=0.3,label="",normalize=:probability,bins=1:2:200,subplot=2)

# Plot distribution
x = collect(0:2:200)
plot!(x,2*p.(x),label="",linewidth=2,subplot=2)

plot!(xlabel=L"v / ms^{-1}",ylabel=L"p~(v)dv",subplot=2)
plot!(framestyle=:box,subplot=2)

plot!(size=(800,400),framestyle=:box)

savefig("./vsample.pdf")