text:
struct WendlandC6{T} <: SPHKernel
    n_neighbours::Int64
    norm_1D::T
    norm_2D::T
    norm_3D::T
end

"""
    WendlandC6(T::DataType=Float64, n_neighbours::Integer=295)

Set up a `WendlandC6` kernel for a given DataType `T`.
"""
WendlandC6(T::DataType=Float64, n_neighbours::Integer=295) = WendlandC6{T}(n_neighbours, T(55.0/32.0), T(78.0/(7.0*π)), T(1365.0/(64.0*π)))

"""
    kernel_value_1D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

Evaluate WendlandC6 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value_1D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_1D * h_inv
        t1 = 1 - u
        t6 = t1*t1*t1*t1*t1*t1
        u2 = u*u
        return ( t6 * t1 * (1 + 7u + 19u2 + 21u2 * u)) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_1D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

Evaluate the derivative of the WendlandC6 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_deriv_1D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

     if u < 1
        n = kernel.norm_1D * h_inv^2
        t1 = 1 - u
        t6 = t1*t1*t1*t1*t1*t1
        u2 = u*u
        return ( -6t6 * u * (35u2 + 18u + 3 )) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_1D(kernel::WendlandC6, density::Real, m::Real, h_inv::Real)

Corrects the density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
function bias_correction_1D(kernel::WendlandC6{T}, density::Real, m::Real, h_inv::Real) where T

    n = kernel.norm_1D * h_inv^3
    wc_correction = T(0.0116) * ( kernel.n_neighbours * T(0.01) )^T(-2.236) * m * n

    if wc_correction < T(0.2)*density
        density -= wc_correction
    end

    return density |> T
end

"""
    (kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

Evaluate WendlandC6 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value_2D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_2D * h_inv^2
        u_m1 = 1 - u
        u_m1 = u_m1 * u_m1  # (1.0 - u)^2
        u_m1 = u_m1 * u_m1  # (1.0 - u)^4
        u_m1 = u_m1 * u_m1  # (1.0 - u)^8
        u2 = u*u
        return ( u_m1 * ( 1 + 8u + 25u2 + 32u2*u )) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_2D(kernel::WendlandC6, u::Real, h_inv::Real)

Evaluate the derivative of the WendlandC6 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_deriv_2D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_2D * h_inv^3
        u_m1 = 1 - u
        u_m1_7 = u_m1 * u_m1 * u_m1 * u_m1 * u_m1 * u_m1 * u_m1
        return ( -22u_m1_7 * u * ( 16*u^2 + 7u + 1 )) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_2D(kernel::WendlandC6{T}, density::Real, m::Real, h_inv::Real) where T

Corrects the density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
function bias_correction_2D(kernel::WendlandC6{T}, density::Real, m::Real, h_inv::Real) where T

    n = kernel.norm_2D * h_inv^3
    wc_correction = T(0.0116) * ( kernel.n_neighbours * T(0.01) )^T(-2.236) * m * n

    if wc_correction < T(0.2)*density
        density -= wc_correction
    end

    return density |> T
end

"""
    kernel_value_3D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

Evaluate WendlandC6 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_value_3D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_3D * h_inv^3
        u_m1 = 1 - u
        u_m1 = u_m1 * u_m1  # (1.0 - u)^2
        u_m1 = u_m1 * u_m1  # (1.0 - u)^4
        u_m1 = u_m1 * u_m1  # (1.0 - u)^8
        u2 = u*u
        return ( u_m1 * ( 1 + 8u + 25u2 + 32u2*u )) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    kernel_deriv_3D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

Evaluate the derivative of the WendlandC6 spline at position ``u = \\frac{x}{h}``.
"""
function kernel_deriv_3D(kernel::WendlandC6{T}, u::Real, h_inv::Real) where T

    if u < 1
        n = kernel.norm_3D * h_inv^2 * h_inv^2
        u_m1 = 1 - u
        u_m1_7 = u_m1 * u_m1 * u_m1 * u_m1 * u_m1 * u_m1 * u_m1
        return ( -22u_m1_7 * u * ( 16u^2 + 7u + 1 )) * n |> T
    else
        return 0.0 |> T
    end

end

"""
    bias_correction_3D(kernel::WendlandC6, density::Real, m::Real, h_inv::Real)

Corrects the density estimate for the kernel bias. See Dehnen&Aly 2012, eq. 18+19.
"""
function bias_correction_3D(kernel::WendlandC6{T}, density::Real, m::Real, h_inv::Real) where T

    n = kernel.norm_3D * h_inv^3
    wc_correction::T = T(0.0116) * ( kernel.n_neighbours * T(0.01) )^T(-2.236) * m * n

    if wc_correction < T(0.2)*density
        density -= wc_correction
    end

    return density |> T
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CSV
using DataFrames
using Dates
using Plots
using BSON
using StatsBase
using Underscores
using ColorSchemes
colors = ColorSchemes.tableau_10.colors
# using Makie
# using GLMakie
# GLMakie.activate!()
plotlyjs()

## ============================================================================
# utils

function parsetime(sraw)
    svec = split(sraw, "-")
    if length(svec) == 3
        d,m,y = svec
        sclean = join([d,m,"20$y"], "-")
        return Date(sclean, "d-u-y")
    else
        d = sraw[1:2]
        m = sraw[3:5]
        y = sraw[6:9]
        sclean = join([d,m,y], "-")
        return Date(sclean, "d-u-y")
    end
end

function readcontacts( fname )
    @_ fname |>
        CSV.read(__, DataFrame) |>
        DataFrames.select(__,
            :dataRound_glb => :batch,
            :edate => ByRow(x->parsetime(x)) => :date,
            "country_string" => :region,
            "qNew_3_2_resp" => :family,
            "qNew_3_3_resp" => :colleagues,
            "qNew_3_4_resp" => :friends,
            "qNew_3_5_resp" => :strangers,
        ) |>
    sort(__, :date)
end

datacols(df) = Symbol.(names(df, Not([:date, :region, :batch])))

function rmfirstday(data)
    dates = data.date
    @assert issorted(dates)
    firstday = first(dates)
    i = findfirst(>(firstday), dates)
    data[i:end, :]
end

function maskoutliers(data, p=99.9)
    outliers = zeros(size(data,1))
    for c in datacols(data)
        threshold = p isa Number ? percentile(data[:,c], p) : p[c]
        @show "threshold for $c: $threshold"
        outliers .+= (threshold .< data[!,c])
        outliers .+= (data[!,c] .< 0)
    end
    return outliers .== 0
end

function filteroutliers(data, p=99.9)
    mask = maskoutliers(data, p)
    @show "remove $(sum(mask .== 0)/size(data,1)*100) percent"
    data = data[mask,:]
    data = rmfirstday(data)#NOTE first day seems unusual
    return data
end

function support_index(stepsize, support_position)
    support_position == :start && return 1
    support_position == :end && return stepsize
    iseven(stepsize) && return stepsize ÷ 2
    return stepsize ÷ 2 + 1
end

function sampledown(df; stepsize=1, support_position=:middle)
    @assert issorted(df.date)
    # df  = copy(df)
    s = Day(support_index(stepsize, support_position)-1)
    day1 = first(df.date)
    df[!, :date] =  @. ( ( df.date - day1 ) ÷ stepsize ) * stepsize + s + day1
    return df
end

function continuous_time(df)
    dates = df.date
    data_cols = names(df, Not(:date))
    @assert issorted(dates)
    ts = dates[1] : Day(1) : dates[end]
    long_df = DataFrame("date" => ts, [c => Vector{Union{Float64, Missing}}(undef, length(ts)) for c in data_cols]...)
    next_date = dates[1]
    j = 1
    for (i,t) in enumerate(ts)
        if t == next_date
            long_df[i, data_cols] = df[j, data_cols]
            j += 1
            j > size(df, 1) && return long_df
            next_date = df.date[j]
        end
    end
end

function rollingmean(xs::AbstractVector, window)
    n  = length(xs)-window+1
    ys = Vector{Float64}(undef, n)
    for i in 1:n
        x = xs[i:i+window-1]
        ys[i] = mean(skipmissing(x))
        @assert isfinite(ys[i]) "infinite at i=$i"
    end
    return ys
end

function rollingmean(df, window)
    @assert isodd(window)
    long_df = continuous_time(df)
    Δ = window ÷ 2
    smoothed = DataFrame( :date => long_df.date[1+Δ:end-Δ] )
    for col in names(df, Not(:date))
        @info col
        smoothed[!,col] = rollingmean(long_df[:,col], window)
        @assert all( isfinite.(smoothed[!,col]) ) "$(findall(x-> !isfinite(x), smoothed[!,col]))"
    end
    smoothed
end

function preprocess(df, param, region)
    @_ df |>
        filter(:region => ==(region), __) |>
        filteroutliers(__, param.threshold) |>
        groupby(__, :batch, sort = true) |>
        combine(__, param.ctypes .=> sum, :date => first, nrow, renamecols = false) |>
        select(__, Not(:batch), AsTable(param.ctypes) => ByRow(sum) => :total) |>
        transform(__, [param.ctypes..., :total] .=> ( x -> x ./ __.nrow ), renamecols = false)
end

function plotweeks!(p, dates)
    dates = @_ DataFrame( :date => dates[1] : Day(1) : dates[end] ) |>
        select(__, :date,
            :date =>( x-> week.(x) )=> :week,
            :date =>( x-> year.(x) )=> :year ) |>
        groupby(__, [:year, :week]) |>
        combine(__, :date .=> [first, last])|>
        filter(:week => x->x%2 == 0, __)

    for row in eachrow(dates)
        vspan!(p, [row.date_first, row.date_last], c=:darkblue, α=0.1, label=false)
    end
    p
end

function firefox(p; fname="tmp.html", fdir=normpath(homedir(), ".tmp"))
    # mkpath(dir)
    filepath = normpath(fdir, fname)
    savefig(p, filepath)
    run(`firefox $(filepath)`, wait=false)
    return p
end
## ===========================================================================
# main

param = (
    fname = projectdir("data/contacts/int/international_contacts.csv")
    , ctypes = [:family, :colleagues, :friends, :strangers]
    , threshold = (family=50, colleagues=100, friends=100, strangers=1000)
)

function main(param)
    rawcontacts = readcontacts( param.fname )
    regions  = unique(rawcontacts.region)
    contacts = Dict{String,DataFrame}()
    for region in regions
        regional = preprocess(rawcontacts, param, region)
        contacts[region] = regional
        CSV.write(projectdir("data/contacts/int/mean_contact_rates_region=$(region)_aggregation=batch.csv"), regional)
    end
    bson(projectdir("data/contacts/int/mean_contact_rates_aggregation=batch.bson"), contacts)
    return contacts
end
contacts = main(param)

## plot

let
    for (region, data) in pairs(contacts)
        p = scatter(data.date, data.total, title=region)
        plotweeks!(p, data.date) |> display
    end
end

# plot raw data
let
    rawcontacts = readcontacts( param.fname )
    region = "Denmark"
    rawregional = filter(:region => ==(region), rawcontacts)
    scatter(rawregional.date, rawregional.strangers, mα=0.1, mc=colors[2], msc=:match)
end

for (region, data) in pairs(contacts)
    Δt = unique(diff(data.date))
    @show region, Δt
end

## ===========================================================================
#              Smoothed Contacts
# ============================================================================
using Plots
using AbstractGPs
using Distributions
using ParameterHandling
using ParameterHandling: value, flatten
using Optim
using Zygote

## ===========================================================================
# utils

struct Loss{X,Y,F,X0}
    x::X
    y::Y
    unpack::F
    x0::X0
end

function Loss(X,Y,θ::NamedTuple)
    x0, unflatten = flatten(θ);
    unpack = value ∘ unflatten
    return Loss(X,Y,unpack,x0)
end

function (nlml::Loss)(ps::AbstractVector)
    @unpack x, y, unpack = nlml
    θ = unpack(ps)
    f = build_model(θ)
    return -logpdf(f(x, θ.s), y)
end

function plot_smoothed(x,y,ŷ,sŷ; title = "")
    p = plot(; title)
    scatter!(x,y; label = false)
    plotweeks!(p, x)
    tf = x[1]:Day(1):x[end]
    plot!(tf, ŷ; label=false, ribbon=3sŷ, hover=tf)
end

function smooth(X::AbstractVector, Y::AbstractVector)
    θ = initθ()
    loss = Loss(X,Y,θ)

    results = Optim.optimize(
        loss,
        ps -> gradient(loss, ps)[1],
        loss.x0,
        LBFGS(),
        Optim.Options(
            allow_f_increases=true,
            show_trace=true
        );
        inplace=false
    )
    @show results.ls_success
    ps = loss.unpack(results.minimizer)
    @show ps
    f = build_model(ps)
    fp = AbstractGPs.posterior(f(X, ps.s), Y)
    marg = marginals(fp(0:X[end]))
    (;
        mean = mean.(marg),
        std  = std.(marg)
    )
end

function smooth(df::DataFrame, param::NamedTuple; title = "")
    dates = df.date
    X = [d.value for d in dates .- dates[1]]

    plots = Vector{Plots.Plot}()
    smoothed = Dict{String, AbstractVector}()
    smoothed["date"] = dates[1]:Day(1):dates[end]

    for ctype in names(df, Not([:date, :nrow]))
        @show ctype
        m, s = mean_and_std(df[!,ctype])
        Y = ( df[:,ctype] .- m ) / s
        ŷ, sŷ = smooth(X,Y)

        smoothed[ctype] = ŷ * s .+ m
        smoothed["$(ctype)_std"] = sŷ * s
        p = plot_smoothed(dates, Y, ŷ, sŷ; title=ctype)
        push!(plots, p)
    end

    p = plot(plots..., layout=(5,1), size=(1000, 5*250), sharex=true, link=:x)
    firefox(p, fname="$title.html")
    return DataFrame(smoothed...)
end

function initθ()
    (
        l1 = positive(90.),
        l2 = bounded(30., 5., 50.),
        s1 = positive(1.5),
        s2 = positive(.7),
        α = positive(2.),
        s = bounded(1., 0.1, 2.),
    )
end

function build_model(θ::NamedTuple)
    GP(
        θ.s1 * SqExponentialKernel() ∘ ScaleTransform(1/θ.l1) +
        θ.s2 * SqExponentialKernel() ∘ ScaleTransform(1/14)
        # θ.s2 * RationalQuadraticKernel(α = θ.α) ∘ ScaleTransform(1/θ.l2)
    )
end

## ===========================================================================
function smooth(param::NamedTuple)
    rawcontacts = readcontacts( param.fname )
    regions  = String.( unique(rawcontacts.region) )
    contacts = Dict{String, DataFrame}()
    for region in regions
        @info region
        aggregated_contacts = preprocess(rawcontacts, param, region)[2:end,:]
        smoothed = smooth(aggregated_contacts, param, title=region)
        contacts[region] = smoothed
        CSV.write(projectdir("data/contacts/int/smoothed_contact_rates_region=$(region).csv"), smoothed)
    end
    bson(projectdir("data/contacts/int/smoothed_contact_rates.bson"), contacts)
    return contacts
end

smoothed_contacts = smooth(param)

## ==========================================================================
# replace Denmark with high resolution data
let
    den_lofi = @_ CSV.read(projectdir("data/contacts/int/smoothed_contact_rates_region=Denmark.csv"), DataFrame) |>
        rename(__, Dict(
            col => "$(col)_low_resolution" for col in filter(x->eltype( __[!,x] ) <: Number, names(__))
        ))
    den_hifi = CSV.read(projectdir("data/contacts/dk","smoothed_contact_rates_window=7.csv"), DataFrame)
    den = leftjoin(den_hifi, den_lofi; on=:date, validate=(true, true), )
    CSV.write(projectdir("data/contacts/int/smoothed_contact_rates_region=Denmark.csv"), den)
end

ext: jl

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
import BinaryProvider: detect_libgfortran_abi

"""
    detect_libgfortran_abi(oh::ObjectHandle, platform::Platform)

Given an ObjectFile, examine its dynamic linkage to discover which (if any)
`libgfortran` it's linked against.  The major SOVERSION will determine which
GCC version we're restricted to.
"""
function detect_libgfortran_abi(oh::ObjectHandle, platform::Platform)
    # We look for linkage to libgfortran
    libs = basename.(path.(DynamicLinks(oh)))
    fortran_libs = filter(l -> occursin("libgfortran", l), libs)
    if isempty(fortran_libs)
        return :gcc_any
    end

    # If we find one, pass it off to BinaryProvider.detect_libgfortran_abi()
    return detect_libgfortran_abi(first(fortran_libs), platform)
end

function check_gcc_version(oh::ObjectHandle, platform::Platform; io::IO = stdout, verbose::Bool = false)
    # First, check the GCC version to see if it is a superset of `platform`.  If it's
    # not, then we have a problem!
    gcc_version = :gcc_any

    try
        gcc_version = detect_libgfortran_abi(oh, platform)
    catch e
        if isa(e, InterruptException)
            rethrow(e)
        end
        warn(io, "$(path(oh)) could not be scanned for libgfortran dependency!")
        warn(io, e)
        return true
    end

    if verbose && gcc_version != :gcc_any
        info(io, "$(path(oh)) locks us to $(gcc_version)")
    end

    if compiler_abi(platform).gcc_version == :gcc_any && gcc_version != :gcc_any
        msg = strip(replace("""
        $(path(oh)) links to libgfortran!  This causes incompatibilities across
        major versions of GCC.  To remedy this, you must build a tarball for
        each major version of GCC.  To do this, immediately after your `platforms`
        definition in your `build_tarballs.jl` file, add the line:
        """, '\n' => ' '))
        msg *= "\n\n    platforms = expand_gcc_versions(platforms)"
        warn(io, msg)
        return false
    end
    return true
end


text:
# Bibliotecas
using Distributions
using Plots
using DataFrames
using Random
using ForwardDiff
using LinearAlgebra
using Interpolations
using DataFrames
using Optim
Random.seed!(1234)
using NBInclude

# Importando código
@nbinclude("utils.ipynb")
@nbinclude("RBC.ipynb")
@nbinclude("EGM.ipynb")
@nbinclude("SteadyState.ipynb")
@nbinclude("GenBKM.ipynb")

# Estado Estácionario
# Encontrando o valor do estado estacionario do capital
p = Params()
H_ss = 1.0
@time oo = optimize(K -> eq_K(K,p), 10, 100, Brent())
K_star = oo.minimizer;
println("O valor do capital o estado estacionario é = $(K_star)")
g_star, c_star, g_low_star, g_high_star, success_flag= solve_EGM(x->log(x), x->log(x), R(K_star, H_ss, p), W(K_star, H_ss, p), p); #solve for policy functions
t_star = make_trans_mat(g_star, p)
d_star = get_stationary_dist(t_star);

# Visualizando a convergência em direção ao estado estacionário
K_grid = collect(range(oo.minimizer-0.5, stop=oo.minimizer+0.5, length=20))
K_implied_grid = similar(K_grid)
R_grid = similar(K_grid)

for (K_index, K_value) in enumerate(K_grid)
    R_s, W_s = R(K_value, H_ss, p), W(K_value, H_ss, p)
    gg, c_star, g_low, g_high, success_flag= solve_EGM(x -> log(x), x -> 2*log(x), R_s, W_s, p)
    tt = make_trans_mat(gg, p)
    dd = get_stationary_dist(tt)
    K_implied = aggregate_K(dd, p)
    R_grid[K_index] = R_s
    K_implied_grid[K_index] = K_implied
    K_grid[K_index] = K_value
end

# Gerando os graficos da oferta e demanda do capital

plot(K_grid, R_grid, label = "Demanda", ylabel="Taxa de Juros", xlabel="Capital")
plot!(K_implied_grid, R_grid, label = "Oferta")
savefig("/home/alexandre/Documentos/Mestrado/Macroeconomia/AiyagariAggregateUncertainty")
# Choque mit
# Atualização posterior e anterior
# Função Atualização anterior
function backward_update(g_low_ss::Function, g_high_ss::Function, K_path_guess::Array{Float64,1}, H_path::Array{Float64,1}, p::Params)

    nT = length(H_path)
    g_low_path = Array{Function}(undef,nT)
    g_high_path =  Array{Function}(undef,nT)
    g_low_path[nT] = g_low_ss
    g_high_path[nT] = g_high_ss
    a_path = zeros(p.nI, p.grid_size, nT)
    R_path = zeros(nT)
    W_path = zeros(nT)

    for t=nT:-1:2
        R_path[t], W_path[t] = R(K_path_guess[t], H_path[t], p), W(K_path_guess[t], H_path[t], p)
        R_path[t-1], W_path[t-1] = R(K_path_guess[t-1], H_path[t-1], p), W(K_path_guess[t-1], H_path[t-1], p)
        a_path[:,:,t-1], c_new, g_low_path[t-1], g_high_path[t-1] = euler_back(g_low_path[t], g_high_path[t], R_path[t-1], W_path[t-1], R_path[t], W_path[t], p)
    end

    return a_path, g_low_path, g_high_path
end

# Função Atualização Posterior
function forward_update(K_star::Float64, a_path::Array{Float64,3}, d_ss::Array{Float64,1}, p::Params)
    nT = length(H_path)
    K_path_forward = zeros(nT)
    K_path_forward[1] = K_star
    dd_path_forward = zeros(size(d_ss,1), nT)
    dd_path_forward[:,1] = d_ss
    for t=2:nT
        tt = make_trans_mat(a_path[:,:,t-1], p)
        dd_path_forward[:,t] = tt*dd_path_forward[:,t-1]
        K_path_forward[t] = aggregate_K(dd_path_forward[:,t], p)
    end

    return dd_path_forward, K_path_forward
end

# Encontrando o caminho de transição
function solve_mit!(K_path, g_low_ss::Function, g_high_ss::Function, d_ss::Array{Float64,1},
                    K_ss::Float64, H_path::Array{Float64,1}, p::Params; convex_combination::Float64 = 0.2,
                    shrink_factor::Float64 = 0.5, expand_factor::Float64 = 1.05,
                    max_iter::Int64 = 1000, tol::Float64=1e-6, verbose::Bool=true, display_iter::Int64 = 20)
    diff = Inf
    diff_old = Inf
    convergence_flag = 0
    damp = convex_combination
    for i_mit=1:max_iter

        a_path, g_low_path, g_high_path = backward_update(g_low_ss, g_high_ss, K_path[i_mit], H_path, p);

        dd_path_forward, K_path_forward = forward_update(K_ss, a_path, d_star, p);

        diff = maximum(abs.(K_path_forward - K_path[i_mit]))

        if verbose==true
            if mod(i_mit,display_iter) == 0
                println("Interação $(i_mit). diferença = $(diff)")
            end
        end

        if diff < tol
            if verbose==true
                println("Convergence rejeitada após $(i_mit) Interações.")
            end
            convergence_flag = 1
            break
        else

            if diff > diff_old
                damp = max(min(damp * shrink_factor, 1.0-eps()), eps())
            else
                damp = max(min(damp * expand_factor, 1.0-eps()), eps())
            end
            if mod(i_mit, 10) == 0
                if verbose==true
                    println("damp = $(damp); diff = $(diff)")
                end
            end
            push!(K_path, damp.*K_path_forward .+ (1.0 - damp).*K_path[i_mit])
            diff_old = diff

        end
    end

    return K_path, convergence_flag
end

# Encontrando o caminho para {K_t} para um choque de 1 padrão
max_t = 300
H_ss = 1.0
H_shock = 2.0

H_path = ones(max_t)
H_path[1] = H_ss*H_shock

for t_index=2:max_t
    H_path[t_index] = H_path[t_index-1]^p.rho
end

K_path = []
push!(K_path, repeat([K_star], max_t))

@time K_path, convergence_flag = solve_mit!(K_path, g_low_star, g_high_star, d_star, K_star, H_path, p, convex_combination=0.25)

R_path = zeros(length(H_path))
W_path = zeros(length(H_path))
for t=length(H_path):-1:1
    R_path[t], W_path[t] = R(K_path[end][t], H_path[t], p), W(K_path[end][t], H_path[t], p)
end
# Visualizando a convergência do caminho de transição

p0 = plot(1:max_t, K_path[1], label= "Interação 0", title="Convergencia de K(t)")
plot!(p0, 2:max_t, K_path[2][2:end], label = "Interação 1")
show_every = 5
for k in 2:length(K_path)
    if mod(k,show_every) == 0
        plot!(p0, 2:max_t, K_path[k][2:end], xlabel="t", label = "Interação $(k)", title="Convergencia de K(t)", legend=:best)
    end
end

p0
savefig("/home/alexandre/Documentos/Mestrado/Macroeconomia/AiyagariAggregateUncertainty/2.png")

p1 = plot(1:max_t, K_path[end], label= "K_t Aiyagari", title="IRF Aiyagari versus RBC")
plot!(p1, xx[RBCp.iK,2:end] .+ K_star, label = "K_t RBC", color = "black", xlabel="t")
savefig("/home/alexandre/Documentos/Mestrado/Macroeconomia/AiyagariAggregateUncertainty/3.png")

p1 = plot(1:max_t, H_path./H_path[end] .-1, label = "H(t)", xlabel= "t")
p2 = plot(1:max_t, R_path./R_path[end] .-1, label= "R(t)", xlabel= "t")
p3 = plot(1:max_t, W_path./W_path[end].-1 , label= "W(t)", xlabel= "t")
p4 = plot(1:max_t, K_path[end]./K_path[end][end] .-1, label= "K(t)", xlabel= "t" )

p5 = plot(p1, p2, p3, p4)
savefig("/home/alexandre/Documentos/Mestrado/Macroeconomia/AiyagariAggregateUncertainty/5.png")
# checando a linearidade do modelo
max_t = 300
H_ss = 1.0

array_sigma = collect(range(-0.75, stop=0.75, step=0.25))

array_sigma = array_sigma[array_sigma .!= 0.]

x_mit_scaled_sigma = zeros(max_t, length(array_sigma))

H_path_sigma = zeros(max_t, length(array_sigma))
H_path_sigma_dev = zeros(max_t, length(array_sigma))

for (index_sigma, sigma) in enumerate(array_sigma)

    H_path = ones(max_t)
    H_path[1] = H_ss + H_ss*sigma

    for t_index=2:max_t
        H_path[t_index] = H_path[t_index-1]^p.rho
    end

    K_path = []
    push!(K_path, repeat([K_star], max_t))

    @time K_path, convergence_flag = solve_mit!(K_path, g_low_star, g_high_star, d_star, K_star, H_path, p, convex_combination=0.2, verbose=false);

    if convergence_flag!=1
        error("No convergence for H(1) = $(H_path[1]).")
    end

    H_path_sigma[:, index_sigma] = H_path

    H_path_sigma_dev[:, index_sigma] = H_path./H_ss .- 1.0

    # Scaled IRF: how a percentage deviation in z_t from its steady-state results in a % deviation of k_t
    x_mit_scaled_sigma[:, index_sigma] = (K_path[end]./K_star .- 1.0)./H_path_sigma_dev[1, index_sigma]

end

p0 = plot()
p1 = plot()
p2 = plot()
p3 = plot()

for (index_sigma, sigma) in enumerate(array_sigma)
    if index_sigma == 1
        p0 = plot(100 .*H_path_sigma_dev[:, index_sigma], label="H(t) com H(1) = $(round(100 .*H_path_sigma_dev[1, index_sigma], digits=2))%")
        p1 = plot(H_path_sigma[:, index_sigma], label="H(t) H(1) = $(round(100 .*H_path_sigma_dev[1, index_sigma], digits=2))%")
        p2 = plot(x_mit_scaled_sigma[:, index_sigma], label="H(1) = $(round(100 .*H_path_sigma_dev[1, index_sigma], digits=2))%")
        p3 = plot(sign(H_path_sigma[1, index_sigma] - H_ss)*x_mit_scaled_sigma[:, index_sigma], label="H(1) = $(round(100 .*H_path_sigma_dev[1, index_sigma], digits=2))%")
    else
        plot!(p0, 100 .*H_path_sigma_dev[:, index_sigma], label="H(t) with H(1) = $(round(100 .*H_path_sigma_dev[1, index_sigma], digits=2))%", title = "% Desvio Padrão do Capital Humano", xlabel="t")
        plot!(p1, H_path_sigma[:, index_sigma], label="H(t) with H(1) = $(round(100 .*H_path_sigma_dev[1, index_sigma], digits=2))%", title = "Agragação do Capital Humano", xlabel="t")
        plot!(p2, x_mit_scaled_sigma[:, index_sigma], label="H(1) = $(round(100 .*H_path_sigma_dev[1, index_sigma], digits=2))%", title = "Impact of shock's size on scaled IRF: x", xlabel="t")
        plot!(p3, sign(H_path_sigma[1, index_sigma] - H_ss)*x_mit_scaled_sigma[:, index_sigma], label="H(1) = $(round(100 .*H_path_sigma_dev[1, index_sigma], digits=2))%", title = "Impact of shock's sign on scaled IRF: sign(x - x_ss)*x", xlabel="t")
    end
end

plot(p2,p3, fg_legend = :transparent, legend=:best, layout=(2,1))

# Dinamica fora do estado estacionario
max_t = 2000
shocks_t = rand(Normal(0,0.005), max_t) # Series of aggregate shocks
# Let's generate a path for the aggregate shock
H_path = ones(max_t)
H_path[1] = H_ss

# Evolution of aggregate productivity in level:
for t_index=2:max_t
    H_path[t_index] = H_path[t_index-1]^p.rho + shocks_t[t_index]
end

# Calculation of GenBKM path:
XT_GenBKM = zeros(max_t);# Initialization
@time GenBKM_path!(XT_GenBKM, max_t, x_mit_scaled_sigma, H_path./H_ss .- 1.0, array_sigma)

p1 = plot(100 .*(H_path./H_ss .- 1.0), label="H_t", xlabel="t")
p2 = plot(100 .*XT_GenBKM, label = "K_t", xlabel="t")
plot(p1,p2, fg_legend = :transparent, legend=:best, layout=(2,1))


text:
include("basic_setup.jl")

using Statistics

function cost(num_responses, integration_nodes, chain_length)
    algorithm = TIEstimate(1024, integration_nodes, chain_length)
    result = Trajectories.marginal_entropy(gen, algorithm=algorithm; num_responses=num_responses, duration=100.0)
    estimator_var = var(result["marginal_entropy"].Sample) / (num_responses - 1)
    time = sum(result["marginal_entropy"].TimeElapsed)
    estimator_var * time
end

using Hyperopt

ho = @hyperopt for i in 50,
            sampler in GPSampler(Min),
            integration_nodes in 4.0:8.0,
            chain_length in 2.0.^(9:0.125:13)
    @show integration_nodes chain_length
    @show cost(200, round(Int, integration_nodes), round(Int, chain_length))
end

using Plots
p = plot(ho, dpi=300)
savefig(p, projectdir("plots", "hyperopt.png"))

using JSON

filename = projectdir("data", "hyperopt.json")

open(filename, "w") do io
    JSON.print(io, Dict(
    "iterations" => ho.iterations,
    "params" => ho.params,
    "candidates" => ho.candidates,
    "history" => ho.history,
    "results" => ho.results
    ))
end


text:
module LTK_universe
    using LinearAlgebra
    using SparseArrays
    using Random
    #export Universe

    import ..Gaugefields:GaugeFields,GaugeFields_1d,
                        IdentityGauges,RandomGauges,
                        set_wing!,
                        substitute!,
                        calc_GaugeAction,
                        SU3GaugeFields,SU3GaugeFields_1d,
                        SU2GaugeFields,SU2GaugeFields_1d,
                        SUNGaugeFields,SUNGaugeFields_1d,
                        Oneinstanton,
                        evaluate_wilson_loops!,
                        U1GaugeFields,U1GaugeFields_1d,
                        apply_smearing,calc_smearingU
    import ..Gaugefields

    import ..Fermionfields:FermionFields,WilsonFermion,StaggeredFermion,substitute_fermion!,gauss_distribution_fermi!,set_wing_fermi!
    import ..Fermionfields
    import ..Actions:GaugeActionParam,FermiActionParam,
                Setup_Gauge_action,Setup_Fermi_action,
                GaugeActionParam_standard,FermiActionParam_Wilson,
                show_parameters_action,
                FermiActionParam_WilsonClover,
                FermiActionParam_Staggered,
                GaugeActionParam_autogenerator,SmearingParam_multi,SmearingParam_single,Nosmearing,
                FermiActionParam_Domainwall
    import ..LieAlgebrafields:LieAlgebraFields,clear!,add_gaugeforce!
    import ..LieAlgebrafields
    import ..Rand:Random_LCGs
    import ..System_parameters:Params
    import ..Diracoperators:DdagD_operator,Domainwall_operator,make_densematrix,
                D5DW_Domainwall_operator
    import ..Diracoperators
    import ..Wilsonloops:make_loopforactions,Wilson_loop_set,make_originalactions_fromloops,
                make_cloverloops
    import ..Verbose_print:Verbose_level,Verbose_3,Verbose_2,Verbose_1
    import ..IOmodule:loadU
    import ..SUN_generator:Generator

    import ..RationalApprox:calc_exactvalue,calc_Anϕ
    import ..ILDG_format:ILDG,load_gaugefield,load_gaugefield!,save_binarydata
    import ..Bridge_format:load_BridgeText!
    import ..Othermethods:tdlogdet
    import ..CGmethods:bicg,cg,shiftedcg
    import ..Rhmc:get_order,get_β,get_α,get_α0,get_β_inverse,get_α_inverse,get_α0_inverse
    #import ..Smearing:apply_smearing

    """
    Your universe is described in this type.
    """
    struct Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}
        NX::Int64
        NY::Int64
        NZ::Int64
        NT::Int64
        NV::Int64
        NC::Int64
        Nwing::Int64
        Dirac_operator::Union{Nothing,String}
        U::Array{Gauge,1}
        Uold::Array{Gauge,1}
        p::Array{Lie,1}
        φ::Union{Nothing,Fermi}
        ξ::Union{Nothing,Fermi}
        η::Union{Nothing,Fermi}
        gparam::GaugeP
        fparam::Union{Nothing,FermiP}
        BoundaryCondition::Array{Int8,1}
        initial::String
        quench::Bool
        NDFALG::Int64

        _temporal_gauge::Array{Gauge_temp,1}
        _temporal_fermi::Union{Nothing,Array{Fermi,1}}
        _temporal_algebra::Array{Lie,1}
        ranf::Random_LCGs
        verboselevel::Int8
        kind_of_verboselevel::Verbose_level
        #rng::MersenneTwister

    end

    function set_β!(univ::Universe,β)
        if typeof(univ.gparam) == GaugeActionParam_standard
            univ.gparam.β = β
            #univ.gparam = GaugeActionParam_standard(β,univ.gparam.NTRACE)
        elseif typeof(univ.gparam) == GaugeActionParam_autogenerator
            #univ.gparam =GaugeActionParam_autogenerator(univ.gparam.βs * (univ.gparam.β/β),β,univ.gparam.numactions,univ.gparam.NTRACE,univ.gparam.loops,univ.gparam.staples)
            univ.gparam.β = β
            univ.gparam.βs .= 0#univ.gparam.βs * (univ.gparam.β/β)
            univ.gparam.βs[1] = β
        else
            error("$(typeof(univ.gparam)) is not supported!")
        end
    end

    function set_β!(univ::Universe,βs::Array{T,1}) where T <: Number
        if typeof(univ.gparam) == GaugeActionParam_standard
            univ.gparam.β = βs[1]
            #univ.gparam = GaugeActionParam_standard(β,univ.gparam.NTRACE)
        elseif typeof(univ.gparam) == GaugeActionParam_autogenerator
            @assert length(univ.gparam.βs) == length(βs) "univ.gparam.βs = $(univ.gparam.βs) and βs = $βs"
            #univ.gparam =GaugeActionParam_autogenerator(univ.gparam.βs * (univ.gparam.β/β),β,univ.gparam.numactions,univ.gparam.NTRACE,univ.gparam.loops,univ.gparam.staples)
            univ.gparam.β = βs[1]
            univ.gparam.βs .= βs
            #univ.gparam.βs = univ.gparam.βs * (univ.gparam.β/β)
        else
            error("$(typeof(univ.gparam)) is not supported!")
        end
    end

    function get_β(univ::Universe)
        if typeof(univ.gparam) == GaugeActionParam_standard
            return univ.gparam.β
        elseif typeof(univ.gparam) == GaugeActionParam_autogenerator
            return univ.gparam.β
        else
            error("$(typeof(univ.gparam)) is not supported!")
        end
    end

    function set_βs!(univ::Universe,βs)
        if typeof(univ.gparam) == GaugeActionParam_standard
            univ.gparam.β = βs[1]
        elseif typeof(univ.gparam) == GaugeActionParam_autogenerator
            univ.gparam.β = βs[1]
            univ.gparam.βs = βs[:]
        else
            error("$(typeof(univ.gparam)) is not supported!")
        end
    end

    include("default.jl")

    function Universe()
        file = "default.jl"
        return Universe(file)
    end

    function Universe(p::Params)
        L = p.L
        β = p.β
        NTRACE = p.NC
        #rng = MersenneTwister(p.randomseed, (0, 1002, 0, 1))

        if p.use_autogeneratedstaples
            #loops = make_loopforactions(p.couplinglist,L)

            #println(loops)
            loops = make_originalactions_fromloops(p.coupling_loops)

            gparam = GaugeActionParam_autogenerator(p.couplingcoeff,loops,p.NC,p.couplinglist)
        else
            gparam =  GaugeActionParam_standard(β,NTRACE)
        end

        if p.Dirac_operator == nothing
            fparam = nothing
        else
            if p.Dirac_operator == "Wilson"
                if p.smearing_for_fermion == "nothing"
                    fparam = FermiActionParam_Wilson(p.hop,p.r,p.eps,p.Dirac_operator,p.MaxCGstep,p.quench)
                else
                    fparam = FermiActionParam_Wilson(p.hop,p.r,p.eps,p.Dirac_operator,p.MaxCGstep,p.quench,
                                                        smearingparameters = "stout",
                                                        loops_list = p.stout_loops,
                                                        coefficients  = p.stout_ρ,
                                                        numlayers = p.stout_numlayers,
                                                        L = p.L)
                end
            elseif p.Dirac_operator == "WilsonClover"
                #if p.NC == 2
                #    error("You use NC = 2. But WilsonClover Fermion is not supported in SU2 gauge theory yet. Use NC = 3.")
                #end
                NV = prod(p.L)
                #CloverFμν = zeros(ComplexF64,p.NC,p.NC,NV,6)
                inn_table= zeros(Int64,NV,4,2)
                internal_flags = zeros(Bool,2)
                _ftmp_vectors = Array{Array{ComplexF64,3},1}(undef,6)
                for i=1:6
                    _ftmp_vectors[i] = zeros(ComplexF64,p.NC,NV,4)
                end

                _is1 = zeros(Int64,NV)
                _is2 = zeros(Int64,NV)

                if p.NC ≥ 4
                    SUNgenerator = Generator(p.NC)
                else
                    SUNgenerator = nothing
                end

                _cloverloops = Array{Wilson_loop_set,2}(undef,3,4)
                for μ=1:3
                    for ν=μ+1:4
                        _cloverloops[μ,ν] = make_cloverloops(μ,ν)
                    end
                end

                #fparam = FermiActionParam_WilsonClover(p.hop,p.r,p.eps,p.Dirac_operator,p.MaxCGstep,p.Clover_coefficient,CloverFμν,
                #                internal_flags,inn_table,_ftmp_vectors,_is1,_is2,
                #                p.quench)
                if p.smearing_for_fermion == "nothing"
                    fparam = FermiActionParam_WilsonClover(p.hop,p.r,p.eps,p.Dirac_operator,p.MaxCGstep,p.Clover_coefficient,
                                    internal_flags,inn_table,_ftmp_vectors,_is1,_is2,
                                    p.quench,SUNgenerator,_cloverloops)
                else
                    error("stout for Wilson with clover is not supported yet!")
                    fparam = FermiActionParam_WilsonClover(p.hop,p.r,p.eps,p.Dirac_operator,p.MaxCGstep,p.Clover_coefficient,
                                    internal_flags,inn_table,_ftmp_vectors,_is1,_is2,
                                    p.quench,SUNgenerator,_cloverloops,
                                    smearingparameters = "stout",
                                    loops_list = p.stout_loops,
                                    coefficients  = p.stout_ρ,
                                    numlayers = p.stout_numlayers,
                                    L = p.L)
                end
            elseif p.Dirac_operator == "Staggered"
                if p.smearing_for_fermion == "nothing"
                    fparam = FermiActionParam_Staggered(p.mass,p.eps,p.Dirac_operator,p.MaxCGstep,p.quench,p.Nf)
                else
                    fparam = FermiActionParam_Staggered(p.mass,p.eps,p.Dirac_operator,p.MaxCGstep,p.quench,p.Nf,
                                                            smearingparameters = "stout",
                                                            loops_list = p.stout_loops,
                                                            coefficients  = p.stout_ρ,
                                                            numlayers = p.stout_numlayers,
                                                            L = p.L
                                                        )
                end
            elseif p.Dirac_operator == "Domainwall"
                if p.smearing_for_fermion == "nothing"
                    fparam = FermiActionParam_Domainwall(p.Domainwall_N5,p.r,p.Domainwall_M,p.Domainwall_m,p.Domainwall_ωs,
                    p.Domainwall_b,p.Domainwall_c,
                    p.eps,p.Dirac_operator,
                    p.MaxCGstep,p.quench)
                else
                    fparam = FermiActionParam_Domainwall(p.Domainwall_N5,p.r,p.Domainwall_M,p.Domainwall_m,p.Domainwall_ωs,
                                                            p.Domainwall_b,p.Domainwall_c,
                                                            p.eps,p.Dirac_operator,
                                                            p.MaxCGstep,p.quench,
                                                            smearingparameters = "stout",
                                                            loops_list = p.stout_loops,
                                                            coefficients  = p.stout_ρ,
                                                            numlayers = p.stout_numlayers,
                                                            L = p.L
                                                        )
                end
                #error("Not implemented! here is in LTK_universe.jl")
            else
                error(p.Dirac_operator," is not supported!")
            end
        end

        univ = Universe(L,gparam,p.Nwing,fparam,p.BoundaryCondition,p.initial,p.NC,p.verboselevel,p.load_fp,p.loadU_format,p.julian_random_number)

    end

    """
    ```Universe(file)```
    - file: file name of the input file.

    Make your universe. The input file is loaded.

    Undefined parameters in your input file are defined with the default values.

    The default values are as follows.

    ```julia
        L = (4,4,4,4)
        β = 6
        NTRACE = 3
        #gparam = Setup_Gauge_action(β)
        gparam =  GaugeActionParam_standard(β,NTRACE)

        BoundaryCondition=[1,1,1,-1]
        Nwing = 1
        initial="cold"
        NC =3

        hop= 0.141139 #Hopping parameter
        r= 1 #Wilson term
        eps= 1e-19
        Dirac_operator= "Wilson"
        MaxCGstep= 3000

        fparam = FermiActionParam_Wilson(hop,r,eps,Dirac_operator,MaxCGstep)

    ```
    """
    function Universe(file)
        if file != "default.jl"
            include(pwd()*"/"*file)
        end
        univ = Universe(L,gparam,Nwing,fparam,BoundaryCondition,initial,NC)

    end

    function show_parameters(univ::Universe)
        println("""
            L = $((univ.NX,univ.NY,univ.NZ,univ.NT))
            """)

        show_parameters_action(univ.gparam)

        println("""
            BoundaryCondition= $(univ.BoundaryCondition)
            Nwing = $(univ.Nwing)
            initial= $(univ.initial)
            NC =$(univ.NC)
            """)

        if univ.fparam != nothing
            show_parameters_action(univ.fparam)
        end
    end

    """
    ```Universe(L::Tuple,gparam::GaugeActionParam,initial="cold",fparam=nothing)```

    - L: system size (NX,NY,NZ,NT)
    - gparam: parameters for gauge actions
    - [initial]: initial Gauge configuration
    - [fparam]: parameters for fermion actions
    """
    function Universe(L::Tuple,gparam::GaugeActionParam,initial="cold",fparam=nothing)
        Nwing = 1
        BoundaryCondition=[1,1,1,-1]
        NC =3
        Universe(L,gparam,Nwing,fparam,BoundaryCondition,initial,NC)
    end

    function Universe(L,gparam,Nwing,fparam,BoundaryCondition,initial,NC,verboselevel,load_fp,loadU_format,julian_random_number)
        #(L::Tuple,gparam::GaugeActionParam;
        #   Nwing = 1,fparam=nothing,
        #  BoundaryCondition=[1,1,1,-1],initial="cold",NC =3)

        if verboselevel == 1
            kind_of_verboselevel = Verbose_1(load_fp)
        elseif verboselevel == 2
            kind_of_verboselevel = Verbose_2(load_fp)
        elseif verboselevel == 3
            kind_of_verboselevel = Verbose_3(load_fp)
        end

        ranf = Random_LCGs(1,julian_random_number=julian_random_number)

        NX = L[1]
        NY = L[2]
        NZ = L[3]
        NT = L[4]
        NV = NX*NY*NZ*NT
        NDFALG  = ifelse(NC==1,1,NC^2-1)

        num_tempfield_g = 5
        num_tempfield_f = 4

        if fparam != nothing && fparam.Dirac_operator == "WilsonClover"
            numbasis = NC^2-1
            num_tempfield_g += 4+4+2numbasis
        elseif fparam != nothing # && fparam.Dirac_operator == "Staggered"
            num_tempfield_f += 4
        end

        if fparam != nothing && fparam.Dirac_operator == "Staggered"
            if fparam.Nf != 4 && fparam.Nf != 8
                N_action = get_order(fparam.rhmc_action)
                N_MD = get_order(fparam.rhmc_MD)
                num_tempfield_f += maximum((N_action,N_MD))+1
            end
        end

        if fparam != nothing && fparam.smearing != nothing
            num_tempfield_g += 4
        end

        if typeof(gparam) == GaugeActionParam_autogenerator
            num_tempfield_g += 1
        end

        if NC == 3
            U = Array{SU3GaugeFields,1}(undef,4)
            _temporal_gauge = Array{SU3GaugeFields_1d,1}(undef,num_tempfield_g)
        elseif NC == 2
            U = Array{SU2GaugeFields,1}(undef,4)
            _temporal_gauge = Array{SU2GaugeFields_1d,1}(undef,num_tempfield_g)
        elseif NC ≥ 4
            U = Array{SUNGaugeFields{NC},1}(undef,4)
            _temporal_gauge = Array{SUNGaugeFields_1d{NC},1}(undef,num_tempfield_g)
        elseif NC == 1
            U = Array{U1GaugeFields,1}(undef,4)
            _temporal_gauge = Array{U1GaugeFields_1d,1}(undef,num_tempfield_g)
        end

#        Uold = Array{GaugeFields,1}(undef,4)

        p = Array{LieAlgebraFields,1}(undef,4)
        for μ=1:4
            p[μ] = LieAlgebraFields(NC,NX,NY,NZ,NT)
        end

        if initial == "cold"
            println(".....  Cold start")
            for μ=1:4
                U[μ] = IdentityGauges(NC,NX,NY,NZ,NT,Nwing)
            end
        elseif initial == "hot"
            println(".....  Hot start")
            for μ=1:4
                U[μ] = RandomGauges(NC,NX,NY,NZ,NT,Nwing)
            end
        elseif initial == "one instanton"
            @assert NC == 2 "From one instanton start, NC should be 2!"
            U = Oneinstanton(NC,NX,NY,NZ,NT,Nwing)

        else #if initial == "file"
            println(".....  File start")
            println("File name is $initial")
            if loadU_format == "JLD"
                U = loadU(initial,NX,NY,NZ,NT,NC)
            elseif loadU_format == "ILDG"
                ildg = ILDG(initial)
                i = 1
                for μ=1:4
                    U[μ] = IdentityGauges(NC,NX,NY,NZ,NT,Nwing)
                end
                load_gaugefield!(U,i,ildg,L,NC)
            elseif loadU_format == "BridgeText"
                for μ=1:4
                    U[μ] = IdentityGauges(NC,NX,NY,NZ,NT,Nwing)
                end
                load_BridgeText!(initial,U,L,NC)
            else
                error("loadU_format should be JLD or ILDG")
            end

            #error("not supported yet.")
        end

        set_wing!(U)
        Uold = similar(U)
        substitute!(Uold,U)

        #=
        for μ=1:4
            set_wing!(U[μ])
            Uold[μ] = similar(U[μ])
            substitute!(Uold[μ],U[μ])
        end
        =#

        #_temporal_gauge = Array{GaugeFields,1}(undef,4)
        for i=1:length(_temporal_gauge)
            _temporal_gauge[i] = GaugeFields_1d(NC,NX,NY,NZ,NT) #similar(U[1])
            #_temporal_gauge[i] = similar(U[1])
        end

        _temporal_algebra = Array{LieAlgebraFields,1}(undef,1)
        for i=1:length(_temporal_algebra)
            _temporal_algebra[i] = similar(p[1])
        end

        if fparam == nothing
            Dirac_operator = nothing
            φ = nothing
            η = nothing
            ξ = nothing
            _temporal_fermi = nothing
            quench = true
        else
            Dirac_operator = fparam.Dirac_operator
            φ = FermionFields(NC,NX,NY,NZ,NT,fparam,BoundaryCondition)
            η = similar(φ)
            ξ = similar(φ)
            _temporal_fermi = Array{FermionFields,1}(undef,num_tempfield_f )
            for i=1:length(_temporal_fermi)
                _temporal_fermi[i] = similar(φ)
            end
            quench = fparam.quench

            if Dirac_operator =="Domainwall"
                gauss_distribution_fermi!(ξ)
                #A = D5DW_Domainwall_operator(U,ξ,fparam)
                A = Domainwall_operator(U,ξ,fparam)
                A_dense_PV = make_densematrix(A.D5DW_PV)
                A_dense = make_densematrix(A.D5DW)
                M = A_dense*inv(A_dense_PV)
                #e,v = eigen(A_dense)
                e,v = eigen(M)
                fp = open("eigenvalues_M.txt","w")
                for ene in e
                    println(fp,real(ene),"\t",imag(ene))
                end
                close(fp)
                #println(e)
                error("stop!")
                bicg(η,A,ξ)
                error("stop!")
            end

        end

        Gauge =eltype(U)
        Fermi = typeof(φ)
        Lie = eltype(p)
        GaugeP = typeof(gparam)
        FermiP = typeof(fparam)
        Gauge_temp = eltype(_temporal_gauge)

        return Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}(
            NX,
            NY,
            NZ,
            NT,
            NV,
            NC,
            Nwing,
            Dirac_operator,
            U,
            Uold,
            p,
            φ,
            ξ,
            η,
            gparam,
            fparam,
            BoundaryCondition,
            initial,
            quench,
            NDFALG,
            _temporal_gauge,
            _temporal_fermi,
            _temporal_algebra,
            ranf,
            verboselevel,
            kind_of_verboselevel
        )

    end

    function gauss_distribution(univ::Universe,nv)
        variance = 1
        nvh = div(nv,2)
        granf = zeros(Float64,nv)
        for i=1:nvh
            rho = sqrt(-2*log(univ.ranf())*variance)
            theta = 2pi*univ.ranf()
            granf[i] = rho*cos(theta)
            granf[i+nvh] = rho*sin(theta)
        end
        if 2*nvh == nv
            return granf
        end

        granf[nv] = sqrt(-2*log(univ.ranf())*variance) * cos(2pi*univ.ranf())
        return granf
    end

    function expF_U!(U::Array{T,1},F::Array{N,1},Δτ,univ::Universe) where {T<: GaugeFields, N <: LieAlgebraFields}
        LieAlgebrafields.expF_U!(U,F,Δτ,univ._temporal_gauge,univ._temporal_algebra[1])
    end

    function calc_gaugeforce!(F::Array{N,1},univ::Universe) where N <: LieAlgebraFields
        clear!(F)
        calc_gaugeforce!(F,univ.U,univ)
        #add_gaugeforce!(F,univ.U,univ._temporal_gauge,univ._temporal_algebra[1])
        return
    end

    function calc_gaugeforce!(F::Array{N,1},U::Array{T,1},univ::Universe) where {N<: LieAlgebraFields, T<: GaugeFields}
        clear!(F)
        calc_gaugeforce!(F,U,univ,univ.gparam)
        #add_gaugeforce!(F,U,univ._temporal_gauge,univ._temporal_algebra[1],univ.gparam)
        return
    end

    function calc_gaugeforce!(F::Array{N,1},U::Array{T,1},univ::Universe,gparam::GP) where {N<: LieAlgebraFields,
                                                                                T<: GaugeFields,
                                                                                GP <: GaugeActionParam_autogenerator}
        clear!(F)
        add_gaugeforce!(F,U,univ._temporal_gauge,univ._temporal_algebra[1],gparam)
        return
    end

    function calc_gaugeforce!(F::Array{N,1},U::Array{T,1},univ::Universe,gparam::GP) where {N<: LieAlgebraFields,
                                                                                T<: GaugeFields,
                                                                                GP}
        clear!(F)
        add_gaugeforce!(F,U,univ._temporal_gauge,univ._temporal_algebra[1])
        return
    end

    function Gaugefields.calc_GaugeAction(univ::Universe)
        Sg,plaq = calc_GaugeAction(univ.U,univ.gparam,univ._temporal_gauge)
        return real(Sg),real(plaq)
    end

    function calc_Action(univ::Universe)
        Sg,plaq = calc_GaugeAction(univ.U,univ.gparam,univ._temporal_gauge)

        SP = univ.p*univ.p/2
        #println("ek = $SP")
        S = Sg + SP
        #println("eym = $plaq")
        #println("e = $S")

        if univ.quench == false
        #if univ.fparam != nothing
            Sf = univ.η*univ.η
            S += Sf
        end
        return real(S),real(plaq)
    end

    function calc_IntegratedFermionAction(univ::Universe;debug = false,M=16,m=100,nc=20,nonc=false)
        println("Making W^+W matrix...")
        NX = univ.U[1].NX
        NY = univ.U[1].NY
        NZ = univ.U[1].NZ
        NT = univ.U[1].NT
        NC = univ.U[1].NC
        T = eltype(univ._temporal_fermi)
        if T == StaggeredFermion
            NG = 1
        else
            NG = 4
        end
        Nsize = NX*NY*NZ*NT*NC*NG

        WdagW = zeros(ComplexF64,Nsize,Nsize)
        @time  make_WdagWmatrix!(WdagW,univ)
        #@time WdagW = make_WdagWmatrix(univ)
        Sfnew = calc_IntegratedFermionAction(univ,WdagW,debug = debug,M=M,m=m,nc=nc,nonc=nonc)
        return Sfnew
    end

    function calc_IntegratedFermionAction(univ::Universe,WdagW;debug = false,M=16,m=100,nc=20,nonc=false)
        if debug

            #WdagW = DdagD_operator(univ.U,univ.η,univ.fparam)
            #println("Calculating approximated logdet")
            #@time Sfnew = -tdlogdet(WdagW,M,m,nc=nc,nonc=nonc)
            #exit()

            WdagWsp = sparse(WdagW)
            N,_ = size(WdagW)
            tempvecs = Array{Array{ComplexF64,1}}(undef,3)
            for i=1:3
                tempvecs[i] = zeros(ComplexF64,N)
            end
            println("Calculating approximated logdet")

            #@time Sfnew = -tdlogdet(WdagWsp,M2,m2,tempvecs,nc=nc2,nonc=nonc)

            #M2 = 100
            #m2 = 100
            #nc2 = 100
            @time Sfnew = -tdlogdet(WdagWsp,M,m,tempvecs,nc=nc,nonc=nonc)

            #exit()
            #=

            M = 8*2
            #M = 4#8*2
            m = 100
            NX = univ.U[1].NX
            #m = 20
            N,_ = size(WdagW)
            tempvecs = Array{Array{ComplexF64,1}}(undef,3)
            for i=1:3
                tempvecs[i] = zeros(ComplexF64,N)
            end

            nc = 20
            nonc = false

            @time ldet = logdet(WdagW)
            if nonc
                fp = open("deltavalue_m$(m)_nonc_NX$(NX)_6.dat","w")
            else
                fp = open("deltavalue_m$(m)_nc$(nc)_NX$(NX)_6.dat","w")
            end

            for M in [4,8,16,32,64,128]
                if nonc
                    filename="delta_$(NX)$(NX)$(NX)$(NX)_$(M)_m$(m)_nonc_6.dat"
                else
                    filename="delta_$(NX)$(NX)$(NX)$(NX)_$(M)_m$(m)_nc$(nc)_6.dat"
                end
                a = tdlogdet(WdagWsp,M,m,tempvecs,filename=filename,nc=nc,nonc=nonc)
                println("$M $a $ldet")
                println(fp,"$M $a $ldet")
            end
            close(fp)

            exit()
            @time Sfnew = -tdlogdet(WdagWsp,M,m,tempvecs)
            =#
        else
            println("Calculating logdet")
            @time Sfnew = -real(logdet(WdagW))
        end

        if univ.Dirac_operator == "Staggered"
            #if univ.fparam.Nf == 4
            #    Sfnew /= 2
            #end
            Sfnew /= (8/univ.fparam.Nf)
        end
        return Sfnew
    end

    struct Wilsonloops_actions{Gauge_temp}
        loops::Array{Wilson_loop_set,1}
        loopaction::Gauge_temp
        temps::Array{Gauge_temp,1}
        numloops::Int64
        β::Float64
        NC::Float64
        couplinglist::Array{String,1}

        function Wilsonloops_actions(univ)
            couplinglist = ["plaq","rect","polyx","polyy","polyz","polyt"]
            L = (univ.U[1].NX,univ.U[1].NY,univ.U[1].NZ,univ.U[1].NT)
            loops = make_loopforactions(couplinglist,L)
            loopaction = similar(univ._temporal_gauge[1])
            sutype = typeof(loopaction)
            temps = Array{sutype,1}(undef,3)
            for i=1:3
                temps[i] = similar(loopaction)
            end
            numloops = length(couplinglist)
            β = univ.gparam.β
            NC = univ.NC

            return new{sutype}(loops,loopaction,temps,numloops,β,NC,couplinglist)
        end
    end

    function calc_looptrvalues(w::Wilsonloops_actions,univ)
        numloops = w.numloops
        trs = zeros(ComplexF64,numloops)
        for i=1:numloops
            evaluate_wilson_loops!(w.loopaction,w.loops[i],univ.U,w.temps[1:3])
            trs[i] = tr(w.loopaction)
        end
        return trs
    end

    function calc_looptrvalues_site(w::Wilsonloops_actions,univ)
        return
        numloops = w.numloops
        NX = univ.U[1].NX
        NY = univ.U[1].NY
        NZ = univ.U[1].NZ
        NT = univ.U[1].NT
        NC = univ.NC
        V = zeros(ComplexF64,NC,NC)

        trs = zeros(ComplexF64,numloops,NX,NY,NZ,NT)
        for it = 1:NT
            for iz = 1:NZ
                for iy = 1:NY
                    for ix=1:NX
                        for i=1:numloops
                            evaluate_wilson_loops!(V,w.loops[i],univ.U,ix,iy,iz,it)
                            trs[i,ix,iy,iz,it] = tr(V)
                        end
                    end
                end
            end
        end

        return trs
    end

    function get_looptrvalues(w::Wilsonloops_actions)
        return w.trs
    end

    function calc_IntegratedSf(w::Wilsonloops_actions,univ)
        Sf = calc_IntegratedFermionAction(univ)
        return Sf
    end

    function calc_gaugeSg(w::Wilsonloops_actions,univ)
        trs = calc_looptrvalues(w,univ)
        Sg = (-w.β/w.NC)*trs[1]
        return real(Sg),trs
    end

    function calc_trainingdata(w::Wilsonloops_actions,univ)
        Sf = calc_IntegratedSf(w,univ)
        Sg,trs = calc_gaugeSg(w,univ)
        return trs,Sg,Sf
    end

    function make_WdagWmatrix(univ::Universe)
        return make_WdagWmatrix(univ.U,univ._temporal_fermi,univ.fparam)
    end

    function make_WdagWmatrix!(WdagW,univ::Universe)
        make_WdagWmatrix!(WdagW,univ.U,univ._temporal_fermi,univ.fparam)
        return
    end

    function make_WdagWmatrix(univ::Universe,U::Array{T,1}) where T <: GaugeFields
        return make_WdagWmatrix(U,univ._temporal_fermi,univ.fparam)
    end

    function make_WdagWmatrix(U::Array{G,1},temps::Array{T,1},fparam) where {G <: GaugeFields,T <:FermionFields}
        x0 = temps[7]
        xi = temps[8]
        return make_WdagWmatrix(U,x0,xi,fparam)
    end

    function make_WdagWmatrix!(WdagW,U::Array{G,1},temps::Array{T,1},fparam) where {G <: GaugeFields,T <:FermionFields}
        x0 = temps[7]
        xi = temps[8]
        make_WdagWmatrix!(WdagW,U,x0,xi,fparam)
        return
    end

    function make_WdagWmatrix(U::Array{G,1},x0::T,xi::T,fparam) where {G <: GaugeFields,T <:FermionFields}
        NX = x0.NX
        NY = x0.NY
        NZ = x0.NZ
        NT = x0.NT
        NC = x0.NC
        if T == StaggeredFermion
            NG = 1
        else
            NG = 4
        end
        Nsize = NX*NY*NZ*NT*NC*NG

        WdagW = zeros(ComplexF64,Nsize,Nsize)
        make_WdagWmatrix!(WdagW,U,x0,xi,fparam)
        return WdagW
    end

    function make_WdagWmatrix!(WdagW,U::Array{G,1},x0::T,xi::T,fparam) where {G <: GaugeFields,T <:FermionFields}

        Fermionfields.clear!(x0)
        NX = x0.NX
        NY = x0.NY
        NZ = x0.NZ
        NT = x0.NT
        NC = x0.NC
        if T == StaggeredFermion
            NG = 1
        else
            NG = 4
        end
        Nsize = NX*NY*NZ*NT*NC*NG
        #WdagW = zeros(ComplexF64,Nsize,Nsize)
        WdagWoperator = DdagD_operator(U,x0,fparam)

        j = 0
        for α=1:NG
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            for ic=1:NC
                                Fermionfields.clear!(x0)
                                j += 1
                                x0[ic,ix,iy,iz,it,α] = 1
                                #set_wing_fermi!(x0)

                                #WdagWx!(xi,U,x0,temps,fparam)
                                #println(xi*xi)

                                #Wx!(xi,U,x0,temps,fparam)
                                #println(xi*xi)

                                #Wx!(xi,U,x0,temps,fparam,(ix,iy,iz,it,α))
                                #println(xi*xi)

                                mul!(xi,WdagWoperator,x0,(ix,iy,iz,it,α))
                                #WdagWx!(xi,U,x0,temps,fparam,(ix,iy,iz,it,α))
                                #println("fast ",xi*xi)
                                #exit()
                                #display(xi)
                                #exit()
                                substitute_fermion!(WdagW,j,xi)

                            end
                        end
                    end
                end
            end
        end

        #return WdagW

    end

    function construct_fermion_gauss_distribution!(univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp})  where {Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}
        gauss_distribution_fermi!(univ.η,univ.ranf)
    end

    function construct_fermion_gauss_distribution!(univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp})  where {Gauge,Lie,Fermi <: StaggeredFermion,GaugeP,FermiP,Gauge_temp}
        gauss_distribution_fermi!(univ.η,univ.ranf)
        if univ.fparam.Nf == 4
            U,_... = calc_smearingU(univ.U,univ.fparam.smearing)

            evensite = false
            W = Diracoperators.Dirac_operator(U,univ.η,univ.fparam)
            mul!(univ.φ,W',univ.η)

            #Wdagx!(univ.φ,univ.U,univ.η,univ._temporal_fermi,univ.fparam)
            Fermionfields.clear!(univ.φ,evensite)

            bicg(univ.η,W',univ.φ,eps = univ.fparam.eps,maxsteps= univ.fparam.MaxCGstep)
        end
    end

    function construct_fermionfield_φ!(univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp})  where {Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}
        U,_... = calc_smearingU(univ.U,univ.fparam.smearing)

        W = Diracoperators.Dirac_operator(U,univ.η,univ.fparam)
        mul!(univ.φ,W',univ.η)
        set_wing_fermi!(univ.φ)
    end

    function construct_fermionfield_φ!(univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp})  where {Gauge,Lie,Fermi <: StaggeredFermion,GaugeP,FermiP,Gauge_temp}
        U,_... = calc_smearingU(univ.U,univ.fparam.smearing)

        if univ.fparam.Nf == 4 || univ.fparam.Nf == 8

            W = Diracoperators.Dirac_operator(U,univ.η,univ.fparam)
            mul!(univ.φ,W',univ.η)
        else
            WdagW = Diracoperators.DdagD_operator(U,univ.η,univ.fparam)
            N = get_order(univ.fparam.rhmc_action)

            x = univ.φ
            vec_x = univ._temporal_fermi[end-N+1:end]
            for j=1:N
                Fermionfields.clear!(vec_x[j])
            end

            # eta =  (MdagM)^{-alpha/2} phi -> phi = (MdagM)^{alpha/2} eta
            #(MdagM)^{alpha/2} eta ~ α0 eta + sum_k αk (MdagM + βk)^{-1} eta
            vec_β = get_β(univ.fparam.rhmc_action)
            vec_α = get_α(univ.fparam.rhmc_action)
            α0 = get_α0(univ.fparam.rhmc_action)
            shiftedcg(vec_x,vec_β,x,WdagW,univ.η,eps = univ.fparam.eps,maxsteps= univ.fparam.MaxCGstep)
            Fermionfields.clear!(univ.φ)
            Fermionfields.add!(univ.φ,α0,univ.η)
            for j=1:N
                αk = vec_α[j]
                Fermionfields.add!(univ.φ,αk,vec_x[j])
            end

        end

        set_wing_fermi!(univ.φ)
    end

    function construct_fermionfield_η!(univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp})  where {Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp}
        U,_... = calc_smearingU(univ.U,univ.fparam.smearing)

        W = Diracoperators.Dirac_operator(U,univ.η,univ.fparam)
        bicg(univ.η,W',univ.φ,eps = univ.fparam.eps,maxsteps= univ.fparam.MaxCGstep,verbose = univ.kind_of_verboselevel)
    end

    function construct_fermionfield_η!(univ::Universe{Gauge,Lie,Fermi,GaugeP,FermiP,Gauge_temp})  where {Gauge,Lie,Fermi <: StaggeredFermion,GaugeP,FermiP,Gauge_temp} # eta = Wdag^{-1} phi
        U,_... = calc_smearingU(univ.U,univ.fparam.smearing)

        if univ.fparam.Nf == 4 || univ.fparam.Nf == 8
            W = Diracoperators.Dirac_operator(U,univ.η,univ.fparam)
            bicg(univ.η,W',univ.φ,eps = univ.fparam.eps,maxsteps= univ.fparam.MaxCGstep,verbose = univ.kind_of_verboselevel)
        else
            WdagW = Diracoperators.DdagD_operator(U,univ.η,univ.fparam)
            N = get_order(univ.fparam.rhmc_action)

            x = univ.η
            vec_x = univ._temporal_fermi[end-N+1:end]
            for j=1:N
                Fermionfields.clear!(vec_x[j])
            end

            # eta =  (MdagM)^{-alpha/2} phi -> phi = (MdagM)^{alpha/2} eta
            #(MdagM)^{alpha/2} eta ~ α0 eta + sum_k αk (MdagM + βk)^{-1} eta
            vec_β = get_β_inverse(univ.fparam.rhmc_action)
            vec_α = get_α_inverse(univ.fparam.rhmc_action)
            α0 = get_α0_inverse(univ.fparam.rhmc_action)
            shiftedcg(vec_x,vec_β,x,WdagW,univ.φ,eps = univ.fparam.eps,maxsteps= univ.fparam.MaxCGstep)
            Fermionfields.clear!(univ.η)
            Fermionfields.add!(univ.η,α0,univ.φ)
            for j=1:N
                αk = vec_α[j]
                Fermionfields.add!(univ.η,αk,vec_x[j])
            end
        end
    end

end


text:
const JServiceFactory = @jimport loci.common.services.ServiceFactory
const JService = @jimport loci.common.services.Service
const JIFormatReader = @jimport loci.formats.IFormatReader
const JImageReader = @jimport loci.formats.ImageReader
const JOMEXMLMetadata = @jimport loci.formats.ome.OMEXMLMetadata
const JMetadataStore = @jimport loci.formats.meta.MetadataStore

const datatypes = [Int8,UInt8,Int16,UInt16,Int32,UInt32,Float32,Float64,Bool]

function create_service(classname::String)
    factory = JServiceFactory(())
    cls = classforname(classname)
    service = jcall(factory, "getInstance", JService, (JClass,), cls)
    return convert(JavaObject{Symbol(classname)}, service)
end

struct OMEXMLReader
    reader::JavaObject
    meta::JavaObject

    function OMEXMLReader()
        service = create_service("loci.formats.services.OMEXMLService")
        meta = jcall(service, "createOMEXMLMetadata", JOMEXMLMetadata, ())
        reader = convert(JIFormatReader, JImageReader(()))

        jcall(reader, "setMetadataStore", Nothing, (JMetadataStore,), meta)
        new(reader, meta)
    end
end

function OMEXMLReader(f::Function, filename::AbstractString)
    oxr = OMEXMLReader(filename)
    try
        f(oxr)
    finally
        close(oxr)
    end
end

Base.close(oxr::OMEXMLReader) = jcall(oxr.reader, "close", Nothing, ())

function OMEXMLReader(filename::AbstractString)
    oxr = OMEXMLReader()
    set_id!(oxr, filename)
    oxr
end

function local_frame(f::Function; capacity=16)
    JNI.PushLocalFrame(jint(capacity))
    try
        f()
    finally
        JNI.PopLocalFrame(C_NULL)
    end
end

function set_id!(oxr::OMEXMLReader, filename::AbstractString)
    jcall(oxr.reader, "setId", Nothing, (JString,), string(filename))
end

function openbytes(oxr::OMEXMLReader, index::Int)
    local_frame() do
        jcall(oxr.reader, "openBytes", Vector{jbyte}, (jint,), index)
    end
end

function set_series!(oxr::OMEXMLReader, index::Int)
    jcall(oxr.reader, "setSeries", Nothing, (jint,), index)
end

function getpixeltype(oxr::OMEXMLReader)
    id = jcall(oxr.reader, "getPixelType", jint, ())
    return datatypes[id+1]
end

function islittleendian(oxr::OMEXMLReader)
    Bool(jcall(oxr.reader, "isLittleEndian", jboolean, ()))
end

function isinterleaved(oxr::OMEXMLReader)
    Bool(jcall(oxr.reader, "isInterleaved", jboolean, ()))
end

function get_dimension_order(oxr)
    jcall(oxr.reader, "getDimensionOrder", JString, ())
end

function get_size(oxr, d)
    jcall(oxr.reader, "getSize$d", jint, ())
end

function get_image_count(oxr)
    jcall(oxr.reader, "getImageCount", jint, ())
end

function get_xml(oxr)
    parse_string(jcall(oxr.meta, "dumpXML", JString, ()))
end


text:
using InteractBulma, InteractBase, Sass
using Test

const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_assets = joinpath(_pkg_root, "assets")

@testset "ijulia" begin
    @test !InteractBase.isijulia()
end

@testset "maketheme" begin
    variables_file = joinpath(InteractBulma.examplefolder, "flatly", "_variables.scss")
    mytheme = InteractBulma.compile_theme(_pkg_assets, variables = variables_file)
    settheme!(mytheme)
    @test isfile(joinpath(_pkg_assets, "main.min.css"))
    @test isfile(joinpath(_pkg_assets, "main_confined.min.css"))
    libs = InteractBase.libraries(mytheme)
    @test all(isfile, libs)
end


text:
using FixedEffectModels

tests = ["fit.jl",
		 "predict.jl",
		 "partial_out.jl"
		 ]

println("Running tests:")

for test in tests
	try
		include(test)
		println("\t\033[1m\033[32mPASSED\033[0m: $(test)")
	 catch e
	 	println("\t\033[1m\033[31mFAILED\033[0m: $(test)")
	 	showerror(stdout, e, backtrace())
	 	rethrow(e)
	 end
end

using Test
@test Vcov.pinvertible(Symmetric([1.0 1.0; 1.0 1.0])) ≈ [1.0 1.0; 1.0 1.0]


text:
module MLSuiteBase

using JSON, SHA, Glob
import ScikitLearnBase: is_classifier

export @grid, gridparams, reset!, istree, isrnn, is_ranker, support_multiclass, signone, modelhash, available_memory, myparam, gendir

reset!(m) = nothing
istree(m) = false
isrnn(m) = false
support_multiclass(m) = is_classifier(m)
is_ranker(m) = false
modelhash(m) = hash(m)

gridparams(grid, combine = Iterators.product) =
    unique(vec([merge(v...) for v in combine(gridparams.(reverse(grid))...)]))

gridparams(pair::Pair) = [Dict{String, Any}(pair[1] => v) for v in pair[2]]

gridparams(grid::Vector{<:Vector}) = mapreduce(gridparams, vcat, grid)

Base.sign(x::Real, Θ) = ifelse(x < -Θ, oftype(x, -1), ifelse(x > Θ, one(x), zero(x)))

function signone(x::Real, Θ = zero(x))
    y = sign(x, Θ)
    !iszero(y) && return y
    rand([-one(x), one(x)])
end

function embed_include(file, mod)
    if isfile(file)
        lines = readlines(file)
    else
        lines = split(file, '\n')
    end
    for (n, line) in enumerate(lines)
        startswith(line, '#') && continue
        m = match(r"(?<=include\().*?(?=\))", line)
        isnothing(m) && continue
        src = mod.eval(Meta.parse(m.match))
        src = isabspath(src) ? src : joinpath(dirname(file), src)
        lines[n] = embed_include(src, mod)
    end
    join(lines, '\n')
end

macro grid(n, ex)
    params = __module__.eval(ex)
    file = string(__source__.file)
    if !isfile(file) || isinteractive() ||
        get(ENV, "GRID_ENABLE", "1") == "0"
        return params[1]
    end
    line = __source__.line
    lines = split(embed_include(file, __module__), '\n')
    root = "job_" * splitext(basename(file))[1]
    project = "#!/bin/bash\nexport JULIA_PROJECT=$(Base.load_path()[1])\n"
    for param in rand(params, n)
        dir = mkpath(joinpath(root, bytes2hex(sha1(json(param)))))
        open(joinpath(dir, "param.json"), "w") do io
            JSON.print(io, param, 4)
        end
        write(joinpath(dir, "julia.sh"), project, raw"""
        source ~/.bashrc
        cd $(dirname $0)
        [ -f done ] && exit 0
        julia main.jl 2>&1 | tee julia.out
        [ $? == 0 ] && touch done
        """)
        write(joinpath(dir, "mpi.sh"), project, raw"""
        cd $(dirname $0)
        [ -f done ] && exit 0
        mpirun julia main.jl 2>&1 | tee mpi.out
        [ $? == 0 ] && touch done
        """)
        write(joinpath(dir, "srun.sh"), project, raw"""
        [ -f done ] && exit 0
        srun julia main.jl
        [ $? == 0 ] && touch done
        """)
        var_name = match(r"(.+?)(?=\s+\=)", lines[line]).match
        lines[line] = "$var_name = $param"
        write(joinpath(dir, "main.jl"), join(lines, '\n'))
        foreach(sh -> chmod(sh, 0o775), glob("*.sh", dir))
    end
    return :(exit())
end

function myparam(params)
    if haskey(ENV, "SLURM_ARRAY_TASK_ID")
        i = parse(Int, "0" * ENV["SLURM_ARRAY_TASK_ID"]) + 1
    elseif haskey(ENV, "PBS_ARRAYID")
        i = parse(Int, "0" * ENV["PBS_ARRAYID"]) + 1
    else
        i = 1
    end
    i <= length(params) ? params[i] : exit(0)
end

function gendir(param)
    dir = mkpath(bytes2hex(sha1(json(param))))
    open(joinpath(dir, "param.json"), "w") do io
        JSON.print(io, param, 4)
    end
    return dir
end

function available_memory()
    if Sys.iswindows()
        Sys.free_memory() / 1024^3
    else
        regex = r"MemAvailable:\s+(\d+)\s+kB"
        meminfo = read("/proc/meminfo", String)
        mem = match(regex, meminfo).captures[1]
        parse(Int, mem) / 1024^2
    end
end

function JSON.lower(a)
    if nfields(a) > 0
        JSON.Writer.CompositeTypeWrapper(a)
    else
        string(a)
    end
end

end


text:
immutable NDBasis{T,M} <: Space{T}
    geo::M
    fns::Vector{Vector{Shape{T}}}
end

refspace(s::NDBasis) = NDRefSpace{scalartype(s)}()

function nedelec(surface, edges)

    T = coordtype(surface)
    num_edges = numcells(edges)

    C = connectivity(edges, surface, identity)
    rows = rowvals(C)
    vals = nonzeros(C)

    fns = Vector{Vector{Shape{T}}}(num_edges)
    for (i,edge) in enumerate(cells(edges))

        fns[i] = Vector{Shape{T}}()

        for k in nzrange(C,i)

            j = rows[k] # j is the index of a cell adjacent to edge
            s = vals[k] # s contains the oriented (signed) local index of edge[i] in cell[j]

            i == 3 && @show s
            push!(fns[i], Shape{T}(j, abs(s), sign(s)))
        end
    end

    NDBasis(surface, fns)
end


text:
using NamedDims
using NamedDims:
    names,
    unify_names,
    unify_names_longest,
    unify_names_shortest,
    dim_noerror,
    expand_dimnames,
    tuple_issubset,
    tuple_cat,
    names_are_unifiable,
    order_named_inds,
    permute_dimnames,
    remaining_dimnames_from_indexing,
    remaining_dimnames_after_dropping
using Test

@testset "dim" begin
    @testset "small case" begin
        @test dim((:x, :y), :x) == 1
        @test dim((:x, :y), :y) == 2
        @test_throws ArgumentError dim((:x, :y), :z)  # not found
    end
    @testset "large case" begin
        @test dim((:x, :y, :a, :b, :c, :d), :x) == 1
        @test dim((:x, :y, :a, :b, :c, :d), :a) == 3
        @test dim((:x, :y, :a, :b, :c, :d), :d) == 6
        @test_throws ArgumentError dim((:x, :y, :a, :b, :c, :d), :z) # not found
    end
end
@testset "allocations: dim" begin
    @test 0 == @ballocated (()->dim((:a, :b), :b))()
    @test 0 == @ballocated (()->dim_noerror((:a, :b, :c), :c))()
    if VERSION >= v"1.1"
        @test 0 == @ballocated (()->dim((:a,:b), (:a,:b)))()
    else
        @test_broken 0 == @ballocated (()->dim((:a,:b), (:a,:b)))()
    end
end

@testset "unify_names_*" begin
    @test_throws DimensionMismatch unify_names((:a,), (:a, :b,))

    @test unify_names_longest((:a,), (:a, :b,)) == (:a, :b)
    @test unify_names_longest((:a,), (:a, :_)) == (:a, :_)
    @test unify_names_longest((:a, :b), (:a, :_, :c)) == (:a, :b, :c)
    @test unify_names_longest((:a,), (:a, :b), (:a, :_, :c)) == (:a, :b, :c)

    @test_throws DimensionMismatch unify_names_longest((:a, :b, :c), (:b, :a))

    @test unify_names_shortest((:a,), (:a, :b,)) == (:a,)
    @test unify_names_shortest((:a,), (:a, :_)) == (:a,)
    @test unify_names_shortest((:a, :b), (:a, :_, :c)) == (:a, :b)
    @test unify_names_shortest((:a,), (:a, :b), (:a, :_, :c)) == (:a,)

    @test_throws DimensionMismatch unify_names_shortest((:a, :b, :c), (:b, :a))

    for unify in (unify_names, unify_names_longest, unify_names_shortest)
        @test unify((:a,), (:a,)) == (:a,)
        @test unify((:a, :b), (:a, :b)) == (:a, :b)
        @test unify((:a, :_), (:a, :b)) == (:a, :b)
        @test unify((:a, :_), (:a, :_)) == (:a, :_)

        @test unify((:a, :b, :c), (:_, :_, :_)) == (:a, :b, :c)
        @test unify((:a, :_, :c), (:_, :b, :_)) == (:a, :b, :c)
        @test unify((:_, :_, :_), (:_, :_, :_)) == (:_, :_, :_)

        @test_throws DimensionMismatch unify((:a,), (:b,))
        @test_throws DimensionMismatch unify((:a,:b), (:b, :a))
        @test_throws DimensionMismatch unify((:a, :b, :c), (:_, :_, :d))
    end

    # vararg version
    @test unify_names((:a, :_), (:a, :b,), (:_, :b)) == (:a, :b)
    @test unify_names((:a, :b,)) == (:a, :b)
    @test_throws DimensionMismatch unify_names((:a, :_), (:a, :b,), (:_, :c))
end
@testset "allocations: unify_names_*" begin
    for unify in (unify_names, unify_names_longest, unify_names_shortest)
        if VERSION >= v"1.1"
            @test 0 == @ballocated (()->$unify((:a, :b), (:a, :_)))()
        else
            @test_broken 0 == @ballocated (()->$unify((:a, :b), (:a, :_)))()
        end
    end
    if VERSION >= v"1.1"
        @test 0 == @ballocated (()->unify_names_longest((:a, :b), (:a, :_, :c)))()
        @test 0 == @ballocated (()->unify_names_shortest((:a, :b), (:a, :_, :c)))()
        @test 0 == @ballocated (()->unify_names_longest((:a,), (:a, :b), (:a, :_, :c)))()
        @test 0 == @ballocated (()->unify_names_shortest((:a,), (:a, :b), (:a, :_, :c)))()
        @test 0 == @ballocated (()->names_are_unifiable((:a, :b), (:a, :_)))()
        @test 0 == @ballocated (()->names_are_unifiable((:a, :b), (:a, :c)))()
    else
        @test_broken 0 == @ballocated (()->unify_names_longest((:a, :b), (:a, :_, :c)))()
        @test_broken 0 == @ballocated (()->unify_names_shortest((:a, :b), (:a, :_, :c)))()
        @test_broken 0 == @ballocated (()->unify_names_longest((:a,), (:a, :b), (:a, :_, :c)))()
        @test_broken 0 == @ballocated (()->unify_names_shortest((:a,), (:a, :b), (:a, :_, :c)))()
        @test_broken 0 == @ballocated (()->names_are_unifiable((:a, :b), (:a, :_)))()
        @test_broken 0 == @ballocated (()->names_are_unifiable((:a, :b), (:a, :c)))()
    end
    @test 0 == @ballocated (()->names_are_unifiable((:a, :b), (:a, :b)))()
end

@testset "expand_dimnames" begin
    @testset "single symbol" begin
        @test expand_dimnames((), :x) == (:x,)
        @test expand_dimnames((:x,), ()) == (:x,)

        @test expand_dimnames((:x, :y), :x) == (:x, :y)
        @test expand_dimnames((:x, :y), :z) == (:x, :y, :z)
    end

    @testset "multiple symbols" begin
        @test expand_dimnames((), (:x, :y)) == (:x, :y)
        @test expand_dimnames((:x, :y), (:z, :w)) == (:x, :y, :z, :w)
        @test expand_dimnames((:x, :y), (:x, :w)) == (:x, :y, :w)
    end

    @testset "numbers et al" begin
        @test expand_dimnames((:x, :y), 1) == (:x, :y)
        @test expand_dimnames((:x, :y), 3) == (:x, :y, :_)
        @test expand_dimnames((:x, :y), (1, 3)) == (:x, :y, :_)
        @test expand_dimnames((:x, :y), :) == (:x, :y)
        @test expand_dimnames((:x, :y), ()) == (:x, :y)
    end
end

@testset "allocations: expand_dims" begin
    @testset "names" begin
        @test 0 == @ballocated (()->NamedDims.expand_dimnames((:x, :y), :x))()
        @test 0 == @ballocated (()->NamedDims.expand_dimnames((:x, :y), :z))()
        @test 0 == @ballocated (()->NamedDims.expand_dimnames((:x, :y), (:x, :z)))()
    end

    @testset "numbers et al" begin
        @test 0 == @ballocated (()->NamedDims.expand_dimnames((:x, :y), 1))()
        @test 0 == @ballocated (()->NamedDims.expand_dimnames((:x, :y), 5))()
        @test 0 == @ballocated (()->NamedDims.expand_dimnames((:x, :y), (1, 5)))()
        @test 0 == @ballocated (()->NamedDims.expand_dimnames((:x, :y), :))()
        @test 0 == @ballocated (()->NamedDims.expand_dimnames((:x, :y), ()))()
    end
end

@testset "order_named_inds" begin
    @test order_named_inds(Val((:x,))) == (:,)
    @test order_named_inds(Val((:x,)); x=2) == (2,)

    @test order_named_inds(Val((:x, :y))) == (:, :)
    @test order_named_inds(Val((:x, :y)); x=2) == (2, :)
    @test order_named_inds(Val((:x, :y)); y=2, ) == (:, 2)
    @test order_named_inds(Val((:x, :y)); y=20, x=30) == (30, 20)
    @test order_named_inds(Val((:x, :y)); x=30, y=20) == (30, 20)
end
@testset "allocations: order_named_inds" begin
    @test 0 == @ballocated (()->order_named_inds(Val((:a, :b, :c)); b=1, c=2))()
    @test 0 == @ballocated (()->order_named_inds(Val((:a, :b, :c)), (b=1, c=2)))()
end

@testset "remaining_dimnames_from_indexing" begin
    @test remaining_dimnames_from_indexing((:a, :b, :c), (10,20,30)) === ()
    @test remaining_dimnames_from_indexing((:a, :b, :c), (10, :,30)) == (:b,)
    @test remaining_dimnames_from_indexing((:a, :b, :c), (1:1, [true], [20])) == (:a, :b, :c)

    @test remaining_dimnames_from_indexing((:a, :b, :c), (:, [CartesianIndex()], :, :)) == (:a, :_, :b, :c)
    @test remaining_dimnames_from_indexing((:a, :b, :c), (1, [CartesianIndex()], 2, :)) == (:_, :c)
    @test remaining_dimnames_from_indexing((:a, :b, :c), (CartesianIndex(1,1), :)) == (:c,)

    # Cases which drop dimension names
    logical = Base.to_indices(ones(2,3), (rand(2,3).>0.5,))
    @test remaining_dimnames_from_indexing((:a, :b), logical) === ()
    cart_vec = findall(x -> x>0.5, rand(2,3))
    @test remaining_dimnames_from_indexing((:a, :b), (cart_vec,)) === ()
end
@testset "allocations: remaining_dimnames_from_indexing" begin
    @test 0 == @ballocated (()->remaining_dimnames_from_indexing((:a, :b, :c), (:,390,:)))()
end

@testset "remaining_dimnames_after_dropping" begin
    @test remaining_dimnames_after_dropping((:a, :b, :c), 1) == (:b, :c)
    @test remaining_dimnames_after_dropping((:a, :b, :c), 3) == (:a, :b)
    @test remaining_dimnames_after_dropping((:a, :b, :c), (1,3)) == (:b,)
    @test remaining_dimnames_after_dropping((:a, :b, :c), (3,1)) == (:b,)
end
@testset "allocations: remaining_dimnames_after_dropping" begin
    @test 0 == @ballocated remaining_dimnames_after_dropping((:a,:b,:c,:d,:e), 4)
    @test 0 == @ballocated (()->remaining_dimnames_after_dropping((:a,:b,:c,:d,:e), (1,3)))()
end

@testset "permute_dimnames" begin
    @test permute_dimnames((:a, :b, :c), (1, 2, 3)) == (:a ,:b, :c)
    @test permute_dimnames((:a, :b, :c), (3, 2, 1)) == (:c ,:b, :a)

    # permute_dimnames allows non-bijective "permutations"
    @test permute_dimnames((:a, :b, :c), (3, 3, 1)) == (:c ,:c, :a)
    @test permute_dimnames((:a, :b, :c), (3, 1))== (:c, :a)

    @test_throws BoundsError permute_dimnames((:a, :b, :c), (30, 30, 30))
    @test_throws BoundsError permute_dimnames((:a, :b), (1, 0))
end
@testset "allocations: permute_dimnames" begin
    if VERSION >= v"1.1"
        @test 0 == @ballocated permute_dimnames((:a,:b,:c), (1,3,2))
    end
end

@testset "tuple_issubset" begin
    @test tuple_issubset((:a, :c), (:a, :b, :c)) == true
    @test tuple_issubset((:a, :b, :c), (:a, :c)) == false
end
@testset "allocations: tuple_issubset" begin
    @test 0 == @ballocated tuple_issubset((:a, :c), (:a, :b, :c))
    @test 0 == @ballocated tuple_issubset((:a, :b, :c), (:a, :c))
end

@testset "tuple_cat" begin
    @test tuple_cat((1, 2), (3, 4, 5), (6,)) == (1, 2, 3, 4, 5, 6)
    @test tuple_cat((1, 2)) == (1, 2)
end
@testset "allocations: tuple_cat" begin
    @test 0 == @ballocated tuple_cat((1, 2), (3, 4, 5), (6,))
end


text:
"""
    countmap(iter; datatype=Any, dicttype=Dict)

Return an `AbstractDict` counting the number of occurences of elements
in `iter`.

Unlike `StatsBase.countmap`, `DataUtils.countmap`
does not require that one first calls `collect` on `iter` if
`iter` is not an `AbstractArray`. `DataUtils.countmap` accepts
any iterator, including one that may generate more values than can
be stored in memory.

`DataUtils.countmap` does not implement a fast, radix sort
as `StatsBase.countmap` does.

# Example
```julia
julia> import StatsBase; import DataUtils;

julia> d = Dict(rand(1:1000) => rand(1:5) for i in 1:100);

julia> StatsBase.countmap(values(d))
ERROR: MethodError: no method matching countmap(::Base.ValueIterator{Dict{Int64,Int64}})
...

julia> DataUtils.countmap(values(d))
Dict{Any,Int64} with 5 entries:
  4 => 19
  2 => 21
  3 => 20
  5 => 11
  1 => 25
```
"""
function countmap(iter; datatype=Any, dicttype=Dict)
    return countmap!(dicttype{datatype,Int}(), iter)
end

"""
    countmap(iter::AbstractArray{T}; dicttype=Dict) where {T}

Equivalent to calling `countmap(iter; datatype=Any, dicttype=Dict)`,
with `datatype=T`.
"""
function countmap(a::AbstractArray{T}; dicttype=Dict) where {T}
    return countmap!(dicttype{T,Int}(), a)
end

function countmap!(cm::Dict, iter)
    for v in iter
        index = Base.ht_keyindex2!(cm, v)
        if index > 0
            @inbounds cm.vals[index] += 1
        else
            @inbounds Base._setindex!(cm, 1, v, -index)
        end
    end
    return cm
end

function countmap!(cm::AbstractDict, iter)
    for v in iter
        if haskey(cm, v)
            @inbounds cm[v] += 1
        else
            @inbounds cm[v] = 1
        end
    end
    return cm
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    GMSH

Ordering convention of the GMSH project:
https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering
"""
struct GMSH <: OrderingConvention end

# Triangle
# v
# ^
# |
# 2
# |`\
# |  `\
# |    `\
# |      `\
# |        `\
# 0----------1 --> u

connectivities(::Type{<:Triangle}, ::Val{1}, ::Type{GMSH}) =
  connect.([(1, 2), (2, 3), (3, 1)], Segment)

# Quadrangle
#       v
#       ^
#       |
# 3-----------2
# |     |     |
# |     |     |
# |     +---- | --> u
# |           |
# |           |
# 0-----------1

connectivities(::Type{<:Quadrangle}, ::Val{1}, ::Type{GMSH}) =
  connect.([(1, 2), (2, 3), (3, 4), (4, 1)], Segment)

# Pyramid
#                4
#              ,/|\
#            ,/ .'|\
#          ,/   | | \
#        ,/    .' | `.
#      ,/      |  '.  \
#    ,/       .' w |   \
#  ,/         |  ^ |    \
# 0----------.'--|-3    `.
#  `\        |   |  `\    \
#    `\     .'   +----`\ - \ -> v
#      `\   |    `\     `\  \
#        `\.'      `\     `\`
#           1----------------2
#                     `\
#                        u

function connectivities(::Type{<:Pyramid}, ::Val{2}, ::Type{GMSH})
  base = connect((1,2,3,4), Quadrangle)
  side = connect.([(1,2,5),(2,3,5),(3,4,5),(4,1,5)], Triangle)
  [base; side]
end

# Tetrahedron
#                    v
#                  .
#                ,/
#               /
#            2
#          ,/|`\
#        ,/  |  `\
#      ,/    '.   `\
#    ,/       |     `\
#  ,/         |       `\
# 0-----------'.--------1 --> u
#  `\.         |      ,/
#     `\.      |    ,/
#        `\.   '. ,/
#           `\. |/
#              `3
#                 `\.
#                    ` w

connectivities(::Type{<:Tetrahedron}, ::Val{2}, ::Type{GMSH}) =
  connect.([(1,2,3),(1,4,3),(1,4,2),(2,3,4)], Triangle)

# Hexahedron
#        v
# 3----------2
# |\     ^   |\
# | \    |   | \
# |  \   |   |  \
# |   7------+---6
# |   |  +-- |-- | -> u
# 0---+---\--1   |
#  \  |    \  \  |
#   \ |     \  \ |
#    \|      w  \|
#     4----------5

connectivities(::Type{<:Hexahedron}, ::Val{2}, ::Type{GMSH}) =
  connect.([(1,2,3,4),(5,6,7,8),(1,5,6,2),
            (2,6,7,3),(3,4,8,7),(1,5,8,4)], Quadrangle)


text:
using StructJuMP

m = StructuredModel()

n = 10
@defVar(m, x[1:n] >= 0)
@defVar(m, y[1:n] >= 0)
for i in 1:n
    @addConstraint(m, x[i] - y[i] <= i)
end
setObjective(m, :Min, x[1] - y[3] + 2x[6] + x[1]*x[2] - 2y[4]*x[3])

numScen = 8

for i in 1:numScen
    bl = StructuredModel(parent=m)
    @defVar(bl, w[1:n] >= 0)
    @addConstraint(bl, sum{i*w[i], i=1:n; iseven(i)} == 1)
    @addConstraint(bl, sum{x[i] - 2w[n-i+1], i=1:n} <= 0)
    setObjective(bl, :Min, w[5] - w[1]*w[1] + 2w[2]*w[3])
end


text:
@node Bernoulli Stochastic [ out, (p, aliases = [ θ ]) ]

conjugate_type(::Type{ <: Bernoulli }, ::Type{ Val{ :out } }) = Bernoulli
conjugate_type(::Type{ <: Bernoulli }, ::Type{ Val{ :θ } })   = Beta

@average_energy Bernoulli (q_out::Any, q_p::Any) = -mean(q_out) * logmean(q_p) - (1.0 - mean(q_out)) * mirroredlogmean(q_p)


text:
using GasSteadySim
using Distributions
using DelimitedFiles
using Random

function run_gf_solver(data::Dict{String,Any}, eos::Symbol; kwargs...)
    ss = initialize_simulator(data, eos=eos)
    solver_return = run_simulator!(ss; kwargs...)
    # if Int(solver_return.status) == 2
    #     @show ss.sol["compressor_flow"][1] # this compressor is not in cycle
    # end
    return solver_return
end

function perturb_compressor_ratios!(data::Dict{String,Any})
    for (_, compressor) in get(data, "boundary_compressor", [])
        compressor["value"] = rand(Uniform(1.1, 1.4))
    end
end

function perturb_injections!(data::Dict{String,Any})
    for (i, nonslack_flow) in get(data, "boundary_nonslack_flow", [])
        scaling = rand(Uniform(0.75, 1.25))
        data["boundary_nonslack_flow"][i] = scaling * nonslack_flow
    end
end

function make_nominal_values_unity!(data::Dict{String,Any})
    data["simulation_params"]["nominal_length"] = 1.0
    data["simulation_params"]["nominal_velocity"] = 1.0
    data["simulation_params"]["nominal_pressure"] = 1.0
    data["simulation_params"]["nominal_density"] = 1.0
end

function write_results(results::Dict{Int,Any}, num_runs, output_path)
    to_write_header = ["case_name", "run_id", "eos", "solver_status",
        "num_iterations", "run_time", "feasibility"]
    to_write = Array{Any,2}(undef, num_runs, 7)
    for i in 1:num_runs
        result = results[i]
        case_name = result["case_name"]
        eos = result["eos"]
        status = result["solver_status"]
        iter = result["num_iterations"]
        time = result["time"]
        feas = result["feasibility"]
        to_write[i, :] = [case_name, i, eos, status, iter, time, feas]
    end
    open(output_path, "w") do io
        writedlm(io, [permutedims(to_write_header); to_write], ',')
    end
end

function paper_runs(folder::AbstractString, eos::Symbol; num_runs = 500,
    output_folder = "./output/", output_file = "8-node-ideal.csv", case_name = "8-node",
    dimensional = false)
    data_sample = GasSteadySim._parse_data(folder)
    # julia compiling run
    solver_return = run_gf_solver(data_sample, eos, method = :newton, show_trace = true)
    Random.seed!(2022)
    results = Dict{Int,Any}()

    for i in 1:num_runs
        results[i] = Dict{String,Any}()
        data = GasSteadySim._parse_data(folder)
        (dimensional == true) && (make_nominal_values_unity!(data))
        perturb_compressor_ratios!(data)
        perturb_injections!(data)
        solver_return = run_gf_solver(data, eos, method = :newton, show_trace = false)
        results[i]["solver_status"] = Int(solver_return.status)
        results[i]["num_iterations"] = solver_return.iterations
        results[i]["time"] = solver_return.time
        results[i]["feasibility"] =
            if Int(solver_return.status) in [2, 3]
                "false"
            elseif Int(solver_return.status) == 1
                "unknown"
            else
                "true"
            end
        results[i]["case_name"] = case_name
        results[i]["eos"] = string(eos)
    end
    write_results(results, num_runs, output_folder * output_file)
    return
end

function get_max_deviation(ideal, cnga)
    p_ideal = ideal.sol["nodal_pressure"]
    p_cnga = cnga.sol["nodal_pressure"]
    p_deviation =  maximum([abs(p_ideal[i] - p_cnga[i])/p_cnga[i]*100.0 for i in keys(p_ideal)])
    rho_ideal = ideal.sol["nodal_density"]
    rho_cnga = cnga.sol["nodal_density"]
    rho_deviation =  maximum([abs(rho_ideal[i] - rho_cnga[i])/rho_cnga[i]*100.0 for i in keys(p_ideal)])
    return p_deviation, rho_deviation
end

function ideal_vs_non_ideal_runs()
    cases = ["GasLib-11", "GasLib-24", "GasLib-40", "GasLib-134"]
    output_file_p = "./output/ideal_vs_nonideal_max_deviation_pressure.csv"
    output_file_rho = "./output/ideal_vs_nonideal_max_deviation_density.csv"
    data_folder = "./data/"
    num_runs = 500
    data_sample = GasSteadySim._parse_data("./data/2-node/")
    # julia compiling run
    _ = run_gf_solver(data_sample, :ideal, method = :newton, show_trace = true)

    Random.seed!(2000)
    results_p = Dict{String,Any}(case => [] for case in cases)
    results_rho = Dict{String,Any}(case => [] for case in cases)

    for case in cases
        for _ in 1:num_runs
            data_i = GasSteadySim._parse_data(data_folder * case * "/")
            perturb_compressor_ratios!(data_i)
            perturb_injections!(data_i)
            data_n = deepcopy(data_i)
            ss_i = initialize_simulator(data_i, eos=:ideal)
            ss_n = initialize_simulator(data_n, eos=:full_cnga)
            _ = run_simulator!(ss_i)
            _ = run_simulator!(ss_n)
            deviation_p, deviation_rho = get_max_deviation(ss_i, ss_n)
            push!(results_p[case], deviation_p)
            push!(results_rho[case], deviation_rho)
        end
    end

    header = cases
    data_p = Array{Any,2}(undef, num_runs, 4)
    data_rho = Array{Any,2}(undef, num_runs, 4)
    for i in 1:num_runs
        data_p[i, :] = [results_p[case][i] for case in cases]
        data_rho[i, :] = [results_rho[case][i] for case in cases]
    end

    open(output_file_p, "w") do io
        writedlm(io, [permutedims(header); data_p], ',')
    end

    open(output_file_rho, "w") do io
        writedlm(io, [permutedims(header); data_rho], ',')
    end

end

max_stars_repo_name: kaarthiksundar/GasFlowRuns

max_issues_repo_name: kaarthiksundar/GasFlowRuns

max_forks_repo_name: kaarthiksundar/GasFlowRuns


text:
using FilingForest
using FilingForest.Tokens
using FilingForest.Parser
import FilingForest.Parser: tokenstring, word, delimiter, inline

@test tokenize(tokenstring, "") == Token[]
@test tokenize(tokenstring, "a") == Token[l"a"]

@test convert(NamedTuple{(:a,),Tuple{Int}},(a=2,b=1)) == (a=2,)

@test tokenize("") == Token[]
@test tokenize("test eines") ==  Token[l"test", delim" ",l"eines"]
@test tokenize(tokenstring,"test eines") ==  Token[l"test", delim" ",l"eines"]

@test tokenize("a") == [ Token(:literal, "a") ]
@test tokenize("a b") == [ l"a",  delim" ", l"b" ]
@test tokenize(".a b") == [ delim".", l"a",  delim" ", l"b" ]
@test tokenize(".aü öb") == [ delim".", l"aü",  delim" ", l"öb" ]

@test tokenize("/a") == [ delim"/", l"a" ]
@test tokenize("/a/b/") == [ delim"/", l"a", delim"/", l"b", delim"/" ]

## TODO

## ?? @test tokenize("<a=b>") == [ Token(:a, "b") ]
## ?? @test tokenize("<b>") == [ Token(:b) ]
## ?? @test tokenize("<b=test>") == [ Token(:b,"test") ]

## ?? @test tokenize("<b=test>.pdf", V) == [ Token(:b,"test"), (@T ext ".pdf") ]

## ?? @test tokenize("/a/<b>") == [ delim"/", l"a", delim"/", Token(:b) ]

if false
    (tokenize("years like 2012 are parsed"))

    @test tokenize(l"a.b") == [ l"a", Token(:ext,".b") ]
    @test tokenize("years are parsed like 2012")[end] == Token(:yyyy,"2012")

    @test (@T yyyy 2012) == Token(:yyyy,"2012")

    @test tokenize("invoice 2012.pdf") == [ l"invoice", delim" ",  (@T yyyy 2012), (@T ext ".pdf") ]
    @test tokenize("invoice 2012.pdf.pdf") ==
        [ l"invoice", delim" ",  (@T yyyy 2012), delim".", l"pdf", (@T ext ".pdf") ]

    @test tokenize("1.2.3.7z",V) == [ (@T version "1.2.3"), @T ext ".7z" ]
    @test tokenize("draft 99.pdf",V) ==
        [ l"draft", delim" ", (@T version "99"), @T ext ".pdf" ]
end
# adding pattern to Filing applies new pattern to Files
# adding File to Filing applies old pattern to File

# query Files matching a pattern

# query Files matching a pattern

import FilingForest.Parser: bracket_number, bracket_reference

@test (tokenize(bracket_number,"[1]")) ==
    Token(:number, "[1]")

@test (tokenize(bracket_reference,"[1, 2]")) ==
    Token(:reference, "[1, 2]")


text:
#Reference
#岑松, 龙志飞, 对转角场和剪应变场进行合理插值的厚板元[J]. 工程力学, 1998,15(3): 1―14.
function K_TMT(elm::Tria)::Matrix{Float64}
    E₀,ν₀=elm.material.E,elm.material.ν
    G₀=elm.material.G
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc-center)
    x₂,y₂,z₂=T*(elm.node2.loc-center)
    x₃,y₃,z₃=T*(elm.node3.loc-center)
    K=Matrix{Float64}(undef,12,12)
    A=elm.A
    a₁,a₂,a₃=[xᵢ[j]*yᵢ[m]-xᵢ[m]*yᵢ[j] for (j,m) in zip([2,3,1],[3,1,2])]
    b₁,b₂,b₃=[yᵢ[j]-yᵢ[m] for (j,m) in zip([2,3,1],[3,1,2])]
    c₁,c₂,c₃=[-xᵢ[j]+xᵢ[m] for (j,m) in zip([2,3,1],[3,1,2])]
    l₁=norm(elm.node3.loc-elm.node2.loc)
    l₂=norm(elm.node1.loc-elm.node3.loc)
    l₃=norm(elm.node1.loc-elm.node2.loc)

    δ₁=(t/d₁)^2/(5/6*(1-ν₀)+2*(t/d₁)^2)
    δ₂=(t/d₂)^2/(5/6*(1-ν₀)+2*(t/d₂)^2)
    δ₃=(t/d₃)^2/(5/6*(1-ν₀)+2*(t/d₃)^2)

    F=-3/2/A*[c₁/l₁^2*(b₂*L₃+b₃*L₂) c₂/l₂^2*(b₃*L₁+b₁*L₃) c₃/l₃^2*(b₁*L₂+b₂*L₁);
             -b₁/l₁^2*(c₂*L₃+c₃*L₂) -b₂/l₂^2*(c₃*L₁+c₁*L₃) -b₃/l₃^2*(c₁*L₂+c₂*L₁);
             M₁ M₂ M₃]

    Δ=[1-2δ₁ 0 0;
       0 1-2δ₂ 0;
       0 0 1-2δ₃]

    G=[0 0 0 -2 -c₁ b₁ 2 -c₁ b₁;
       2 -c₂ b₂ 0 0 0 -2 -c₂ b₂;
       -2 -c₃ b₃ 2 -c₃ b₃ 0 0 0]

    Bb⁰=-1/2/A*[0 b₁ 0 0 b₂ 0 0 b₃ 0;
               0 0 c₁ 0 0 c₂ 0 0 c₃;
               0 c₁ b₁ 0 c₂ b₂ 0 c₃ b₃]

    FDG=F*Δ*G

    Bb=Bb⁰+FDG

    Δʼ=[δ₁ 0 0;
        0 δ₂ 0;
        0 0 δ₃]

    D₀=E₀*t^3/12/(1-ν₀^2)
    D=D₀*[1 ν₀ 0;
          ν₀ 1 0;
          0  0 (1-ν₀)/2]
    k=5/6

    function BtDB(x)
        x,y=x[1],x[2]
        L₁=1/2/A*(a₁+b₁*x+c₁*y)
        L₂=1/2/A*(a₂+b₂*x+c₂*y)
        L₃=1/2/A*(a₃+b₃*x+c₃*y)

        M₁=1/l₁^2*((c₁*c₂-b₁*b₂)*L₃+(c₃*c₁-b₃*b₁)*L₂)
        M₂=1/l₂^2*((c₂*c₃-b₂*b₃)*L₁+(c₁*c₂-b₁*b₂)*L₃)
        M₁=1/l₃^2*((c₃*c₁-b₃*b₁)*L₂+(c₂*c₃-b₂*b₃)*L₁)

        H=1/2/A*[b₃*L₂-b₂*L₃ b₁*L₃-b₃*L₁ b₂*L₁-b₁*L₂;
                 c₃*L₂-c₂*L₃ c₁*L₃-c₃*L₁ c₂*L₁-c₁*L₂]

        Bs=H*Δʼ*G

        C=k*G₀*t*[1 0;0 1]

        Kb=transpose(Bb)*D*Bb
        Ks=transpose(Bs)*C*Bs
        K=Kb+Ks
    end
    K=hcubature(BtDB,[-1,-1],[1,1])[1]
    # #left-hand system to right-hand system
    # I=1:9
    # J=[1,3,2,4,6,5,7,9,8,10,12,11]
    # L=sparse(I,J,1.,12,12)
    # for i in [1,3,4,6,7,9,10,12]
    #     for j in [2,5,8,11]
    #         K[i,j]=K[j,i]=-K[i,j]
    #     end
    # end
    # K=L'*sparse(K)*L
    #9x9 to 18x18
    I=1:9
    J=[3,4,5,9,10,11,15,16,17]
    L=sparse(I,J,1.,9,18)
    return L'*K*L
end