text:
#=
Functions to compute likelihood based on DataFrames or JuliaDB tables
=#

# find a prefix that is not used by any columns in the table by appending
# numbers to prefix
function find_unused_prefix(table, prefix::String)
    names = string.(Tables.columnnames(table))
    return find_unused_prefix(names, prefix)
end

# function find_unused_prefix(table::JuliaDB.AbstractIndexedTable, prefix::String)
#     names = string.(colnames(table))
#     return find_unused_prefix(names, prefix)
# end

function find_unused_prefix(names::Union{Vector{String}, NTuple{N, String}}, prefix::String) where N
    name = prefix
    idx = 0
    while any(startswith.(names, [name]))
        name = @sprintf "%s_%d" prefix idx
        idx += 1
    end
    return Symbol(name)
end

# common code to process availability into a vector of which column
# indicates availability for each numbered choice
function index_availability(availability, alt_numbers)
    if isnothing(availability)
        avail_cols = nothing
    else
        avail_cols = convert(Vector{Union{Missing, keytype(alt_numbers)}}, fill(missing, length(alt_numbers)))
        for (alt, avail_col) in availability
            idx = alt_numbers[alt]
            ismissing(avail_cols[idx]) || error("availability for $alt multiply defined, or duplicate alt indices")
            avail_cols[idx] = avail_col
        end

        any(ismissing.(avail_cols)) && error("incomplete availability matrix")
    end

    return avail_cols
end

# TODO change signature to generic Tables.jl table
function prepare_data(table::DataFrame, chosen, alt_numbers, availability)
    # copying input data, not great
    output_table = DataFrame(table)

    # find an unused name
    choice_col = find_unused_prefix(output_table, "enumerated_choice")

    output_table[!, choice_col] = get.([alt_numbers], output_table[:, chosen], [-1])
    any(output_table[!, choice_col] .== -1) && error("not all alternatives appear in utility functions")

    avail_cols = index_availability(availability, alt_numbers)
    !isnothing(avail_cols) && check_availability(output_table, avail_cols, choice_col)

    return output_table, choice_col, avail_cols
end

function prepare_data(table::DTable, chosen, alt_numbers, availability)
    # find an unused name
    choice_col = find_unused_prefix(table, "enumerated_choice")

    table = map(table) do r
        merge(NamedTuple{(choice_col,), Tuple{Int64}}(alt_numbers[r[chosen]]), pairs(r))
    end

    avail_cols = index_availability(availability, alt_numbers)
    !isnothing(avail_cols) && check_availability(table, avail_cols, choice_col)

    return table, choice_col, avail_cols
end

function check_availability(table, avail_cols, choice_col)
    for (i, row) in enumerate(Tables.rows(table))
        reduce(|, map(c -> row[c], avail_cols)) || @error "At row $i, no alternatives are available"
        row[avail_cols[row[choice_col]]] || @error "At row $i, chosen alternative is not available"
    end
end

#=
Version of prepare_data to work with JuliaDB IndexedTables
Does exactly the same as the above function, using multiple dispatch
to seamlessly handle different data sources. Should work out of core
automatically on a Julia cluster.
function prepare_data(table::JuliaDB.AbstractIndexedTable, chosen, alt_numbers, availability)
    # find an unused name
    choice_col = find_unused_prefix(table, "enumerated_choice")

    # TODO ensure that alt_numbers, chosen are available on remote workers
    # smaller concern - make sure they don't leak memory by remaining in memory
    # (they're small, so not a huge concern).
    # TODO does this involve a lot of copying?
    # TODO probably don't need an Int64 here
    # Constructing NamedTuple manually since the field names are determined at runtime
    output_table = @transform table NamedTuple{(choice_col,), Tuple{Int64}}((alt_numbers[cols(chosen)],))
    reduce(min, output_table; select=(choice_col => Int64)) < 0 && error("not all alternatives appear in utility functions")

    avail_cols = index_availability(availability, alt_numbers)
    return output_table, choice_col, avail_cols
end
=#

function rowwise_loglik(loglik_for_row, table::DataFrame, params::Vector{T}, args...)::T where T <: Number
    #mapreduce(r -> loglik_for_row(r, params, args...), +, Tables.namedtupleiterator(table), init=zero(T))::T

    #loglik_for_row!(f, out, r, params, args...) = out[Threads.threadid()] += f(r, params, args...)

    #reduce(+, asyncmap(r -> loglik_for_row(r, params, args...), Tables.namedtupleiterator(table)); init=zero(T))

    @floop ThreadedEx() for row in Tables.namedtupleiterator(table)
        ll_row = loglik_for_row(row, params, args...)::T
        @reduce(ll += ll_row)
    end

    ll
end

function rowwise_loglik(loglik_for_row, table, params::Vector{T}, args...)::T where T <: Number
    # make the vector the same as the element type of params so ForwardDiff works
    mapreduce(r -> loglik_for_row(r, params, args...), +, Tables.rows(table), init=zero(T))::T
    #@time _rowwise_loglik(loglik_for_row, Tables.rows(table), params, args...)
end

#=
as above, but for a (possibly-distributed) IndexedTable
=#
function rowwise_loglik(loglik_for_row, table::DTable, params::Vector{T}, args...)::T where T <: Number
    # fingers crossed forwarddiff can handle distributed functions, I think it should
    # most function calls use first defn, first call on each worker uses second
    # TODO ensure enclosure here is type-stable
    # reducer(x::T, y::NamedTuple)::T = x + loglik_for_row(y, params, args...)::T
    # reducer(x::NamedTuple, y::NamedTuple)::T = loglik_for_row(x, params, args...)::T + loglik_for_row(y, params, args...)::T
    # reducer(x::T, y::T)::T = (x + y)::T

    ll_parts = Vector{Dagger.EagerThunk}()
    sizehint!(ll_parts, length(table.chunks))
    for chunk in table.chunks
        push!(ll_parts, Dagger.spawn(rowwise_loglik, loglik_for_row, chunk, params, args...))
    end

    fetch(Dagger.spawn(+, ll_parts...))
end


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("..")  # CytofRepFAM
using CytofRepFAM, Random, BSON, Test, Distributions
include("repfamFS_tests.jl")
=#

printstyled("Test fitting repFAM via trained priors...\n", color=:yellow)
@testset "repFAM PT" begin
  config = init_state_const_data(N=[3,2]*500, L=Dict(0=>1, 1=>1),
                                 LMCMC=Dict(0=>3, 1=>3),
                                 mus=Dict(0 => [-2.0], 1 => [2.0]),
                                 sig2=fill(2.0, 2), seed=0,
                                 allow_repeated_Z_columns=false)
  cfs = CytofRepFAM.Model.ConstantsFS(config[:c])
  # cfs.omega_prior = Normal(0, 1)
  dfs = CytofRepFAM.Model.DataFS(config[:d], config[:X])
  sfs = CytofRepFAM.Model.StateFS{Float64}(config[:s], dfs)
  tfs = CytofRepFAM.Model.TunersFS(config[:t], config[:s], config[:X])

  # CytofRepFAM.Model.printConstants(cfs)

  # Fit model.
  # For algo tests
  nmcmc = 200
  nburn = 100

  # For compile tests
  # nmcmc = 3
  # nburn = 3

  Nsum = sum(dfs.data.N)
  alpha = 1.0
  @time out = CytofRepFAM.Model.fit_fs_tp!(sfs, cfs, dfs,
                                           nmcmc=nmcmc, nburn=nburn,
                                           batchprop=0.05,
                                           prior_thin=4,
                                           # iMCMC settings
                                           temper=(alpha + Nsum) / alpha,
                                           anneal=false,
                                           mb_update_burn_prop=0.6,
                                           #
                                           Z_marg_lamgam=1.0,
                                           Z_marg_lamgam_decay_rate=100.0,
                                           # Z_marg_lamgam_min=0.05,
                                           Z_marg_lamgam_min=1.0,
                                           printFreq=1, seed=6,
                                           computedden=true,
                                           computeDIC=true,
                                           computeLPML=true,
                                           thin_dden=2,
                                           time_updates=false,
                                           verbose=4)
  println("Writing Output ...")

  outdir = "results/repfam-tp"
  mkpath(outdir)

  BSON.bson("$(outdir)/out_fs_pt.bson", out)
  BSON.bson("$(outdir)/data_fs_pt.bson", Dict(:simdat => config[:simdat]))
end

#= Test
outdir = "results/repfam-tp"
out = BSON.load(joinpath(outdir, "out_fs_pt.bson"))
simdat = BSON.load(joinpath(outdir, "data_fs_pt.bson"))[:simdat]
samps = out[:samples][1]
extract(sym) = [s[sym] for s in samps]
Zs = extract(:theta__Z)
mean(Zs)
Ws = extract(:theta__W)
mean(Ws)
etas = extract(:theta__eta)

deltas = extract(:theta__delta)
acceptance_rate = length(unique(deltas)) / length(deltas)
mus0s = -[cumsum(d[0]) for d in deltas]
mus1s = [cumsum(d[1]) for d in deltas]
mus0 = hcat(mus0s...)'
mus1 = hcat(mus1s...)'
sig2s = extract(:theta__sig2)
sig2 = hcat(sig2s...)'

using RCall
@rimport rcommon
@rimport grDevices
@rimport graphics
begin
  nburn = out[:nburn]
  grDevices.pdf("$(outdir)/results.pdf");
  rcommon.plotPosts(mus0, cname="mus0");
  rcommon.plotPosts(mus1, cname="mus1");
  rcommon.plotPosts(sig2, cname="sig2");
  graphics.plot(out[:ll][end-nburn:end], main="loglike", type="l",
                xlab="iter", ylab="loglike");
  grDevices.dev_off();
end

include("../runs/PlotUtils/PlotUtils.jl")
PlotUtils.plot_dden(ddens=out[:dden],
                    etas=etas, Ws=Ws, Zs=Zs, sig2s=sig2s, deltas=deltas,
                    ygrid=out[:c].constants.y_grid,
                    imgdir=outdir, simdat=simdat)
=#

ext: jl


text:
function inverse(side::HorizonSide)
    return HorizonSide(mod(Int(side) + 2, 4))
end

function f()
    steps1 = 0
    steps2 = 0
    while !isborder(robot, Nord)
        move!(robot, Nord)
        steps1 += 1
    end
    while !isborder(robot, West)
        move!(robot, West)
        steps2 += 1
    end
    putmarker!(robot)
    side = Ost
    while !isborder(robot, Sud)
        putmarker!(robot)
        while !isborder(robot, side)
            move!(robot, side)
            putmarker!(robot)
        end
        side = inverse(side)
        move!(robot, Sud)
    end
    putmarker!(robot)
    while !isborder(robot, side)
        move!(robot, side)
        putmarker!(robot)
    end
    side = inverse(side)
    while !isborder(robot, West)
        move!(robot, West)
    end
    while !isborder(robot, Nord)
        move!(robot, Nord)
    end
    while steps1 > 0
        move!(robot, Sud)
        steps1 -= 1
    end
    while steps2 > 0
        move!(robot, Ost)
        steps2 -= 1
    end
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
using Test
using Expronicon
project_dir = dirname(dirname(pathof(Expronicon)))
test_dir = joinpath(project_dir, "test")
build_dir = joinpath(test_dir, "build")
rm(build_dir; force=true, recursive=true)

cd(project_dir) do
    expand_project(;
        mod=Expronicon,
        build_dir=build_dir,
        uuid="55351af7-c7e9-48d6-89ff-24e801d99491",
        macronames=[Symbol("@match"), Symbol("@switch"), Symbol("@λ")],
        exclude_src=["match.jl", "expand.jl", "patches.jl"],
        src_dont_touch=["types.jl", "codegen.jl"],
        exclude_paths=[
            "README.md",
            joinpath("test", "match.jl"),
            joinpath("test", "expand.jl"),
            "generate.jl",
            joinpath("docs", "Manifest.toml"),
            joinpath("docs", "build"),
            "Manifest.toml", "build",
            ".git", "bin", ".vscode",
            ".github",
        ],
        exclude_modules=[:TOML, :Pkg, :MLStyle]
    )
end

cd(build_dir) do
    run(`$(Base.julia_cmd()) --project -e 'using Pkg; Pkg.test()'`)
end

rm(build_dir; force=true, recursive=true)

options = ExpandOptions(
    mod=Expronicon,
    uuid="55351af7-c7e9-48d6-89ff-24e801d99491",
    macronames=[Symbol("@match"), Symbol("@switch"), Symbol("@λ")],
    exclude_src=["match.jl", "expand.jl", "patches.jl"],
    src_dont_touch=["types.jl", "codegen.jl"],
    exclude_paths=[
        "README.md",
        joinpath("test", "match.jl"),
        joinpath("test", "expand.jl"),
        "generate.jl",
        joinpath("docs", "Manifest.toml"),
        joinpath("docs", "build"),
        "Manifest.toml", "build",
        ".git", "bin", ".vscode",
        ".github",
    ],
    exclude_modules=[:TOML, :Pkg, :MLStyle]
)


text:
struct TextLine{I, S, F <: AbstractFont} <: AbstractShape
    position::Point{I}
    text::S
    font::F
end

is_valid(shape::TextLine) = all(x -> has_char(shape.font, x), shape.text)

get_drawing_optimization_style(::TextLine) = PUT_PIXEL

function draw!(f::F, image::AbstractMatrix, shape::TextLine, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    position = shape.position
    text = shape.text
    font = shape.font

    bitmap = font.bitmap

    height = size(bitmap, 1)
    width = size(bitmap, 2)

    char_position = position

    for char in text
        draw!(f, image, Character(char_position, char, font), color)
        char_position = Point(char_position.i, char_position.j + width)
    end

    return nothing
end


text:
module InputParser

using TypesParticles, TypesBasis
using ReadBasis

export getcalcsetup,parseinputfile

"""
    Use: getcalcsetup(filename::String)

Input:
    filename - input file name and location

Return:
    system - TypesParticles.SystemOfAtoms for calculation
    basis - TypesBasis for calculation

Description: this function is called by the main driver program HatreeFock.jl to
setup the calculation.

TODO: Unit test functions needed

""" function getcalcsetup(filename::String)
    results = parseinputfile(filename);

    system = SystemOfAtoms(results["natoms"],
                           results["positions"],
                           results["cell"]);
    btype = results["basistype"];
    if btype == "6-31"
        basis = retreadbasis(results["basisfile"],btype=btype);
    else
        error("Basis type $btype is not implemented.")
    end

    return (system,basis)
end #getcalcsetup

"""
        Use: parseinputfile(filename)
Inputs:
        filename - input file name and location
Returns:
        parserresults - dictionary of keywords => function evaluations

Description: This function parses the input file and then calls the respective keyword
get function to parse the line. The function returns a dictionary whos keys are
the input keywords and values are the results from the function call.

Input syntax: keyword list-of-commands
                   & - line read continuation character
                   # - line read comment character, included in function call return.

Notes: All errors associated with keyword lines are produced in get_xxxx calls

TODOs: unit test to ensure all gotchas

""" function parseinputfile(filename::String)

    if ! isfile(filename)
        error("The file $filename does not exist")
    end

    parseresults = Dict();
    infile = open(filename,"r")
    keywords = getreqkeywords();
    comments=Array{String,1}();

    while !eof(infile)
        line = readline(infile);
        if (occursin("#",line) || isempty(line))
            push!(comments,line);
        else
            for k in keywords
            #Construct function name to call based on keyword
            funcname = Symbol(:get,k);

            if occursin(k,line)
                linestring = split(line,"$k ")[end];

                # Read keyword multiple lines using \& symbol
                while occursin(r" * &\s?$",linestring)
                    nextline=readline(infile);
                    linestring = linestring*nextline;
                end

                #CHECK - do I want the line below, requires complex string parsing if so.
                linestring=replace(linestring, "&" => "");

                #Evaluate and insert results into dictionary
                parseresults[k] = eval(funcname)(linestring);
            end
            end #keywords
        end
    end #eof

   #Need to loop over any required keywords that were not found in file.
   #TODO - I assume a more clever approach could be used  here.
    for k in keywords
        try
            parseresults[k]
        catch
            funcname = Symbol(:get,k);
            parseresults[k] = eval(funcname)("");
        end
    end

    parseresults["comments"] = comments;
    close(infile)

    return parseresults

end #parseinputfile

"""
        Use: get_cell(line)

Inputs:
    line - String that can be unpacked to a (3,1) Array
    default - Specify if keyword is found was found in read file.
Returns:
    default = [10.0; 10.0; 10.0] or
    3-element Array{Float64,1}

Description: parse simulation x, y, and z lengths in Angstroms. Only
supports cubic cells.

""" function getcell(line::String)
    splitline = split(line);

    if isempty(line)
        println("NOTICE: default cell 10.0 10.0 10.0 is used!")
        return [10.00e0 0.00e0 0.00e0;
                0.00e0 10.00e0 0.00e0;
                0.00e0 0.00e0 10.00e0]
     elseif length(splitline) != 3
            error("The length of the parsed cell string does not correspond to 3-element array.")
    else
        #Only supports cubic cell
        diag = parse.(Float64,split(line))
        cell = zeros(3,3);
        cell[1,1],cell[2,2],cell[3,3] = diag[1], diag[2], diag[3];
        return cell
    end
end #get_cell

"""
        Use: get_positions(lines::String
Inputs: lines - is the atomic number and positions for 1 or more atoms
Returns: atoms - Array of data type Atom

""" function getpositions(lines::String)

    nslots = 4;
    atoms = Array{Atom,1}();
    if isempty(lines)
        println("NOTICE: using default position 5.0 5.0 5.0 for 1 hydrogen atom!")
        return push!(atoms,Atom(5.00e0,5.00e0,5.00e0,1))
    else
        splitlines = split(lines);
        natoms = Int(length(splitlines)/nslots);
        nentries = reshape(splitlines,(nslots,natoms));
        for n=1:natoms
            Z = parse(Int,nentries[1,n]);
            x,y,z = parse.(Float64,nentries[2:4,n]);
            push!(atoms,Atom(x,y,z,Z));
        end #natoms
        return atoms
    end

end #get_positions

"""
        Use: get_natoms(line::String)
Inputs: line - number of atoms
Returns: Integer data type of line

""" function getnatoms(line::String)
    if isempty(line)
        println("NOTICE: using 1 atom as default!")
        return 1
    elseif length(replace(line," " => "")) != 1
        println(line)
        error("The length of the parsed natoms string does not correspond to a single value.")
    else
        return parse(Int,line)
    end

end #get_natoms

"""
returns allowed keywords in input file.
""" function getreqkeywords()
    keywords = Set{String}(["basisfile",
                            "basistype",
                            "cell",
                            "natoms",
                            "positions"
                            ])
end #getkeywords

#Inline functions
getbasisfile(line::String)=line;
getbasistype(line::String)=replace(line," " => "");

end #InputParser


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

# ---------------------------------------------------------
# methods for active clocks
# ---------------------------------------------------------
tau(ac::ActiveClock) = tau(ac.clock)
busy(ac::ActiveClock) = ac.clock.state == Busy()
sync!(ac::ActiveClock, clk::Clock) = sync!(ac.clock, clk)

delay!(ac::ActiveClock, args...) = delay!(ac.clock, args...)
wait!(ac::ActiveClock, args...; kwargs...) = wait!(ac.clock, args...; kwargs...)
now!(ac::ActiveClock, ex::Action) = now!(ac.clock, ex)

# ---------------------------------------------------------
# state machine operations of active clocks
# ---------------------------------------------------------
"Run an active clock for a given duration."
step!(A::ActiveClock, ::Idle, σ::Run) = do_run!(A.clock, σ.duration)

step!(A::ActiveClock, ::Union{Idle, Busy}, ::Sync) = nothing

step!(A::ActiveClock, ::ClockState, ::Query) = put!(A.back, Response(A))

step!(A::ActiveClock, ::Union{Idle, Busy}, σ::Register) = assign(A, σ.x, A.id)

step!(A::ActiveClock, ::Union{Idle, Busy}, ::Reset) = nothing

step!(A::ActiveClock, q::ClockState, σ::ClockEvent) = error("transition q=$q, σ=$σ not implemented")

"""
    activeClock(ch::Channel)

Operate an active clock on a given channel. This is its event loop.
"""
function activeClock(cmd::Channel, ans::Channel)
    info = take!(cmd).x # get a pointer to the master clock and id
    ac = ActiveClock(Clock(), info[1], cmd, ans, info[2], threadid())
    ac.clock.id = info[2]
    sf = Array{Base.StackTraces.StackFrame,1}[]
    exp = nothing
    ac.clock.state = Idle()
    sync!(ac.clock, ac.master[])

    while true
        try
            σ = take!(cmd)
            if σ isa Stop
                break
            elseif σ isa Diag
                put!(ans, Response((exp, sf)))
            else
                step!(ac, ac.clock.state, σ)
            end
        catch exp
            sf = stacktrace(catch_backtrace())
            @warn "clock $(ac.id), thread $(ac.thread) exception: $exp"
#            throw(exp)
        end
    end
    # stop task
    # close channel
end

# ---------------------------------------------------------
# starting and destroying active clocks
# ---------------------------------------------------------
"""
    start_threads(f::Function, mul::Int=3)

Start a task on each available thread (other than 1).

# Arguments
- `f::Function`: function to start, has to take two channels as arguments,
- `mul::Int=3`: startup multiplication factor,
"""
function start_threads(f::Function) :: Vector{ClockChannel}
    ac = ClockChannel[]
    @threads for i in 1:nthreads()
        if threadid() > 1
            ai = ClockChannel(Ref{Task}(), Channel{ClockEvent}(8), Channel{ClockEvent}(5),
                    threadid(), false)
            ai.ref = Ref(@async f(ai.forth, ai.back))
            push!(ac, ai)
        end
    end
    sort!(ac, by = x->x.thread)
    return ac
end

"""
    fork!(master::Clock)

Establish copies of a clock on all parallel threads and operate them as active
clocks under control of the master clock.

# Arguments
- `master::Clock`: the master clock, must be on thread 1
"""
function fork!(master::Clock)
    if (master.id == 0) && (threadid() == 1)
        if VERSION >= v"1.3"
            if nthreads() > 1
                if isempty(master.ac)
                    master.ac = start_threads(activeClock)           # startup steps
                    for i in eachindex(master.ac)                    # to all active clocks …
                        put!(master.ac[i].forth, Start(( Ref(master), i) )) # send pointer and id
                    end
                else
                    println(stderr, "clock already has $(length(clk.ac)) active clocks!")
                end
            else
                println(stderr, "no parallel threads available!")
            end
        else
            println(stderr, "you have Julia $VERSION, threading is available for ≥ 1.3!")
        end
    else
        println(stderr, "only the master clock on thread 1 can be multiplied!")
    end
end

"""
    collapse!(master::Clock)

Transfer the schedules of the parallel clocks to master and them stop them.

!!! note
    If there are processes on other threads registered to parallel clocks,
    make sure that they aren't needed anymore before calling `collapse`. They
    are not transferred to and cannot be controlled by master.
"""
function collapse!(master::Clock)
    if (master.id == 0) && (threadid() == 1)
        for ac in master.ac
            put!(ac.forth, Query())
            c = take!(ac.back).x.clock
            for (ev, t) in pairs(c.sc.events)    # transfer events to master
                while any(i->i==t, values(master.sc.events)) # in case an event at that time exists
                    t = nextfloat(float(t))                  # increment scheduled time
                end
                master.sc.events[ev] = t
            end
            append!(master.sc.cevents, c.sc.cevents)
            append!(master.sc.samples, c.sc.samples)
            put!(ac.forth, Stop())
        end
        empty!(master.ac)
    else
        println(stderr, "only the master clock on thread 1 can be collapsed!")
    end
end

"""
    PClock(Δt::Number=0.01; t0::Number=0, unit::FreeUnits=NoUnits)

Setup a clock with parallel clocks on all available threads.

# Arguments

- `Δt::Number=0.01`: time increment. For parallel clocks Δt has to be > 0.
    If given Δt ≤ 0 it is set to 0.01.
- `t0::Number=0`: start time for simulation
- `unit::FreeUnits=NoUnits`: clock time unit. Units can be set explicitely by
    setting e.g. `unit=minute` or implicitly by giving Δt as a time or else setting
    t0 to a time, e.g. `t0=60s`.

!!! note
    Processes on multiple threads are possible in Julia ≥ 1.3 and with
    [`JULIA_NUM_THREADS > 1`](https://docs.julialang.org/en/v1/manual/environment-variables/#JULIA_NUM_THREADS-1).
"""
function PClock(Δt::Number=0.01; t0::Number=0, unit::FreeUnits=NoUnits)
    Δt = Δt > 0 ? Δt : 0.01
    clk = Clock(Δt, t0=t0, unit=unit)
    fork!(clk)
    return clk
end

# ---------------------------------------------------------
# control of active clocks
# ---------------------------------------------------------
"""
```
pclock(clk::Clock, id::Int ) :: AbstractClock
pclock(ac::ActiveClock, id::Int ) :: AbstractClock
```
Get a parallel clock to a given clock.

# Arguments
- `master::Clock`: a master clock or
- `ac::ActiveClock`: an active clock,
- `id::Int=threadid()`: thread id, defaults to the caller's current thread.

# Returns
- the master `Clock` if id==0,
- a parallel `ActiveClock` else
"""
function pclock(clk::Clock, id::Int) :: AbstractClock
    if id == 0
        @assert clk.id == 0 "you cannot get master from a local clock!"
        return clk
    elseif id in eachindex(clk.ac)
        put!(clk.ac[id].forth, Query())
        return take!(clk.ac[id].back).x
    else
        println(stderr, "parallel clock $id not available!")
    end
end
function pclock(ac::ActiveClock, id::Int) :: AbstractClock
    if id == ac.clock.id
        return ac
    else
        return pclock(ac.master[], id)
    end
end

"""
    diag(clk::Clock, id::Int)

Return the stacktrace from a parallel clock.

# Arguments
- `clk::Clock`: a master clock,
- `id::Int`: the id of a parallel clock.
"""
function diag(clk::Clock, id::Int)
    if id in eachindex(clk.ac)
        if istaskfailed(clk.ac[id].ref[])
            return clk.ac[id].ref[]
        else
            put!(clk.ac[id].forth, Diag())
            return take!(clk.ac[id].back).x
        end
    else
        println(stderr, "parallel clock $id not available!")
    end
end


text:
function julia_cmd(cmd)
    return `
        $(Base.julia_cmd())
        --color=$(Base.have_color ? "yes" : "no")
        --compiled-modules=$(Base.JLOptions().use_compiled_modules != 0 ? "yes" : "no")
        --history-file=no
        --startup-file=$(Base.JLOptions().startupfile != 2 ? "yes" : "no")
        --code-coverage=$(["none", "user", "all"][1+Base.JLOptions().code_coverage])
        $cmd
    `
end

# Pkg.test runs with --check_bounds=1, forcing all bounds checks.
# This is incompatible with CUDAnative (see JuliaGPU/CUDAnative.jl#98)
if Base.JLOptions().check_bounds == 1
    @warn "Running with --check-bounds=yes, restarting tests."
    file = @__FILE__
    run(julia_cmd(`$file`))
    exit()
end

# Setup environment
using Pkg

branch = haskey(ENV, "GITLAB_CI") ? ENV["CI_COMMIT_REF_NAME"] : nothing
for package in ("GPUArrays", "CuArrays", "DistributedArrays")
    installed = false
    if branch !== nothing
        try
            Pkg.add(PackageSpec(name=package, rev=String(branch)))
            @info "Installed $package from $branch branch"
            installed = true
        catch ex
            @warn "Could not install $package from same branch, trying master branch" exception=ex
        end
    end
    if !installed
        Pkg.add(PackageSpec(name=package, rev="master"))
    end
end

using Test
include("examples.jl")


text:
# © Przemyslaw Szufel, 2018, run under Julia 1.0.0

function do_item(modelC :: ModelConfig,rng::MersenneTwister,worker_grid::Matrix{Int},worker_locations::Dict{Int,Tuple{Int,Int}},p2p::Bool) :: String
    GRID_SIZE=modelC.GRID_SIZE
    status=zeros(Float64,(GRID_SIZE))
    borders=Matrix{BorderData}(undef,GRID_SIZE...)
    leaveLocData=Matrix{LeaveLocData}(undef,GRID_SIZE...)
    dim_i=modelC.dim_i
    dim_j=modelC.dim_j

    avg_asyn_bord = 0.
    max_asyn_bord = 0.

    start = time()
    if p2p
        @sync for ix in modelC.slaves
            (i,j) = worker_locations[ix]
            @async status[i,j] = remotecall_fetch(step_setborder_nodes_p2p!,worker_grid[i,j])
        end
        avg_asyn_bord = mean(status)
        max_asyn_bord = maximum(status)
    else
        @sync for ix in modelC.slaves
            (i,j) = worker_locations[ix]
            @async borders[i,j] = remotecall_fetch(get_border,worker_grid[i,j])
        end

        @sync for ix in modelC.slaves
            (i,j) = worker_locations[ix]
            j_left = j-1
            if j_left < 1; j_left = GRID_SIZE[2]; end
            j_right = j+1
            if j_right > GRID_SIZE[2]; j_right = 1; end
            i_top = i-1
            if i_top < 1; i_top = GRID_SIZE[1]; end
            i_bottom = i+1
            if i_bottom > GRID_SIZE[1]; i_bottom = 1; end

            sendto(worker_grid[i,j],border_data=BorderData( #left right top bottom
                    borders[i,j_left].right,  borders[i,j_right].left,
                    borders[i_top,j].bottom,  borders[i_bottom,j].top),
                corner_data=CornerData( #topleft topright bottomleft bottomright
                    borders[i_top,j_left].bottom[modelC.dim_j],  borders[i_top,j_right].bottom[1],
                    borders[i_bottom,j_left].top[modelC.dim_j],  borders[i_bottom,j_right].top[1])
                    )
            @async status[i,j] = remotecall_fetch(step_setborder_nodes!,worker_grid[i,j])
        end
    end

    @sync for ix in modelC.slaves
        (i,j) = worker_locations[ix]
        @async begin
            leaveLocData[i,j] = remotecall_fetch(step_leave_location!,worker_grid[i,j])
            status[i,j] = leaveLocData[i,j].elapsed_time
        end
    end
    max_leave_loc_time = maximum(status)
    avg_leave_loc_time = mean(status)

    reloc_start_time = time()

    add_new_citizens=Matrix{Vector{Int}}(undef,GRID_SIZE...)

    # all units to be located on the board of each color type
    relocate_counts = zeros(Int,length(modelC.probs_w)-1)
    free_cells_count=zeros(Int,GRID_SIZE)

    for ix in modelC.slaves
        (i,j) = worker_locations[ix]
        relocate_counts += leaveLocData[i,j].relocate_counts
        add_new_citizens[i,j]=zeros(Int,length(modelC.probs_w)-1)
        free_cells_count[i,j]=leaveLocData[i,j].free_cells_count
    end

    while sum(relocate_counts) > 0
        for t in 1:length(relocate_counts)
            if relocate_counts[t] == 0
                continue
            end
            relocated = 0
            total_free_cells_count = sum(free_cells_count)
            for ix in shuffle(rng,modelC.slaves)
                (i,j) = worker_locations[ix]
                ρ = free_cells_count[i,j]/total_free_cells_count
                μ = ρ * relocate_counts[t]
                if μ < 2.0
                    μ = 2.0
                end
                σ = √(μ * (1-ρ))
                val = round(randn(rng)*σ+μ)
                t_rand = max(0, Int(min(val,relocate_counts[t]-relocated,free_cells_count[i,j]))  )
                free_cells_count[i,j] -= t_rand
                add_new_citizens[i,j][t] += t_rand
                relocated += t_rand
            end
            relocate_counts[t] -= relocated
        end
    end
    reloc_time = time() - reloc_start_time

    @sync for ix in modelC.slaves
        (i,j) = worker_locations[ix]

        sendto(worker_grid[i,j],new_citizens=add_new_citizens[i,j])
        @async status[i,j] = remotecall_fetch(step_set_new_location!,worker_grid[i,j])
    end
    max_new_loc_time = maximum(status)
    avg_new_loc_time = mean(status)
    elapsed_time = time()-start
    return string("$GRID_SIZE x $dim_i x $dim_j t: $(@sprintf("%.3f", elapsed_time  )) = loc $(@sprintf("%.3f", reloc_time  ))  avg/max_asyn_bord $(@sprintf("%.3f", avg_asyn_bord  ))/$(@sprintf("%.3f", max_asyn_bord  )) avg/max_leave $(@sprintf("%.3f", avg_leave_loc_time))/$(@sprintf("%.3f", max_leave_loc_time)) avg/max_new $(@sprintf("%.3f", avg_new_loc_time  ))/$(@sprintf("%.3f", max_new_loc_time  )) ")

end
function do_item() :: String
    return do_item(modelC,rng,worker_grid,worker_locations,p2p)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#Copyright 2015-16 Julia Computing Inc. All rights reserved

#This file is used as an example to demo the functionality of this package.

using Distributions
using JuliaInXL
using Logging
Logging.configure(level=DEBUG)

f(u) = exp(-u^2/2)/√(2pi)

const u1=Uniform(-5.0,+5.0)
const u2=Uniform(0.0, 0.5)

function simulate_pt()
    x=rand(u1)
    y=rand(u2)
    y<f(x)
end

function simulate(num::Number)
     hits=0
     for i in 1:num
        hits = hits + simulate_pt()
    end
    return hits/num*(0.5*10)
end

function simulateTime(num)
    tic()
    r = simulate(num)
    t = toc()
    return (r, t)
end

simulateArray(x::Array) = [simulate(i) for i in x]

#process_async([simulate, simulateTime], "tcp://127.0.0.1:9999"; bind=true)


text:
module MobileSimulation

using ModiaLang

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

const depthMax = 6
const enableVisualization = true
const animationFile = string("Mobile", depthMax, ".json")

const rodLength       = 0.1
const rodDiameter     = 0.05*rodLength
const rodMaterial     = "DryWood"
const rodColor        = VisualMaterial(color="Brown", transparency=0.5)
const rodTranslation1 = [0.0,  rodLength/2, 0.0]
const rodTranslation2 = [0.0, -rodLength/2, 0.0]

const sphereDiameter  = 0.3*rodLength
const sphereMaterial  = "DryWood"
const sphereColor     = VisualMaterial(color="Blue")

const barWidth       = 0.3*rodLength
const barThickness   = 0.3*rodLength
const barMaterial    = "DryWood"
const barColor       = VisualMaterial(color="DeepSkyBlue3", transparency=0.5)
const leaveBarLength = 0.8*rodLength

const revDiameter = rodDiameter
const revLength   = 1.2*barThickness
const revColor    = VisualMaterial(color="Red")

const damping = 0.006

Rod = Model(
    frame0 = Object3D(feature = Solid(shape = Cylinder(axis=2, length=rodLength, diameter=rodDiameter),
                                    solidMaterial  = rodMaterial,
                                    visualMaterial = rodColor)),
    frame1 = Object3D(parent = :frame0, translation = rodTranslation1),
    frame2 = Object3D(parent = :frame0, translation = rodTranslation2)
)
Bar = Model(
    L = missing,
#    frame0 = Object3D(feature = Solid(shape = Beam(axis = 1, length = :L, width=barWidth, thickness=barThickness),
    frame0 = Object3D(feature = Solid(shape = Box(lengthX=:(L + barWidth), lengthY=barWidth, lengthZ=barThickness),
                                      solidMaterial  = barMaterial,
                                      visualMaterial = barColor)),
    frame1 = Object3D(parent = :frame0, translation = :[-L/2, 0.0, 0.0]),
    frame2 = Object3D(parent = :frame0, translation = :[ L/2, 0.0, 0.0])
)
RevoluteWithDamping(;obj1, obj2, axis=3, phi_start=0.0) = Model(
    rev     = RevoluteWithFlange(obj1=obj1, obj2=obj2, axis=axis, phi=Var(init=phi_start)),
    cyl     = Object3D(parent  = obj1,
                    feature = Visual(shape = Cylinder(axis=axis, diameter=revDiameter, length=revLength),
                                        visualMaterial = revColor)),
    damper  = Damper | Map(d=damping),
    fixed   = Fixed,
    connect = :[(damper.flange_b, rev.flange),
                (damper.flange_a, fixed.flange)]
)
function barLength(depth)
    if depth<=2
        return leaveBarLength
    end
    return 2*barLength(depth-1)
end
function createMobile(depth)
    if depth == 1
        Model(
            rod    = Rod,
            sphere = Object3D(parent = :(rod.frame0), translation = rodTranslation2,
                              feature = Solid(shape = Sphere(diameter=sphereDiameter),
                                              solidMaterial  = sphereMaterial,
                                              visualMaterial = sphereColor))
        )
    else
        Model(
            rod  = Rod,
            bar  = Bar | Map(L=barLength(depth)),
            sub1 = createMobile(depth-1),
            sub2 = createMobile(depth-1),
            rev0 = RevoluteWithDamping(obj1=:(rod.frame2), obj2=:(bar.frame0)),
            rev1 = RevoluteWithDamping(obj1=:(bar.frame1), obj2=:(sub1.rod.frame1)),
            rev2 = RevoluteWithDamping(obj1=:(bar.frame2), obj2=:(sub2.rod.frame1))
        )
    end
end
Mobile = Model(
    world      = Object3D(feature=Scene(gravityField=UniformGravityField(g=9.81, n=[0, -1, 0]), enableContactDetection=false, nominalLength=0.2*barLength(depthMax), enableVisualization=enableVisualization, animationFile=animationFile)),
    worldFrame = Object3D(parent = :world, feature = Visual(shape = CoordinateSystem(length=rodLength))),
    top        = createMobile(depthMax),
    rev0       = RevoluteWithDamping(obj1 = :world, obj2 = :(top.rod.frame1), phi_start=0.2)
)

mobile = @instantiateModel(buildModia3D(Mobile), unitless=true, log=false, logDetails=false, logModel=false, logStateSelection=false, logCode=false, logExecution=false, logTiming=false)

stopTime = 5.0
tolerance = 1e-5
requiredFinalStates = [-0.009569246139092647, 0.06206606324720822, 0.0010219666026627296, 0.024678616817829752, 0.004218575209138456, 0.028468399940849164, -0.009906941842068384, 0.04296850265297842, 0.0007318009858593783, 0.02017149964438393, 0.003687261686657736, 0.023502975932150394, -0.0448394459612976, 0.026723967721919494, 0.0537435564833799, 0.12331025494664288, 0.05248795026088127, 0.11675479912295675, -0.009264844274058073, 0.06728855030873694, 0.0011992360615969712, 0.025406524349364322, 0.004501037229020791, 0.02917639504873484, -0.00965181608855167, 0.0477780102217653, 0.0008958989126344595, 0.020867265917878713, 0.003946827273570862, 0.024201038339941326, -0.04462421281684384, 0.03746102813853381, 0.05486746017952349, 0.11577048656288175, 0.05363002444778805, 0.11033674247191452, -0.02865814302511715, -0.46249723347164023, 0.02588030023351807, 0.562896850566723, 0.02821667734365895, 0.5623903988440904, -0.009403295780267986, 0.06162802647976917, 0.0010701734591814242, 0.024348108279934208, 0.004222449370467353, 0.028064091478322103, -0.00978017118379629, 0.042937738418303645, 0.0007712669619712274, 0.02000389926557959, 0.0036870480729875752, 0.02329559286742148, -0.04450699270608273, 0.027584416181856048, 0.053512955516200826, 0.1203910664037745, 0.052242780723209266, 0.1148698118645492, -0.00908973848103805, 0.06683238262833087, 0.0012488875262191405, 0.025054556524096706, 0.004506535939791427, 0.02874309504494291, -0.009516262851887764, 0.047741100768557164, 0.0009367831965132963, 0.0206820605278372, 0.003948288304262993, 0.02396953329142854, -0.0442740599719771, 0.03835719168973073, 0.054628209624966825, 0.1126450558374131, 0.053378321393267744, 0.10825519520092068, -0.029248882022186243, -0.45868389528398645, 0.02661598555301995, 0.5585200081021191, 0.02895380855666197, 0.5578643900507765, 0.11057838176009542, -0.4328898557101326, -0.1137607979103064, 0.5023169087131035, -0.11293209024062373, 0.49924136851309536, -0.00949264483432839, 0.06269840391598419, 0.0010594569792564434, 0.024662133036232943, 0.004286395100088604, 0.02844710490443639, -0.009858798739817745, 0.04356190104934767, 0.0007615897059759054, 0.020185339490969664, 0.0037442343747962708, 0.023525248558572666, -0.04490190529834802, 0.02861187090459901, 0.054093472471048135, 0.12121222952180354, 0.05282381817075915, 0.11524012664035016, -0.00918252751200148, 0.06794970842859763, 0.0012381041193024645, 0.02538680700446271, 0.004571342380041327, 0.029148124149157, -0.009598363758120855, 0.048403550630416084, 0.0009271317194699239, 0.020879696023068203, 0.00400620925514733, 0.02421928779134631, -0.044678930789653376, 0.039416625110429974, 0.05521797625599808, 0.11354470175688848, 0.05396904771937769, 0.1087130515149207, -0.029386149298980365, -0.4625396185956002, 0.026768483822850988, 0.5631828688184961, 0.02911172224703885, 0.5626015954247693, -0.009324704028879006, 0.062258112033935976, 0.0011079870888003323, 0.024328333168236824, 0.004290387257149793, 0.028037530012991814, -0.009730064979351052, 0.043535748038149275, 0.0008014306498257068, 0.02001640392846844, 0.0037442737008724775, 0.023314890388120756, -0.04456469081149338, 0.0294911482051167, 0.05385917279095154, 0.11825295160221237, 0.052576102903529355, 0.11332114444600068, -0.00900535540075688, 0.0674898715499478, 0.0012880524022445207, 0.02503113089861741, 0.0045769255979494095, 0.02870901282603822, -0.009460770397544161, 0.04837019570880829, 0.0009683769175080957, 0.020692819349648818, 0.004007904268091745, 0.023984347424616213, -0.04432382727176435, 0.04033007578307317, 0.054974659993453495, 0.11037837031847329, 0.05371447828202559, 0.10659649333068186, -0.029979463641433897, -0.4586694882625391, 0.027508572008452217, 0.5587394516064685, 0.02985302711571344, 0.5580073793505101, 0.10994893558074191, -0.43345695851264765, -0.113030926843198, 0.5031313427136914, -0.11219753434123993, 0.5000032205945514, 0.17883870418449752, -0.3672711553933069, -0.17900124607493148, 0.423658663329383, -0.17837270297492822, 0.424290817579418, 0.011456378191880315, 0.3654620146635912]
simulate!(mobile, stopTime=stopTime, tolerance=tolerance, log=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(mobile, "rev0.rev.phi")

end


text:
struct MaximumLikelihoodEVA{T} <: fittedEVA{T}
    "Extreme value model definition"
    model::T
    "Maximum likelihood estimate"
    θ̂::Vector{Float64}
end

"""
    hessian(model::MaximumLikelihoodEVA)::Array{Float64, 2}

Calculates the Hessian matrix associated with the MaximumLikelihoodEVA model.
"""
function hessian(model::MaximumLikelihoodEVA)::Array{Float64, 2}

    fobj(θ) = -loglike(model.model, θ)
    return ForwardDiff.hessian(fobj, model.θ̂)

end

"""
    parametervar(fm::MaximumLikelihoodEVA)::Array{Float64, 2}

Compute the covariance parameters estimate of the fitted model `fm`.

"""
function parametervar(fm::MaximumLikelihoodEVA)::Array{Float64, 2}

    # Compute the parameters covariance matrix
    V = inv(hessian(fm))

    return V
end

"""
    loglike(fd::MaximumLikelihoodEVA)::Real

Compute the model loglikelihood evaluated at θ̂ if the maximum likelihood method has been used.

"""
function loglike(fm::MaximumLikelihoodEVA)::Real

    ll = loglike(fm.model, fm.θ̂)

    return ll

end

"""
    getdistribution(fittedmodel::MaximumLikelihoodEVA)::Vector{<:Distribution}

Return the fitted distribution in case of stationarity or the vector of fitted distribution in case of non-stationarity.

"""
function getdistribution(fittedmodel::MaximumLikelihoodEVA)::Vector{<:Distribution}

    model = fittedmodel.model
    θ̂ = fittedmodel.θ̂

    res = getdistribution(model, θ̂)

    return res

end

"""
    quantile(fm::MaximumLikelihoodEVA, p::Real)::Vector{<:Real}

Compute the quantile of level `p` from the fitted model by maximum likelihood. In the case of non-stationarity, the effective quantiles are returned.

"""
function quantile(fm::MaximumLikelihoodEVA, p::Real)::Vector{<:Real}

    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    q = quantile(fm.model, fm.θ̂, p)

    return q

end

"""
    quantilevar(fd::MaximumLikelihoodEVA, level::Real)::Vector{<:Real}

Compute the variance of the quantile of level `level` from the fitted model `fm`.

"""
function quantilevar(fm::MaximumLikelihoodEVA, level::Real)::Vector{<:Real}

    θ̂ = fm.θ̂
    H = hessian(fm)

    q = quantile(fm, level)

    V = zeros(length(q))

    for i=1:length(q)

        f(θ::DenseVector) = quantile(fm.model,θ,level)[i]
        Δf(θ::DenseVector) = ForwardDiff.gradient(f, θ)
        G = Δf(θ̂)

        V[i] = G'/H*G

    end

    return V

end

"""
    returnlevel(fm::MaximumLikelihoodEVA{BlockMaxima}, returnPeriod::Real)::ReturnLevel

Compute the return level corresponding to the return period `returnPeriod` from the fitted model `fm`.

"""
function returnlevel(fm::MaximumLikelihoodEVA{BlockMaxima}, returnPeriod::Real)::ReturnLevel

      @assert returnPeriod > zero(returnPeriod) "the return period should be positive."

      # quantile level
      p = 1-1/returnPeriod

      return ReturnLevel(BlockMaximaModel(fm), returnPeriod, quantile(fm, p))

end

function cint(rl::ReturnLevel{MaximumLikelihoodEVA{BlockMaxima}}, confidencelevel::Real=.95)::Vector{Vector{Real}}

      @assert rl.returnperiod > zero(rl.returnperiod) "the return period should be positive."
      @assert zero(confidencelevel)<confidencelevel<one(confidencelevel) "the confidence level should be in (0,1)."

      # quantile level
      p = 1-1/rl.returnperiod

      q = quantile(rl.model.fm, p)

      # Compute the credible interval

      α = (1 - confidencelevel)

      v = quantilevar(rl.model.fm,p)

      qdist = Normal.(q,sqrt.(v))

      a = quantile.(qdist,α/2)
      b = quantile.(qdist,1-α/2)

      return slicematrix(hcat(a,b), dims=2)

end

"""
    returnlevel(fm::MaximumLikelihoodEVA{ThresholdExceedance}, threshold::Real, nobservation::Int,
        nobsperblock::Int, returnPeriod::Real)::ReturnLevel

Compute the return level corresponding to the return period `returnPeriod` from the fitted model `fm`.

The threshold should be a scalar. A varying threshold is not yet implemented.

"""
function returnlevel(fm::MaximumLikelihoodEVA{ThresholdExceedance}, threshold::Real, nobservation::Int,
    nobsperblock::Int, returnPeriod::Real)::ReturnLevel

    @assert returnPeriod > zero(returnPeriod) "the return period should be positive."

    # Exceedance probability
    ζ = length(fm.model.data.value)/nobservation

    # Appropriate quantile level given the probability exceedance and the number of obs per year
    p = 1-1/(returnPeriod * nobsperblock * ζ)

    return ReturnLevel(PeakOverThreshold(fm, threshold, nobservation, nobsperblock),
        returnPeriod, threshold .+ quantile(fm, p))

end

function cint(rl::ReturnLevel{MaximumLikelihoodEVA{ThresholdExceedance}}, confidencelevel::Real=.95)::Vector{Vector{Real}}

    @assert rl.returnperiod > zero(rl.returnperiod) "the return period should be positive."
    @assert zero(confidencelevel)<confidencelevel<one(confidencelevel) "the confidence level should be in (0,1)."

    # Exceedance probability
    ζ = length(rl.model.fm.model.data.value)/rl.model.nobservation

    # Appropriate quantile level given the probability exceedance and the number of obs per year
    p = 1-1/(rl.returnperiod * rl.model.nobsperblock * ζ)

    q = rl.model.threshold .+ quantile(rl.model.fm, p)

    # Compute the credible interval

    α = (1 - confidencelevel)

    # Computing the variance corresponding to ζ
    f(ζ::Real) = Extremes.quantile(rl.model.fm.model,rl.model.fm.θ̂,1-1/(rl.returnperiod * rl.model.nobsperblock * ζ))
    h = 1e-6 # Step for the finite difference approximation of the gradient
    Δf = (f(ζ + h) - f(ζ - h))/2/h
    v₁ = Δf.^2*ζ*(1-ζ)/rl.model.nobservation # Variance correponding to the estimates of ζ

    # This component seems to be forgoten by Coles (2001) in Section 4.4.1

    # Computing the variance corresponding to the other parameters
    v₂ = quantilevar(rl.model.fm, p)

    v = v₁ .+ v₂

    qdist = Normal.(q,sqrt.(v))

    a = quantile.(qdist,α/2)
    b = quantile.(qdist,1-α/2)

    return slicematrix(hcat(a,b), dims=2)

end

"""
    showfittedEVA(io::IO, obj::MaximumLikelihoodEVA; prefix::String = "")

Displays a MaximumLikelihoodEVA with the prefix `prefix` before every line.

"""
function showfittedEVA(io::IO, obj::MaximumLikelihoodEVA; prefix::String = "")

    println(io, prefix, "MaximumLikelihoodEVA")
    println(io, prefix, "model :")
    showEVA(io, obj.model, prefix = prefix*"\t")
    println(io)
    println(io, prefix, "θ̂  :\t", obj.θ̂)

end

"""
    transform(fm::MaximumLikelihoodEVA{BlockMaxima})::MaximumLikelihoodEVA

Transform the fitted model for the original covariate scales.
"""
function transform(fm::MaximumLikelihoodEVA{BlockMaxima})::MaximumLikelihoodEVA

    locationcovstd = fm.model.location.covariate
    logscalecovstd = fm.model.logscale.covariate
    shapecovstd = fm.model.shape.covariate

    locationcov = Extremes.reconstruct.(locationcovstd)
    logscalecov = Extremes.reconstruct.(logscalecovstd)
    shapecov = Extremes.reconstruct.(shapecovstd)

    # Model on the original covariate scale
    model = BlockMaxima(fm.model.data, locationcov = locationcov, logscalecov = logscalecov, shapecov = shapecov)

    # Transformation of the parameter estimates
    θ̂ = deepcopy(fm.θ̂)
    ind = Extremes.paramindex(fm.model)

    for (var, par) in zip([locationcovstd, logscalecovstd, shapecovstd],[:μ, :ϕ, :ξ])
        if !isempty(var)
            a = getfield.(var, :scale)
            b = getfield.(var, :offset)

            # θ̂[ind[par][1]] = fm.θ̂[ind[par][1]] - sum( fm.θ̂[ind[par][1+i]] * b[i]/a[i] for i=1:length(a) )

            for i=1:length(a)
                θ̂[ind[par][1]] -= fm.θ̂[ind[par][1+i]]*b[i]/a[i]
                θ̂[ind[par][1+i]] = fm.θ̂[ind[par][1+i]]/a[i]
            end
        end
    end

    # Contruction of the fittedEVA structure
    return MaximumLikelihoodEVA(model, θ̂)

end

"""
    transform(fm::MaximumLikelihoodEVA{ThresholdExceedance})

Transform the fitted model for the original covariate scales.
"""
function transform(fm::MaximumLikelihoodEVA{ThresholdExceedance})

    logscalecovstd = fm.model.logscale.covariate
    shapecovstd = fm.model.shape.covariate

    logscalecov = Extremes.reconstruct.(logscalecovstd)
    shapecov = Extremes.reconstruct.(shapecovstd)

    # Model on the original covariate scale
    model = ThresholdExceedance(fm.model.data, logscalecov = logscalecov, shapecov = shapecov)

    # Transformation of the parameter estimates
    θ̂ = deepcopy(fm.θ̂)
    ind = Extremes.paramindex(fm.model)

    for (var, par) in zip([logscalecovstd, shapecovstd],[:ϕ, :ξ])
        if !isempty(var)
            a = getfield.(var, :scale)
            b = getfield.(var, :offset)

            # θ̂[ind[par][1]] = fm.θ̂[ind[par][1]] - sum( fm.θ̂[ind[par][1+i]] * b[i]/a[i] for i=1:length(a) )

            for i=1:length(a)
                θ̂[ind[par][1]] -= fm.θ̂[ind[par][1+i]]*b[i]/a[i]
                θ̂[ind[par][1+i]] = fm.θ̂[ind[par][1+i]]/a[i]
            end
        end
    end

    # Contruction of the fittedEVA structure
    return MaximumLikelihoodEVA(model, θ̂)

end

function cint(fm::MaximumLikelihoodEVA, clevel::Real=.95)::Array{Array{Float64,1},1}

    @assert 0<clevel<1 "the confidence level should be between 0 and 1."

    α = 1 - clevel

    V = parametervar(fm)

    confint = Vector{Vector{Float64}}()

    q = quantile.(Normal(0,1),[α/2, 1 - α/2])

    for i in eachindex(fm.θ̂)
        push!(confint, fm.θ̂[i] .+ q*sqrt(V[i,i]))
    end

    return confint

end


text:
module Simulations

export TimeStepWizard, Simulation, run!, Callback, iteration, stopwatch

import Base: show

using OrderedCollections: OrderedDict
using Oceananigans: AbstractDiagnostic, AbstractOutputWriter, fields

using Oceananigans.Models
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.TimeSteppers
using Oceananigans.Utils

include("callback.jl")
include("time_step_wizard.jl")
include("nan_checker.jl")
include("simulation.jl")
include("run.jl")

end # module


text:
include("window.jl")

myWindow = j3d.Window(800, 600, "JULIA3D")

j3d.Init(myWindow)

while !(j3d.IsClose(myWindow))

    j3d.PollEvent()
    j3d.Clear(myWindow)
    j3d.Display(myWindow)

end

j3d.Destroy(myWindow)


text:
"""
This struct contains the technical data of a battery.
    type: Battery technology type.
    input_active_power_limits: Named tuple of minimum and maximum input power limits.
    output_active_power_limits: Named tuple of minimum and maximum output power limits.
    ramp_limits: Ramp limits (MW/hr).
    storage_capacity: Named tuple of minimum and maximum storage capacity.
    soc: State of charge - Important for first time-step of the iteration.
    efficiency: Charging and discharging efficiencies.
    bus: Name of bus at which the battery is located.
    zone: Name of zone in which the battery is located.
"""

struct BatteryTech
    type::String
    input_active_power_limits::NamedTuple{(:min, :max),Tuple{Float64,Float64}}
    output_active_power_limits::NamedTuple{(:min, :max),Tuple{Float64,Float64}}
    ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    storage_capacity::NamedTuple{(:min, :max),Tuple{Float64,Float64}}
    soc::Float64
    efficiency::NamedTuple{(:in, :out), Tuple{Float64, Float64}}
    bus::Int64
    zone::String
end

get_type(tech::BatteryTech) = tech.type
get_input_active_power_limits(tech::BatteryTech) = tech.input_active_power_limits
get_output_active_power_limits(tech::BatteryTech) = tech.output_active_power_limits
get_ramp_limits(tech::BatteryTech) = tech.ramp_limits
get_storage_capacity(tech::BatteryTech) = tech.storage_capacity
get_soc(tech::BatteryTech) = tech.soc
get_efficiency(tech::BatteryTech) = tech.efficiency
get_bus(tech::BatteryTech) = tech.bus
get_zone(tech::BatteryTech) = tech.zone

"""
This struct contains all the data for battery storage.
    name: Battery's name.
    tech: Struct containing technical data.
    decision_year: The year in which the investment decision is made.
    construction_year: The year in which construction is completed.
    retirement_year: The year in which the generator is retired.
    end_life_year: Age-based retirement year.
    products: Vector containing the data of all products provided by the generator.
    finance_data: Struct containing financial data.
"""
mutable struct BatteryEMIS{T <: BuildPhase} <: StorageEMIS{T}
    name::String
    tech::BatteryTech
    decision_year::Int64
    construction_year::Int64
    retirement_year::Int64
    end_life_year::Int64
    products::Vector{Product}
    finance_data::Finance
end


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds Base
# The underlying abstraction layer for DynamicBounds.
# See https://github.com/PSORLab/DynamicBoundsBase.jl
#############################################################################
# src/utilities.jl
# Defines the utility structure IntegratorStates which stores the current
# condition of an integrator.
#############################################################################

"""
$(TYPEDEF)

A structure that holds the state of the integrator.

$(TYPEDFIELDS)
"""
mutable struct IntegratorStates
    "Has the first point been evaluated"
    first_pnt_eval::Bool
    "Have the box-constraints changed since the last evaluation"
    new_decision_box::Bool
    "Has the decision point changed since the last evaluation"
    new_decision_pnt::Bool
    "The relaxed problem has been updated"
    set_lower_state::Bool
    "The local problem has been updated"
    set_upper_state::Bool
    "How did the integrator terminate"
    termination_status::TerminationStatusCode
end
IntegratorStates() =  IntegratorStates(false, true, true, false, false, RELAXATION_NOT_CALLED)

function seed_duals(::Val{N}, x::AbstractArray{V}) where {V,N}
  seeds = construct_seeds(Partials{N,V})
  duals = [Dual{Nothing}(x[i],seeds[i]) for i in eachindex(x)]
end


text:
#   Author: V. Vitale
#   Feb 2022

include("mps.jl")
using LsqFit

# MPO W-matrix is a 4-index tensor, W[i,j,s,t]
#     s
#     |
#  i -W- j
#     |
#     t
#
# [s,t] act on the local Hilbert space,
# [i,j] act on the virtual bonds

mutable struct MPO <: AbstractTN
  data::Dict
  N::Int
end

MPO() = MPO(Dict(), 0)

function MPO_dot(W::MPO,Q::MPO)
    if isempty(W.data) || isempty(Q.data)
        @warn "Empty MPO."
        return 0
    end
    temp=MPO()
    temp.N=W.N
    for i in 1:W.N
        sW=size(W.data[i])
        sQ=size(Q.data[i])
        @tensor temp.data[i][:] := W.data[i][-1,-3,-5,1]*Q.data[i][-2,-4,1,-6]
        temp.data[i]=reshape(temp.data[i],sW[1]*sQ[1],sW[2]*sQ[2],sW[3],sQ[4])
    end
    return temp
end

Base.:*(W::MPO,Q::MPO)=MPO_dot(W,Q)

function MPO_MPS_dot(W::MPO,Q::MPS)
    if isempty(W.data)
        @warn "Empty MPO."
        #return 0
    end
    if  isempty(Q.data)
        @warn "Empty MPS."
        #return 0
    end

    temp=MPS()
    temp.N=W.N
    for i in 1:W.N
        sW=size(W.data[i])
        sQ=size(Q.data[i])
        @tensor temp.data[i][:] := W.data[i][-1,-4,-3,1]*Q.data[i][-2,1,-5]
        temp.data[i]=reshape(temp.data[i],sW[1]*sQ[1],sW[3],sW[2]*sQ[3])
    end
    return temp
end

Base.:*(W::MPO,A::MPS)=MPO_MPS_dot(W,A)

function trace_MPO(A::MPO)
    if isempty(A.data)
        @warn "Empty MPO."
        return 0
    end
    temp=Base.copy(A.data[1])
    for i in 1:A.N-1
        @tensor temp[-1,-2,-3,-4] := temp[-1,2,1,1]*A.data[i+1][2,-2,-3,-4]
    end
    @tensor result=temp[1,1,2,2]
    return result
end

LinearAlgebra.:tr(A::MPO)= trace_MPO(A)

function adjoint(A::MPO)
    #dagA=copy(A)
    dagA=MPO()
    dagA.N=A.N
    for i in 1:dagA.N
       @tensor dagA.data[i][:] := conj(A.data[i][-1,-2,-4,-3])
    end
    return dagA
end

function Initialize!(s::String,W::MPO,J::Float64,h::Float64,N::Int)
    if s=="Ising"
        d=2
        D=3
        id = [ 1 0 ; 0 1 ]
        sp = [ 0 1 ; 0 0 ]
        sm = [ 0 0 ; 1 0 ]
        sz = [ 1 0 ; 0 -1 ]
        sx = [ 0 1 ; 1 0 ]
        id = [ 1 0 ; 0 1 ]
        Wt = im *  zeros(D,D,d,d)
        W1 = im *  zeros(1,D,d,d)
        W2 = im *  zeros(D,1,d,d)
        Wt[1,1,:,:]=id
        Wt[2,1,:,:]=sz
        Wt[3,1,:,:]=-h*sx
        Wt[3,2,:,:]=-J*sz
        Wt[3,3,:,:]=id

        W1[1,1,:,:]=-h*sx
        W1[1,2,:,:]=-J*sz
        W1[1,3,:,:]=id

        W2[1,1,:,:]=id
        W2[2,1,:,:]=sz
        W2[3,1,:,:]=-h*sx

        W.data[1] = Base.copy(W1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(W2)
        W.N=N
        return "TFIM MPO"
    else
        @warn "Wrong parameters"
    end
end

function Initialize!(s::String,W::MPO,N::Int)
    if s=="Magnetization"
        chi=2
        d=2

        σz = 0.5*[1 0; 0 -1]
        Id2= [1 0; 0 1]
        O2 = [0 0; 0 0]

        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)

        Wt1[1,1,:,:] = σz
        Wt1[1,2,:,:] = Id2
        Wt[1,1,:,:]= Id2
        Wt[1,2,:,:]= O2
        Wt[2,1,:,:]= σz
        Wt[2,2,:,:]= Id2
        Wt2[1,1,:,:] = Id2
        Wt2[2,1,:,:] = σz

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
        return "Magnetization MPO"
    elseif s=="Neel"
        h=100
        d=2
        D=2
        id = [ 1 0 ; 0 1 ]
        sz = [ 1 0 ; 0 -1 ]
        Wp = im *  zeros(D,D,d,d)
        Wd = im *  zeros(D,D,d,d)
        W1 = im *  zeros(1,D,d,d)
        W2 = im *  zeros(D,1,d,d)

        Wp[1,1,:,:]=id
        Wp[2,1,:,:]=-h*sz
        Wp[2,2,:,:]=id
        Wd[1,1,:,:]=id
        Wd[2,1,:,:]=h*sz
        Wd[2,2,:,:]=id

        W1[1,1,:,:]=-h*sz
        W1[1,2,:,:]=id

        W2[1,1,:,:]=id
        W2[2,1,:,:]=h*sz

        W.data[1] = Base.copy(W1)
        for i in 2:2:(N-1)
            W.data[i] = Base.copy(Wd)
            W.data[i+1] = Base.copy(Wp)
        end
        W.data[N] = Base.copy(W2)
        W.N=N
        return "Neel"
    elseif s=="Local_Haar"
        chi=1
        d=2
        dist = Haar(d)
        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)
        Wt[1,1,:,:] =rand(dist, d)
        Wt1[1,1,:,:] = rand(dist, d)
        Wt2[1,1,:,:] = rand(dist, d)
        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
    elseif s=="Id"
        chi=1
        d=2
        id= [1 0; 0 1]
        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)
        Wt[1,1,:,:] = id
        Wt1[1,1,:,:] = id
        Wt2[1,1,:,:] = id

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
    else
        @warn "Wrong parameters"
    end
end

function Initialize!(s::String,W::MPO,d::Int,chi::Int,N::Int)
    if s=="Zeros"
        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
        return "Null MPO"
    elseif s=="Random"
        Wt = im *  rand(chi,chi,d,d)
        Wt1 = im *  rand(1,chi,d,d)
        Wt2 = im *  rand(chi,1,d,d)

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
        return "Random MPO"
    else
        @warn "Wrong parameters"
    end
end

function truncate_MPO(W::MPO,tol::Float64)
    N=length(W)
    for i in 1:N-1
        temp=permutedims(W.data[i],(1,3,4,2))
        sW= size(temp)
        F = qr(reshape(temp,(sW[1]*sW[2]*sW[3],sW[4])))
        q=convert(Matrix{Complex{Float64}}, F.Q)
        r=convert(Matrix{Complex{Float64}}, F.R)
        q = permutedims(reshape(q,(sW[1],sW[2],sW[3],:)),(1, 4, 2, 3))
        W.data[i] = q
        @tensor W.data[i+1][:]:= r[-1,1]*W.data[i+1][1,-2,-3,-4]
    end

    for i in N:-1:2
        temp = permutedims(W.data[i],(1,3,4,2))
        sW = size(temp)
        U,S,V = svd(reshape(temp,(sW[1],sW[2]*sW[3]*sW[4])),full=false)
        s_norm = norm(S)
        S=S/norm(S)
        indices = findall(1 .-cumsum(S.^2) .< tol)
        if length(indices)>0
            chi = indices[1]+1
        else
            chi = size(S)[1]
        end
        if size(S)[1] > chi
            U = U[:,1:chi]
            S = S[1:chi]
            V =  V[:,1:chi]
        end
        V=V'
        S /= norm(S)
        S *= s_norm

        W.data[i] = permutedims(reshape(V,(:,sW[2],sW[3],sW[4])),(1,4,2,3))
        temp=U*diagm(S)
        @tensor W.data[i-1][:] := W.data[i-1][-1,1,-3,-4]*temp[1,-2]
    end
    return W
end

function Initialize!(s::String,W::MPO,alpha::Float64,J::Float64,hz::Float64,k::Int,N::Int)
    chi=k+2
    if s=="LongRangeIsing"
        function expfit(x,p)
            res=0
            for i in 1:div(length(p),2)
                res=res.+p[2*i-1]*exp.(-x./p[2*i])
            end
            return res
        end
        model(x, p) = expfit(x,p)
        xdata=Array(range(1, stop=N,length=2048))
        ydata=xdata.^(-alpha);
        fit =curve_fit(model,xdata,ydata,rand(2*k),lower=1e-9zeros(2*k))
        c = coef(fit)[1:2:end]; λ = coef(fit)[2:2:end]

        σx = [0 1; 1 0]
        σz = [1 0; 0 -1]
        Id2= [1 0; 0 1]

        #Nalpha = 0
        #for i in 1:N
        #    for j in i+1:N
        #        Nalpha += 1/(j-i)^alpha
        #    end
        #end
        #J = J/(N-1)*Nalpha

        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)

        Wt[1, 1, :, :] = Id2
        for i in 1:k
            Wt[1+i, 1, :, :] = σx
            Wt[1+i,1+i,:, :] = exp(-1/λ[i])*Id2
            Wt[end, 1+i,:, :] = -J*exp(-1/λ[i])*c[i]*σx
        end
        Wt[end, 1, :, :] = -hz*σz
        Wt[end,end,:,:] = Id2
        Wt1  = reshape(Wt[end,:,:,:],(1,chi,2,2))
        Wt2 = reshape(Wt[:,1,:,:],(chi,1,2,2))

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
        return "Long Range MPO"
    else
        @warn "Wrong parameters"
    end
end

function Initialize!(s::String,W::MPO,alpha::Float64,J::Float64,hz::Float64,γp::Float64,γm::Float64,γz::Float64,k::Int,N::Int)
    chi=2*k+2
    d=4
    if s=="OpenLongRangeIsing"
        function expfit(x,p)
            res=0
            for i in 1:div(length(p),2)
                res=res.+p[2*i-1]*exp.(-x./p[2*i])
            end
            return res
        end
        model(x, p) = expfit(x,p)
        xdata=Array(range(1, stop=N,length=2048))
        ydata=xdata.^(-alpha);
        fit =curve_fit(model,xdata,ydata,ones(2*k),lower=1e-9zeros(2*k))
        c = coef(fit)[1:2:end]; λ = coef(fit)[2:2:end]
        println(coef(fit))
        σx = [0 1; 1 0]
        σz = [1 0; 0 -1]
        Id2= [1 0; 0 1]
        Id4= kron(Id2,Id2)
        σx1 = kron(σx,Id2)
        σx2 = kron(Id2,σx)
        σz1 = kron(σz,Id2)
        σz2 = kron(Id2,σz)
        Diss= [-(γm+γp)/2 0 0 γp ; 0 -(γp+2*γz) 0 0; 0 0 -(γm +2*γz) 0; γm 0 0 -(γm+γp)/2]

        Nalpha = 0
        for i in 1:N
            for j in i+1:N
                Nalpha += 1/(j-i)^alpha
            end
        end
        J = J/(N-1)*Nalpha

        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)

        Wt[1, 1, :, :] = Id4
        for i in 1:k
            Wt[1+i, 1, :, :] = σx1
            Wt[1+i,1+i,:, :] = exp(-1/λ[i])*Id4
            Wt[end, 1+i,:, :] = -J*exp(-1/λ[i])*c[i]*σx1
            Wt[1+i+k, 1, :, :] = σx2
            Wt[1+i+k,1+i+k,:, :] = exp(-1/λ[i])*Id4
            Wt[end,1+i+k,:, :] = -J*exp(-1/λ[i])*c[i]*σx2
        end

        Wt[end, 1, :, :] += -hz*σz1-hz*σz2+im*Diss
        Wt[end,end,:,:] = Id4

        Wt1  = reshape(Wt[end,:,:,:],(1,chi,d,d))
        Wt2 = reshape(Wt[:,1,:,:],(chi,1,d,d))

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
        return "Open Long Range MPO"
    else
        @warn "Wrong parameters"
    end
end

function Initialize_Brydges!(s::String,W::MPO,N::Int)
    J=[223.394 0.778458; 292.381 0.539682; 208.695 0.182969]

    if s=="Closed"
        d=2
        k=3
        chi=2*k+2
        hz=0.
        c = 0.001*J[:,1]
        λ = J[:,2]

        σp = [0 1; 0 0]
        σm = [0 0; 1 0]
        σz = [1 0; 0 -1]
        Id2= [1 0; 0 1]

        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)

        Wt[1, 1, :, :] = Id2
        for i in 1:k
            Wt[1+i, 1, :, :] = σp
            Wt[1+i,1+i,:, :] = λ[i]*Id2
            Wt[end, 1+i,:, :] = -λ[i]*c[i]*σm

            Wt[k+1+i, 1, :, :] = σm
            Wt[k+1+i,k+1+i,:, :] = λ[i]*Id2
            Wt[end, k+1+i,:, :] = -λ[i]*c[i]*σp
        end
        Wt[end, 1, :, :] = -hz*σz
        Wt[end,end,:,:] = Id2
        Wt1  = reshape(Wt[end,:,:,:],(1,chi,d,d))
        Wt2 = reshape(Wt[:,1,:,:],(chi,1,d,d))

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
        return "Brydges MPO"
    elseif s=="Open"
        k=3
        chi=4*k+2
        hz=0.
        d=4
        γm=(1/1.17)*0.001
        γx=0.69*0.001

        c = 0.001*J[:,1]
        λ = J[:,2]
        σp = [0 1; 0 0]
        σm = [0 0; 1 0]
        σx = [0 1; 1 0]
        σz = [1 0; 0 -1]
        Id2= [1 0; 0 1]

        Id4= kron(Id2,Id2)
        σp1 = kron(σp,Id2)
        σp2 = kron(Id2,σp)
        σm1 = kron(σm,Id2)
        σm2 = kron(Id2,σm)
        σz1 = kron(σz,Id2)
        σz2 = kron(Id2,σz)
        Diss= γm*kron(σm,σm)-0.5*γm*kron(σp*σm,Id2)-0.5*γm*kron(Id2,σp*σm)+γx*kron(σx,σx)-γx*Id4
        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)

        Wt[1, 1, :, :] = Id4
        for i in 1:k
            Wt[1+i, 1, :, :] = σp1
            Wt[1+i,1+i,:, :] = λ[i]*Id4
            Wt[end, 1+i,:, :] = -λ[i]*c[i]*σm1

            Wt[k+1+i, 1, :, :] = σm1
            Wt[k+1+i,k+1+i,:, :] = λ[i]*Id4
            Wt[end, k+1+i,:, :] = -λ[i]*c[i]*σp1

            Wt[2*k+1+i, 1, :, :] = σp2
            Wt[2*k+1+i,2*k+1+i,:, :] = λ[i]*Id4
            Wt[end, 2*k+1+i,:, :] = λ[i]*c[i]*σm2

            Wt[3*k+1+i, 1, :, :] = σm2
            Wt[3*k+1+i,3*k+1+i,:, :] = λ[i]*Id4
            Wt[end, 3*k+1+i,:, :] = λ[i]*c[i]*σp2
        end
        Wt[end, 1, :, :] = im*Diss
        Wt[end,end,:,:] = Id4
        Wt1  = reshape(Wt[end,:,:,:],(1,chi,d,d))
        Wt2 = reshape(Wt[:,1,:,:],(chi,1,d,d))

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
        return "Open Brydges MPO"
    elseif s=="Trajectories"
        γm=(1/1.17)*0.001
        γx=0.69*0.001
        d=2
        k=3
        chi=2*k+2
        hz=0.
        c = 0.001*J[:,1]
        λ = J[:,2]

        σp = [0 1; 0 0]
        σm = [0 0; 1 0]
        σz = [1 0; 0 -1]
        Id2= [1 0; 0 1]

        Wt = im *  zeros(chi,chi,d,d)
        Wt1 = im *  zeros(1,chi,d,d)
        Wt2 = im *  zeros(chi,1,d,d)

        Wt[1, 1, :, :] = Id2
        for i in 1:k
            Wt[1+i, 1, :, :] = σp
            Wt[1+i,1+i,:, :] = λ[i]*Id2
            Wt[end, 1+i,:, :] = -λ[i]*c[i]*σm

            Wt[k+1+i, 1, :, :] = σm
            Wt[k+1+i,k+1+i,:, :] = λ[i]*Id2
            Wt[end, k+1+i,:, :] = -λ[i]*c[i]*σp
        end
        Wt[end, 1, :, :] = -hz*σz-0.5*im*γm*σp*σm
        Wt[end,end,:,:] = Id2
        Wt1  = reshape(Wt[end,:,:,:],(1,chi,d,d))
        Wt2 = reshape(Wt[:,1,:,:],(chi,1,d,d))

        W.N=N
        W.data[1] = Base.copy(Wt1)
        for i in 2:(N-1)
            W.data[i] = Base.copy(Wt)
        end
        W.data[N] = Base.copy(Wt2)
        return "Brydges MPO"
    else
        @warn "Wrong parameters"
    end
end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
using Logging

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

include("../src/ExactODEReduction.jl")
using .ExactODEReduction: load_ODEs, construct_jacobians, many_invariant_subspaces, polynormalize, invariant_subspace_1

# e2.ode
fname = "./data/ODEs/testing/fceri_ji/fceri_ji.ode"

equations = load_ODEs(fname)
As = construct_jacobians(equations)

@info "Loaded system of dimension: " size(As[1])

Vs = many_invariant_subspaces(As, backend_algorithm=invariant_subspace_1)
@info "The number of reductions is $(length(Vs))"

println(equations)

for v in Vs
    @info "-----------------------"
    polyform = polynormalize(v, parent(first(keys(equations))))
    @info "Dimenson $(length(polyform))"
    for p in polyform
        @info p
    end
end


text:
using XGPaint
using Healpix

halo_pos, halo_mass = read_halo_catalog_hdf5(
    "/tigress/zequnl/xgpaint/websky_halos-light.hdf5")

## Load halos from HDF5 files, establish a CIB model and cosmology
cosmo = get_cosmology(h=0.7f0, OmegaM=0.25f0)
radio_model = Radio_Sehgal2009{Float32}(a_0=-1)

##

@time begin
    sources = generate_sources(radio_model, cosmo, halo_pos, halo_mass);
end;

##
m = Map{Float64,RingOrder}(radio_model.nside)
for freq in ["030", "090", "148", "219", "277", "350"]
    @time begin
        paint!(m, parse(Float32, freq) * 1.0f9, radio_model, sources)
        Healpix.saveToFITS(m, "!/tigress/zequnl/xgpaint/jl/radio$(freq).fits")
    end
end
##


text:
# I think this might be in correct
# it claims to be correct by treating the µ vector as a scalar for single vector case
#function ∑matrix(x::Vector)
function sigma(x::Vector)
    y = x * mean(x)
    y * y'
end

function Shapiro_Wilks{T<:Float64}(x::Vector{T})
    sx = sort(x)
    d = Truncated(Normal(mean(xs),1),minimum(xs), maximum(xs))
    m = sort(rand(d, length(x)))
    #V = ∑matrix(sx)
    V = sigma(sx)
    a = m' * inv(V) / sqrt(m' * inv(V) * inv(V) * m)
    y - x.-mean(x)

    # ∑ ai * xi and then ^ 2
    # so a cross product and then square it
    # a'x ^ 2
    num = a'x ^ 2

    # ∑ (xi - xbar)  ^ 2
    # so first get y = xi - xbar vector
    # y'y
    den = y'y
    num/den
end