text:
using Flux, Statistics, Trebuchet
using Flux.Tracker: data, gradient
using Flux.Optimise: update!

#using CuArrays

using BSON: @save, @load

using Distributions: Normal, sample
using DataStructures: CircularBuffer
using Printf

#=
# Description of the problem:
    There is Trebuchet, which throws a mass to a target. The mass is to be
    released at an angle, and at certain velocity so that it lands on the target.
    The velocity of release is determined by the counterweight of the Trebuchet.
    Given conditions of environment we are required to predict the angle of
    release and counterweight.

    The problem is gamified by introducing a threshold. The player gets 99
    attempts. The first attempt has threshold value of 1. Threshold determines
    the tolerable relative error between actual distance travelled by the Projectile
    and the target distance. This relative error has to be less than the threshold
    in order to continue playing the game. The game gets tougher at each attempts
    by reduction of threshold by 0.01.

# Input:  Wind speed,   Target distance
# Output: ReleaseAngle, Weight
=#

# --------------------------------- Parameters ------------------------------- #

STATE_SIZE = 3
ACTION_SIZE = 2
BATCH_SIZE = 128
MEM_SIZE = 1000000
MAX_EP = 100000
MAX_DIST  = 500	# Maximum target distance
MAX_SPEED =  10 # Maximum wind speed
START_TRAINING_AFTER = 10000

# ----------------------------- Hyperparameters ------------------------------ #

γ = 99f-2     	# Discount rate
τ = 1f-2 		# For running average while updating target networks
η_act = 1f-4    # Learning rate of actor
η_crit = 1f-3	# Learning rate of critic
noise = Normal(0, 1)
noise_scale = 1f-1

# ----------------------------- Model Architecture --------------------------- #

actor = Chain(Dense(STATE_SIZE, 256, relu), Dense(256, 512, relu),
	      Dense(512, 256, relu), Dense(256, ACTION_SIZE, relu)) |> gpu
actor_target =  deepcopy(actor)

critic = Chain(Dense(STATE_SIZE + ACTION_SIZE, 256, relu), Dense(256, 512, relu),
	       Dense(512, 256, relu), Dense(256, 1)) |> gpu
critic_target = deepcopy(critic)

# ---------------------------- Param Update Functions ------------------------ #

function nullify_grad!(p)
  if typeof(p) <: TrackedArray
    p.grad .= 0.0f0
  end
  return p
end

function zero_grad!(model)
  model = mapleaves(nullify_grad!, model)
end

function update_target!(target, model; τ = 1.0f0)
  for (p_t, p_m) in zip(params(target), params(model))
    p_t.data .= (1.0f0 - τ) * p_t.data .+ τ * p_m.data
  end
end

# -------------------------------- DDPG -------------------------------------- #

opt_crit = ADAM(η_crit)
opt_act = ADAM(η_act)

function trainer()
  # Getting data in shape
  minibatch = sample(memory, BATCH_SIZE)
  x = hcat(minibatch...)

  s = hcat(x[1, :]...) |> gpu
  a = hcat(x[2, :]...) |> gpu
  s′ = hcat(x[3, :]...) |> gpu
  r = hcat(x[4, :]...) |> gpu
  s_mask = .!hcat(x[5, :]...) |> gpu

  a′ = actor_target(s′).data
  crit_tgt_in = vcat(s′, a′)
  v′ = critic_target(crit_tgt_in).data
  y = r .+ γ * v′ .* s_mask	# set v′ to 0 where s_ is terminal state

  crit_in = vcat(s, a)
  v = critic(crit_in)
  loss_crit = Flux.mse(y, v)

  # Update Actor
  actions = actor(s)
  crit_in = param(vcat(s, actions.data))
  crit_out = critic(crit_in)

  Flux.back!(sum(crit_out))

  act_grads = -crit_in.grad[end-ACTION_SIZE+1:end, :]
  zero_grad!(actor)
  Flux.back!(actions, act_grads)  # Chain rule
  update!(opt_act, params(actor))

  # Update Critic
  zero_grad!(critic)
  Flux.back!(loss_crit)
  update!(opt_crit, params(critic))

end

# --------------------------------- Env Functions ---------------------------- #

memory = CircularBuffer{Any}(MEM_SIZE)

# stores the tuple of state, action, reward, next_state, and done
remember(state, action, next_state, reward, done) =
	push!(memory, [state, action, next_state, reward, done])

function action(state, train=true)
  state = reshape(state, size(state)..., 1)
  act_pred = cpu(data(actor(state |> gpu))) .+
  	     train * noise_scale * Float32.(rand(noise, ACTION_SIZE))
end

function reward(target_dist, threshold, wind_speed, release_angle, weight)
  ws, ra, w = Float64.([wind_speed, release_angle, weight])
  t = TrebuchetState(wind_speed=ws, release_angle=ra, weight=w)
  simulate(t)
  actual_dist = Float32(t.sol.Projectile[end][1])
  r = (target_dist - actual_dist) / target_dist
  return max(-1f0, threshold - abs(r))
end

function episode(train=true)
  total_reward = 0
  threshold = 1f0

  wind_speed  = (2rand(Float32) - 1) * MAX_SPEED
  target_dist = (2rand(Float32) - 1) * MAX_DIST
  s = [wind_speed, target_dist, 100threshold]

  attempts = 0

  while threshold > 0
    a = action(s, train)

    r = reward(target_dist, threshold, wind_speed, a[1], a[2])
    total_reward += r * (r > 0)

    threshold -= 1f-2
    s′ = r < 0 ? [wind_speed, target_dist, 100threshold] : zeros(Float32, 3)
    train && remember(s, a, s′, r, r < 0)

    s .= s′
    attempts += 1

    if train && length(memory) ≥ START_TRAINING_AFTER
      trainer()
      update_target!(actor_target, actor; τ=τ)
      update_target!(critic_target, critic; τ=τ)
    end

    r < 0 && break
  end

  total_reward, attempts
end

# ------------------------------ Saving & Loading ---------------------------- #
#=
function saveModel()
  act_wts = cpu.(Tracker.data.(params(actor)))
  @save "models/actor.bson" act_wts

  crit_wts = cpu.(Tracker.data.(params(critic)))
  @save "models/critic.bson" crit_wts
  println("Model saved")
end

function loadModel()
  @load "models/actor.bson" act_wts
  @load "models/critic.bson" crit_wts

  Flux.loadparams!(actor, gpu.(act_wts))
  Flux.loadparams!(critic, gpu.(crit_wts))
  println("Model loaded")
end
=#
# ------------------------------- Training ----------------------------------- #

scores   = CircularBuffer{Float32}(1000)
attempts = CircularBuffer{Int}(1000)

println("Training starts...")
for e = 1:MAX_EP
  total_reward, att = episode()
  push!(scores, total_reward)
  push!(attempts, att)
  avg_score, avg_attempts = mean.([scores, attempts])

  print("Episode: $(@sprintf("%6d", e)) || ")
  print("Score: $(@sprintf("%8.5f", total_reward)) | ")
  print("Avg score: $(@sprintf("%7.5f", avg_score)) || ")
  print("Attempts: $(@sprintf("%2d", att)) | ")
  println("Avg attempts: $(@sprintf("%7.5f", avg_attempts)) || ")

  #saveModel()
end

println("Training ends...")

#------------------------------------- Testing ------------------------------- #
empty!(scores)
empty!(attempts)

println("Testing starts...")

for e = 1:1000
  total_reward, att = episode(false)
  push!(scores, total_reward)
  push!(attempts, att)
  avg_score, avg_attempts = mean.([scores, attempts])

  print("Episode: $(@sprintf("%6d", e)) || ")
  print("Score: $(@sprintf("%8.5f", total_reward)) | ")
  print("Avg score: $(@sprintf("%7.5f", avg_score)) || ")
  print("Attempts: $(@sprintf("%2d", att)) | ")
  println("Avg attempts: $(@sprintf("%7.5f", avg_attempts))")
end

max_stars_repo_name: LislPisl/model-zoo

max_issues_repo_name: deepaksuresh/model-zoo

max_forks_repo_name: deepaksuresh/model-zoo


text:
using Rethinking
using Documenter

makedocs(;
    modules=[Rethinking],
    authors="andreas koher",
    repo="https://github.com/andreasko/Rethinking.jl/blob/{commit}{path}#L{line}",
    sitename="Rethinking.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://andreasko.github.io/Rethinking.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/andreasko/Rethinking.jl",
)

ext: jl


text:
"""
    KaplanMeier

An immutable type containing survivor function estimates computed
using the Kaplan-Meier method.
The type has the following fields:

* `times`: Distinct event times
* `nevents`: Number of observed events at each time
* `ncensor`: Number of right censored events at each time
* `natrisk`: Size of the risk set at each time
* `survival`: Estimate of the survival probability at each time
* `stderr`: Standard error of the log survivor function at each time

Use `fit(KaplanMeier, ...)` to compute the estimates and construct
this type.
"""
struct KaplanMeier{T<:Real} <: NonparametricEstimator
    times::Vector{T}
    nevents::Vector{Int}
    ncensor::Vector{Int}
    natrisk::Vector{Int}
    survival::Vector{Float64}
    stderr::Vector{Float64}
end

estimator_start(::Type{KaplanMeier}) = 1.0  # Estimator starting point
stderr_start(::Type{KaplanMeier}) = 0.0 # StdErr starting point

estimator_update(::Type{KaplanMeier}, es, dᵢ, nᵢ) = es * (1 - dᵢ / nᵢ) # Estimator update rule
stderr_update(::Type{KaplanMeier}, gw, dᵢ, nᵢ) = gw + dᵢ / (nᵢ * (nᵢ - dᵢ)) # StdErr update rule

"""
    confint(km::KaplanMeier, α=0.05)

Compute the pointwise log-log transformed confidence intervals for the survivor
function as a vector of tuples.
"""
function StatsBase.confint(km::KaplanMeier, α::Float64=0.05)
    q = quantile(Normal(), 1 - α/2)
    return map(km.survival, km.stderr) do srv, se
        l = log(-log(srv))
        a = q * se / log(srv)
        exp(-exp(l - a)), exp(-exp(l + a))
    end
end

"""
    fit(KaplanMeier, times, status) -> KaplanMeier

Given a vector of times to events and a corresponding vector of indicators that
dictate whether each time is an observed event or is right censored, compute the
Kaplan-Meier estimate of the survivor function.
"""
StatsBase.fit(::Type{KaplanMeier}, times, status)


text:
"""
delete_nodes!(mtg::Node,<keyword arguments>)

Delete nodes in mtg following filters rules.

# Arguments

## Mandatory arguments

- `node::Node`: The node to start at.

## Keyword Arguments (filters)

- `scale = nothing`: The scale to delete. Usually a Tuple-alike of integers.
- `symbol = nothing`: The symbol to delete. Usually a Tuple-alike of Strings.
- `link = nothing`: The link with the previous node to delete. Usually a Tuple-alike of Char.
- `all::Bool = true`: Continue after the first deletion (`true`), or stop?
- `filter_fun = nothing`: Any filtering function taking a node as input, e.g. [`isleaf`](@ref)
to decide whether to delete a node or not.
- `child_link_fun = new_child_link`: a function that takes the child node of a deleted node
as input and returns its new link (see details).

# Notes

1. The function is acropetal, meaning it will apply the deletion from leaves to the root to ensure
that one pass is enough and we don't repeat the process of visiting already visited children.
2. The function does not do anything fancy, it let the user take care of its own rules when
deleting nodes. So if you delete a branching node, the whole subtree will be modified and take
the link of the children. This process is left to the user becaue it highly depends on the mtg
structure.
3. The package provides some pre-made functions for filtering. See for example [`is_segment!`](@ref)
to re-compute the mtg at a given scale to have only nodes at branching points. This is often used
to match automatic reconstructions from e.g. LiDAR point cloud with manual measurements.
4. The default function used for `child_link_fun` is [`new_child_link`](@ref), which tries to be
clever considering the parent and child links. See its help page for more information. If the
link shouldn't be modified, use the following function instead:
```julia
node -> node.MTG.link
```

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","A1B1.mtg")
mtg = read_mtg(file)

delete_nodes!(mtg, scale = 2) # Will remove all nodes of scale 2

# Delete the leaves:
delete_nodes!(mtg, symbol = "Leaf")
# Delete the leaves and internodes:
delete_nodes!(mtg, symbol = ("Leaf","Internode"))

# Make the mtg match field measurements made only at branching points for the scales 1 + 2:
mtg = delete_nodes!(mtg, filter_fun = is_segment!, scale = 2)
```
"""
function delete_nodes!(
    node;
    scale = nothing,
    symbol = nothing,
    link = nothing,
    all::Bool = true, # like continue in the R package, but actually the opposite
    filter_fun = nothing,
    child_link_fun = new_child_link
)

    # Check the filters once, and then compute the descendants recursively using `descendants_`
    check_filters(node, scale = scale, symbol = symbol, link = link)
    filtered = is_filtered(node, scale, symbol, link, filter_fun)

    if filtered
        node = delete_node!(node)
        # Don't go further if all == false
        all ? nothing : return nothing
    end

    delete_nodes!_(node, scale, symbol, link, all, filter_fun, child_link_fun)

    return node
end

function delete_nodes!_(node, scale, symbol, link, all, filter_fun, child_link_fun)
    if !isleaf(node)
        # First we apply the algorithm recursively on the children:
        for chnode in ordered_children(node)
            delete_nodes!_(chnode, scale, symbol, link, all, filter_fun, child_link_fun)
        end
    end

    # Then we work on the node itself. This ensures that its children will not be deleted
    # afterwards (the deletion is acropetal, i.e. from leaves to root)

    # Is there any filter happening for the current node? (true is deleted):
    filtered = is_filtered(node, scale, symbol, link, filter_fun)

    if filtered
        node = delete_node!(node, child_link_fun = child_link_fun)
        # Don't go further if all == false
        all ? nothing : return nothing
    end
end

"""
delete_node!(node; child_link_fun = new_child_link)

Delete a node and re-parent the children to its own parent.

If the node is a root and it has only one child, the child becomes the root, if it has
several children, it returns an error.

`child_link_fun` is a function that takes the child node of a deleted node as input and
returns its new link. The default function is [`new_child_link`](@ref), which tries to be
clever considering the parent and child links. See its help page for more information. If the
link shouldn't be modified, use the following function instead:

```julia
node -> node.MTG.link
```

The function returns the parent node (or the new root if the node is a root)
"""
function delete_node!(node; child_link_fun = new_child_link)
    if isroot(node)
        if length(node.children) == 1
            # If it has only one child, make it the new root:
            chnode = children(node)[1]
            chnode.parent = nothing
            # Add to the new root the mandatory root attributes:
            root_attrs = Dict(
                :symbols => node[:symbols],
                :scales => node[:scales],
                :description => node[:description]
            )

            append!(chnode, root_attrs)

            chnode.MTG.link = child_link_fun(chnode)

            node_return = chnode
        else
            error("Can't delete the root node if it has several children")
        end
    else
        parent_node = parent(node)

        if !isleaf(node)
            # We re-parent the children to the parent of the node.
            for chnode in ordered_children(node)
                # Updating the link of the children:
                chnode.MTG.link = child_link_fun(chnode)
                addchild!(parent_node, chnode; force = true)
            end
        end

        # Delete the node as child of his parent:
        pop!(parent_node.children, node.id)

        node_return = parent_node
    end

    node.parent = nothing
    node.children = nothing
    node = nothing

    return node_return
end

"""
    new_child_link(node)

Compute the new link of the child node when deleting a parent node. The rule is to give the child
node link of its parent node that is deleted, except when the parent was following its own
parent.

The node given as input is the child node here.

The rule is summarized in the following table:

|Deleted node link|Child node link|New child node link|warning|
|:---------------:|:-------------:|:-----------------:|:-----:|
|/                |/              |/                  |       |
|/                |+              |+                  |yes (1)|
|/                |<              |/                  |       |
|+                |/              |/                  |yes (2)|
|+                |+              |+                  |       |
|+                |<              |+                  |       |
|<                |/              |/                  |       |
|<                |+              |+                  |       |
|<                |<              |<                  |       |

The warnings happens when there is no satisfactory way to handle the new link, *i.e.* when
mixing branching and change in scale.

Note that in the case (1) of the warning the first child only takes the "/" link, the others
keep their links.
"""
function new_child_link(node)

    deleted_link = parent(node).MTG.link
    child_link = node.MTG.link

    if deleted_link == "+" && child_link == "/"
        new_child_link = child_link
        @warn join(
            [
            "Scale of the child node decomposed but its deleted parent was branching.",
            " Keep decomposition, please check if the branching is still correct."
        ]
        ) deleted_link child_link new_child_link
    elseif deleted_link == "+" && child_link == "<"
        new_child_link = deleted_link
    elseif deleted_link == "/" && child_link == "+"
        new_child_link = child_link
        @warn join(
            [
            "Scale of the child node branched but its deleted parent was decomposing.",
            " Keep branching, please check if the decomposition is still correct."
        ]
        ) deleted_link child_link new_child_link
    elseif deleted_link == "/" && child_link == "<"
        new_child_link = deleted_link
    else
        new_child_link = child_link
    end

    return new_child_link
end


text:
using LSystem
using Test
using Lazy

@testset "LSystem.jl" begin

# Algae

algae_model = @lsys begin
    axiom : A
    rule  : A → AB
    rule  : B → A
end

@test (@> algae_model LState next(0) result) == "A"
@test (@> algae_model LState next(1) result) == "AB"
@test (@> algae_model LState next(2) result) == "ABA"
@test (@> algae_model LState next(3) result) == "ABAAB"

end


text:
mutable struct Explicit <: TuningStrategy end

# models! returns all available models in the range at once:
MLJTuning.models!(tuning::Explicit, model, history::Nothing,
                  state, verbosity) = state
MLJTuning.models!(tuning::Explicit, model, history,
                  state, verbosity) = state[length(history) + 1:end]

function MLJTuning.default_n(tuning::Explicit, range)
    try
        length(range)
    catch MethodError
        10
    end
end


text:
export mach_dis

function optimizeit(model::Union{FFJORDModel, MvnModel}, lss_f::Function, p::Vector{Float64})
    mdl_name = typeof(model).name.name
    res = (u=p,)
    res_l = lss_f(res.u)
    @info "Initial $mdl_name loss is $res_l"
    for optm in model.optms
        m_name = typeof(optm.method).name.name
        t0 = time()
        if optm.method isa Optim.AbstractOptimizer
            res = DiffEqFlux.sciml_train(lss_f, res.u, optm.method, model.adtype;
                progress=false,
                maxiters=optm.maxiters,
                cb=cb,
                x_tol=nothing, f_tol=nothing, g_tol=nothing,
                x_abstol=optm.atol, x_reltol=-Inf64,
                f_abstol=optm.atol, f_reltol=-Inf64,
                g_abstol=optm.atol, g_reltol=-Inf64,
                outer_x_tol=nothing, outer_f_tol=nothing, outer_g_tol=nothing,
                outer_x_abstol=optm.atol, outer_x_reltol=-Inf64,
                outer_f_abstol=optm.atol, outer_f_reltol=-Inf64,
                outer_g_abstol=optm.atol, outer_g_reltol=-Inf64,
                f_calls_limit=0, g_calls_limit=0, h_calls_limit=0,
                allow_f_increases=true, allow_outer_f_increases=true,
                successive_f_tol=typemax(Int64),
                # iterations=optm.maxiters,
                outer_iterations=optm.maxiters,
                store_trace=true,
                trace_simplex=true,
                show_trace=false,
                # extended_trace=true,
                show_every=1,
                # callback=cb,
                time_limit=Inf64,
            )
        else
            res = DiffEqFlux.sciml_train(lss_f, res.u, optm.method, model.adtype;
                progress=false,
                maxiters=optm.maxiters,
                cb=cb,
            )
        end
        t1 = time()
        dt = t1 - t0
        res_l = lss_f(res.u)
        @info "Using $m_name optimization method, $mdl_name loss reached to $res_l, in $dt seconds"
    end
    @info "Final $mdl_name loss is $res_l"
    res
end

function lin_sca(xs::Vector; weights_base::Float64=2.0)
    sum(enumerate(reverse(xs))) do (i, x)
        (weights_base ^ (i - 1)) * x
    end / ((weights_base ^ size(xs, 1)) - 1)
end

function cb(p, l)
    @info "loss reached to $l"
    false
end

function mach_dis(data, pred_mach, base_mach)
    base_pdf = MLJBase.transform(base_mach, data)[!, 1]
    pred_pdf = MLJBase.transform(pred_mach, data)[!, 1]
    (
        tv_dis=totalvariation(pred_pdf, base_pdf) / size(data, 1),
        cos_dis=cosine_dist(pred_pdf, base_pdf),
        mad_dis=meanad(pred_pdf, base_pdf),
    )
end


text:
function util(c::AbstractFloat,γ::AbstractFloat)
  (c^(1.0-γ))/(1.0-γ)
end

function utilh(c::AbstractFloat,h::AbstractFloat,γ::AbstractFloat,η::AbstractFloat)
  ((c^(1.0-η)*h^η)^(1.0-γ))/(1.0-γ)
end


text:
using Documenter
using Pkg

ENV["PYTHON"] = ""
Pkg.build("PyCall")

using IMinuit

makedocs(;
    modules=[IMinuit],
    authors="Feng-Kun Guo",
    repo="https://github.com/fkguo/IMinuit.jl/blob/{commit}{path}#L{line}",
    sitename="IMinuit.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://fkguo.github.io/IMinuit.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "iminuit functions" => "iminuit.md",
        "Example" => "example.md",
        "Data and plot macros" => "data.md",
        "chisq and model_fit" => "modelfit.md",
        "Error analysis" => "contour.md"
    ],
)

deploydocs(;
    repo="github.com/fkguo/IMinuit.jl",
    target = "build",
    deps = nothing,
    make = nothing,
    branch = "gh-pages"
)


text:
using .Mux

function serve_gantt(svg_path, prof_path; port=8000, delay=5)
    # Setup Mux app
    @app gantt_app = (
        Mux.defaults,
        page(req->begin
            data = String(read(svg_path))
            html = """

            <title>
            Dagger: Scheduler Gantt Chart
            </title>
            <meta http-equiv="refresh" content="$delay">

            Go to profile data
            $data

            """
            hdrs = Dict(
                Symbol("Content-Type")=>"text/html",
                Symbol("Cache-Control")=>"no-store",
            )
            Dict(
                :body=>html,
                :headers=>hdrs,
            )
        end),
        page("/profile", req->begin
            prof = String(read(prof_path))
            html = """

            <title>
            Dagger: Scheduler Profile Flamegraph
            </title>
            <meta http-equiv="refresh" content="$delay">

            $prof

            """
            hdrs = Dict(
                Symbol("Content-Type")=>"text/html",
                Symbol("Cache-Control")=>"no-store",
            )
            Dict(
                :body=>html,
                :headers=>hdrs,
            )
        end),
        Mux.notfound()
    )

    # Start serving app
    Threads.@spawn begin
        try
            serve(gantt_app, port)
        catch err
            @error exception=(err,catch_backtrace())
        end
    end
end


text:
# lasso.jl - Lasso solvers based on FISTA and ADMM using ProximalOperators
#
#   minimize 0.5*||A*x - b||^2 + lam*||x||_1
#

using LinearAlgebra
using Random
using ProximalOperators

Random.seed!(0)

# Define solvers

function lasso_fista(A, b, lam, x; tol=1e-3, maxit=50000)
  x_prev = copy(x)
  g = NormL1(lam)
  gam = 1.0/norm(A)^2
  for it = 1:maxit
    # extrapolation step
    x_extr = x + (it-2)/(it+1)*(x - x_prev)
    # compute least-squares residual
    res = A*x_extr - b
    # compute gradient (forward) step
    y = x_extr - gam*(A'*res)
    # store current iterate
    x_prev .= x
    # compute proximal (backward) step
    prox!(x, g, y, gam)
    # stopping criterion
    if norm(x_extr-x, Inf)/gam <= tol*(1+norm(x, Inf))
      break
    end
  end
  return x
end

function lasso_admm(A, b, lam, x; tol=1e-8, maxit=50000)
  u = zero(x)
  z = copy(x)
  f = LeastSquares(A, b)
  g = NormL1(lam)
  gam = 100.0/norm(A)^2
  for it = 1:maxit
    # perform f-update step
    prox!(x, f, z - u, gam)
    # perform g-update step
    prox!(z, g, x + u, gam)
    # stopping criterion
    if norm(x-z, Inf) <= tol*(1+norm(u, Inf))
      break
    end
    # dual update
    u .+= x - z
  end
  return z
end

# Generate random problem

println("Generating random lasso problem")

m, n, k, sig = 500, 2500, 100, 1e-3
A = randn(m, n)
x_true = [randn(k)..., zeros(n-k)...]
b = A*x_true + sig*randn(m)
lam = 0.1*norm(A'*b, Inf)

# Call solvers

println("Calling solvers")

x_fista = lasso_fista(A, b, lam, zeros(n))
println("FISTA")
println("      nnz(x)    = $(norm(x_fista, 0))")
println("      obj value = $(0.5*norm(A*x_fista-b)^2 + lam*norm(x_fista, 1))")

x_admm = lasso_admm(A, b, lam, zeros(n))
println("ADMM")
println("      nnz(x)    = $(norm(x_admm, 0))")
println("      obj value = $(0.5*norm(A*x_admm-b)^2 + lam*norm(x_admm, 1))")

@test x_fista ≈ x_admm rtol=1e-3


text:
"""
altro_solver.jl
"""

"""$(TYPEDEF)
Augmented Lagrangian Trajectory Optimizer (ALTRO) is a solver developed by
the Robotic Exploration Lab at Stanford University.
The solver is special-cased to solve Markov Decision Processes
by leveraging the internal problem structure.

ALTRO consists of two "phases":
1) AL-iLQR: iLQR is used with an Augmented Lagrangian framework to solve the
problem quickly to rough constraint satisfaction
2) Projected Newton: A collocation-flavored active-set solver projects the solution
from AL-iLQR onto the feasible subspace to achieve machine-precision constraint satisfaction.
"""
struct ALTROSolver{T,S,Tp} <: ConstrainedSolver{T}
    opts::SolverOptions{T}
    stats::SolverStats{T}
    solver_al::AugmentedLagrangianSolver{T,S}
    solver_pn::Tp
    solver_name::Symbol
end

function ALTROSolver(prob::Problem{IR,T,Tm,<:Any,Tx,Tix,Tu,Tiu,Tt,TE,TM,TMd,TV},
                     opts::SolverOptions=SolverOptions(),
                     solver_uncon=iLQRSolver) where {IR,T,Tm,Tx,Tix,Tu,Tiu,Tt,TE,TM,TMd,TV}
    stats = SolverStats{T}(parent=:ALTRO)
    # make a new problem that references the old, but has ALObjective
    al_obj = ALObjective(prob, opts)
    To = typeof(al_obj)
    prob_al = Problem{IR,T,Tm,To,Tx,Tix,Tu,Tiu,Tt,TE,TM,TMd,TV}(
        prob.n, prob.m, prob.N, prob.model, al_obj, prob.convals, prob.X,
        prob.X_tmp, prob.ix, prob.U, prob.U_tmp, prob.iu, prob.ts, prob.E,
        prob.M, prob.Md, prob.V
    )
    # construct sub-solvers
    solver_al = AugmentedLagrangianSolver(prob_al, opts, stats; solver_uncon=solver_uncon)
    solver_pn = ProjectedNewtonSolver(prob_al, opts, stats)
    # put it all together
    S = typeof(solver_al.solver_uncon)
    Tp = typeof(solver_pn)
    solver = ALTROSolver{T,S,Tp}(opts, stats, solver_al, solver_pn, :ALTRO)
    return solver
end

# methods
@inline states(solver::ALTROSolver) = solver.solver_pn.X
@inline controls(solver::ALTROSolver) = solver.solver_pn.U
@inline max_violation_info(solver::ALTROSolver) = max_violation_info(solver.solver_pn.convals)

function reset!(solver::ALTROSolver)
    reset!(solver.stats, solver.opts.iterations + solver.opts.n_steps + 1, solver.solver_name)
    reset!(solver.solver_al)
end

# solve
function solve!(solver::ALTROSolver)
    reset!(solver)

    # solve with AL
    solve!(solver.solver_al)

    if status(solver) <= SOLVE_SUCCEEDED
        # Check convergence
        i = solver.solver_al.stats.iterations
        c_max = solver.solver_al.stats.c_max[i]
        if solver.opts.projected_newton && c_max > solver.opts.pn_vtol
            solve!(solver.solver_pn)
        end
        # Back-up check
        if status(solver) <= SOLVE_SUCCEEDED
            # TODO: improve this check
            if max_violation_penalty(solver.solver_pn.convals)[1] < solver.opts.pn_vtol
                solver.stats.status = SOLVE_SUCCEEDED
            end
        end
    end

    terminate!(solver)
    solver
end


text:
module TestGradientAlgebra
using Test
using LinearAlgebra
using WoodburyFactorizations
using BlockFactorizations
using CovarianceFunctions
using CovarianceFunctions: EQ, RQ, Dot, ExponentialDot, NN, GradientKernel,
        ValueGradientKernel, DerivativeKernel, ValueDerivativeKernel, input_trait,
        LazyMatrixSum, DotProductGradientKernelElement, SeparableProduct
using ForwardDiff
const AbstractMatOrFac = Union{AbstractMatrix, Factorization}

@testset "input trait property" begin
    # dot product
    k = CovarianceFunctions.Dot()^5
    h = CovarianceFunctions.Dot()^2 + 1 # quadratic
    @test input_trait(k+h) == DotProductInput()
    @test input_trait(k*h) == DotProductInput()

    # isotropic
    k = CovarianceFunctions.EQ()
    h = CovarianceFunctions.RQ(1.0) # quadratic
    @test input_trait(k+h) == IsotropicInput()
    @test input_trait(k*h) == IsotropicInput()

    # heterogeneous
    k = CovarianceFunctions.EQ()
    h = CovarianceFunctions.Dot()^2 # quadratic
    @test input_trait(k+h) == GenericInput()
    @test input_trait(k*h) == GenericInput()
end

@testset "gradient algebra" begin
    # test data
    d, n = 3, 7
    X = randn(d, n)

    # sum of homogeneous kernel types
    k = CovarianceFunctions.Dot()^5
    h = CovarianceFunctions.Dot()^2 + 1 # quadratic
    gk, gh = GradientKernel(k), GradientKernel(h)
    G = GradientKernel(k + h)
    K, Kk, Kh = gramian(G, X), gramian(gk, X), gramian(gh, X)
    MK, Mkk, Mkh = Matrix(K), Matrix(Kk), Matrix(Kh)
    @test K isa BlockFactorization
    @test K.A[1, 1] isa DotProductGradientKernelElement # this means it correctly consolidated the sum kernel into one dot product kernel
    @test MK ≈ Mkk + Mkh

    # sum of heterogeneous kernel types
    eq = CovarianceFunctions.EQ() # EQ
    line = CovarianceFunctions.Dot() # dot
    quad = line^2 + 1
    cosine = CovarianceFunctions.Cosine(randn(d))
    x, y = X[:, 1], X[:, 2]
    # TODO: Chained, VerticalScaling, Warped, SeparableSum
    heterogeneous_combinations = [eq + line, eq + quad, eq * cosine, eq * quad,
                                SeparableProduct(k, h, k), line * line]
    gradient_kernels = [GradientKernel, ValueGradientKernel]
    for gradient_kernel in gradient_kernels
        for kh in heterogeneous_combinations
            # testing product gradient kernel
            G = gradient_kernel(kh)
            kh_control = (x, y) -> (kh)(x, y)
            G_control = gradient_kernel(kh_control)
            Gxy = G(x, y)
            Gxy_control = G_control(x, y)
            @test typeof(Gxy) != typeof(Gxy_control) # means special structure was discovered
            @test Matrix(Gxy) ≈ Gxy_control
        end
    end
end

end


text:
using Documenter, RandomizedLinAlg

makedocs(
    modules = [RandomizedLinAlg],
    format = :html,
    doctest = false,
    clean = true,
    sitename = "RandomizedLinAlg.jl",
    pages = [
        "Home" => "index.md",
    ]
)

deploydocs(
    repo = "github.com/haampie/RandomizedLinAlg.jl.git",
    target = "build",
    osname = "linux",
    julia  = "1.0",
    deps = nothing,
    make = nothing,
)


text:
module mkbuild

@doc """ copy build file from destine dir to source dir
""" ->
function cpbuild(dst, src)
    build_fl = joinpath(dst, "BUILD")
    @assert isfile(build_fl)

    data = readlines(build_fl)

    srcname = basename(src)
    data = replace(data, r"name = \"[a-z]+\"", srcname)

    dstbuild = joinpath(src, "BUILD")
    #write(dstbuild, data)
end

end


text:
### Fit a full rank model with PRISMA
import FirstOrderOptimization: PRISMA, PrismaParams, PrismaStepsize

export PrismaParams, PrismaStepsize, fit_prisma!, PrismaParams

defaultPrismaParams = PrismaParams(stepsize=PrismaStepsize(Inf),
                                   maxiter=100,
                                   verbose=1,
                                   reltol=1e-5)

### FITTING
function fit_prisma!(gfrm::GFRM, params::PrismaParams = defaultPrismaParams;
			  ch::ConvergenceHistory=ConvergenceHistory("PrismaGFRM"),
			  verbose=true,
			  kwargs...)

    # we're closing over yidxs and gfrm and m and n
    yidxs = get_yidxs(gfrm.losses)
    d = maximum(yidxs[end])
    m,n = size(gfrm.A)

    # W will be the symmetric parameter; U is the upper right block
    U(W) = W[1:m, m+1:end]

    ## Grad of f
    function grad_f(W)
        G = zeros(m,d)
        Umat = U(W)
        for j=1:n
            for i in gfrm.observed_examples[j]
                # there's a 1/2 b/c 1/2 is coming from the upper right block and 1/2 from the lower left block
                G[i,yidxs[j]] = .5*grad(gfrm.losses[j], Umat[i,yidxs[j]], gfrm.A[i,j])
            end
        end
        return [zeros(m,m) G; G' zeros(d,d)]
    end

    ## Prox of g
    prox_g(W, alpha) = prox(gfrm.r, W, alpha)

    ## Prox of h
    # we're going to use a closure over gfrm.k
    # to remember what the rank of prox_h(W) was the last time we computed it
    # in order to avoid calculating too many eigentuples of W
    function prox_h(W, alpha=0; TOL=1e-10)
        while gfrm.k < size(W,1)
            l,v = eigs(Symmetric(W), nev = gfrm.k+1, which=:LR) # v0 = [v zeros(size(W,1), gfrm.k+1 - size(v,2))]
            if l[end] <= TOL
                gfrm.k = sum(l.>=TOL)
                return v*diagm(max(l,0))*v'
            else
                gfrm.k = 2*gfrm.k # double the rank and try again
            end
        end
        # else give up on computational cleverness
        l,v = eig(Symmetric(W))
        gfrm.k = sum(l.>=TOL)
        return v*diagm(max(l,0))*v'
    end

    obj(W) = objective(gfrm, W, yidxs=yidxs)

    # initialize
    # lipshitz constant for f (XXX right now a wild guess that makes sense for unscaled quadratic loss)
    Lf = 2
    if params.stepsizerule.initial_stepsize == Inf
        R = sqrt(obj(gfrm.W)*m*n/sum(map(length, gfrm.observed_examples)))/max(gfrm.r.scale,1) # estimate of distance to solution
        params.stepsizerule.initial_stepsize = 2*R/Lf
    end

    # recover
    t = time()
    gfrm.W = PRISMA(gfrm.W, Lf,
           grad_f,
           prox_g,
           prox_h,
           obj,
           params)

    # t = time() - t
    # update_ch!(ch, t, obj(W))

    gfrm.U = U(gfrm.W)

    return gfrm.U, ch
end


text:
using LightGraphs

"""
solve
"""

export solve!

function solve!(::LPModel, ::JuMP.Model, ::AbstractSolveAttribute) <: Bool end

# Will construct the constraint graph and find a feasible solution or deem the problem infeasible
# NB: Assumes the model to have difference constraints only!
function solve!(lpmodel::LPModel, model::JuMP.Model, ::ShortestPath)
    A, b = get_constraint_matrices(lpmodel)

    edges = fill((0,0), size(A)[1])   # Number of constraints

    num_vars = size(A)[2]
    distmx = zeros(num_vars, num_vars)

    # Get list of edges and distance matrix
    index = 1
    for row in eachrow(A)
        i = findall(x->x==-1, row)[1]
        j = findall(x->x==1, row)[1]
        edge = (i,j)
        edges[index] = edge

        distmx[i,j] = b[index]

        index += 1
    end

    # Directed graph
    digraph = SimpleDiGraph(Edge.(edges))
    # Use all vertices as sources, as we omit the starting source
    sources = collect(1:num_vars)

    try
        bf = bellman_ford_shortest_paths(digraph, sources, distmx)
        distances = bf.dists
        print("\nFeasible solution found:\n")
        show(distances)
        print("\n")
        # Update to feasible unbounded solution
        return true
    catch error
        if isa(error, LightGraphs.NegativeCycleError)
            # Update no feasible solution
            print("\nInfeasible solution.\n")
            return true
        else
            throw(error)
            return false
        end
    end

end

"""
function solve!(model::Model, ::ConstantObjective)
    # Set optimal objective value (if legal)
    MOI.set(model, MOI.ObjectiveValue(), objective_function(model, objective_function_type(model)).constant)

    # Check if there is a legal solution for the variables, and set result status accordingly
    for (F, S) in list_of_constraint_types(model)
        for cref in all_constraints(model, F, S)
            setVariable(model, cref, F, S)
        end
    end
end

function setVariable(model::Model, ::Any, ::VariableRef, ::Union{})

    MOI.set(model, MOI.PrimalStatus(), MOI.INFEASIBLE_POINT)
    MOI.set(model, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
end
"""

max_stars_repo_name: henriean/solver-peeker

max_issues_repo_name: henriean/solver-peeker

max_forks_repo_name: henriean/solver-peeker


text:
"Inertia matrix of a uniformly dense cuboid with mass `m`"
inertia_cuboid(length, width, thickness, m) = m / 12 * Diagonal(@SVector [width^2 + thickness^2, length^2 + thickness^2, length^2 + width^2])

"Inertia matrix of a uniformly dense thick cylinder with inner and outer radii `r1` and `r2`, height `h` and mass `m`"
inertia_thick_cylinder(r1, r2, h, m) = Diagonal(@SVector [ 1 / 12 * m * (3 * (r2^2 + r1^2) + h^2), 1 / 12 * m * (3 * (r2^2 + r1^2) + h^2), m * (r2^2 + r1^2) / 2])


text:
"""
    ncvar_def(ds, name, val, dims::Vector{<:AbstractString}, attrib = Dict();
        compress = 1, kwargs...)
    ncvar_def(ds, name, val, dims::NcDim, attrib = Dict();
        compress = 1, kwargs...)

# Arguments

- `compress`: Compression level: 0 (default) means no compression and 9 means
    maximum compression. Each chunk will be compressed individually.

- `type`: which type to save? Julia variable types (not string), e.g., `Float32`.

- `options`: Dictionary object,

    - `fillvalue`: A value filled in the NetCDF file to indicate missing data.
        It will be stored in the _FillValue attribute.

    - `chunksizes`: Vector integers setting the chunk size. The total size of a
        chunk must be less than 4 GiB.

    - `shuffle`: If true, the shuffle filter is activated which can improve the
        compression ratio.

    - `checksum`: The checksum method can be :fletcher32 or :nochecksum
        (checksumming is disabled, which is the default)

    - `attrib`: An iterable of attribute name and attribute value pairs, for
        example a Dict, DataStructures.OrderedDict or simply a vector of pairs
        (see example below)

# Examples

```julia
f = "f1.nc"
isfile(f) && rm(f)
ds = nc_open(f, "c")
compress = 1

ncdim_def(ds, "lon", lon, Dict("longname" => "Longitude", "units" => "degrees east"))
ncdim_def(ds, "lat", lat, Dict("longname" => "Latitude", "units" => "degrees north"))
ncdim_def(ds, "time", 1:ntime)

ncvar_def(ds, "HI", dat, ["lon", "lat", "time"]; compress = compress, longname = "hello")
close(ds)

# append a variable
ds = nc_open(f, "a")
ncvar_def(ds, "HI2", dat, ["lon", "lat", "time"]; compress = compress)
print(ds)
close(ds)

nc_info(f)
```

@seealso [ncdim_def](ref)
"""
function ncvar_def(ds, name, val, dims::Vector{<:AbstractString}, attrib = Dict();
    compress = 1, type = nothing, kwargs...)

    # attrib["deflatelevel"] = compress
    if name in keys(ds)
        # println(options)
        @warn "Variable `$name`: exist!"
        return
    end
    # change datatype
    if type !== nothing && eltype(val) != type
        val = @.(type(val))
    end
    # defDim(ds, name, length(val))
    defVar(ds, name, val, dims; attrib = attrib, deflatelevel = compress, kwargs...)
end

function ncvar_def(ds, name, val, dims::Vector{NcDim}, attrib = Dict();
    compress = 1, kwargs...)

    # attrib["deflatelevel"] = compress
    ncdim_def(ds, dims) # define dimensions if previous not exist
    dimnames = map(x -> x.name, dims)
    ncvar_def(ds, name, val, dimnames; attrib = attrib, deflatelevel = compress, kwargs...)
    # var = NcVar(varname, dims; t = type, compress = compress)
end


text:
# https://medium.com/coffee-in-a-klein-bottle/developing-your-julia-package-682c1d309507
using Random, Distributions
using Graphs
using Pkg, Revise
using CSV, DataFrames

Pkg.activate(".")
using RRSTExperiments

function prepare_graph(filename::String, α::Float64 = 1.0, β::Float64 = 3.0)
  # Read vertices, edges and first stage costs
  n, E, C = parse_graph_data(filename)

  # Assume second stage costs are equal to first stage costs
  c = copy(C)

  dist = Beta(α, β)
  d = generate_uncertain_costs(dist, c)

  return (n, E, C, c, d)
end

function experiment(seed::UInt32)
  # Generate uncertain costs
  Random.seed!(seed)

  # Prepare graph
  n, E, C, c, d = prepare_graph("data/london/london.gr")
  m = length(E)

  A = [InputEdge(a, b, C[i], c[i], d[i]) for (i, (a,b)) in enumerate(E)]
  ns = []; ks = []; cs = []; ms = []; nums = []; c1s = []; c2s = []; λs = []
  for k in [0, floor(Int64, 0.05 * m)]
    println("n=$n, k=$k: ")
    result1, x₁ = RRSTExperiments.solve_rec_st_with_algorithm(n, copy(A), k)

    for λ in [0.25, 0.75]
      print("    RRST = $result1, ")
      _, result2, x₂ = RRSTExperiments.solve_rec_st_hurwicz(n, copy(A), k, λ)
      println("Hurwicz = $result2.")

      println("    |symdiff(x₁, x₂)| = $(length(symdiff(x₁, x₂)))")
      println("    ∑Cx₁ = $(sum([e.C for e ∈ A if (e.i, e.j) in x₁ || (e.j, e.i) in x₁]))")
      println("    ∑Cx₂ = $(sum([e.C for e ∈ A if (e.i, e.j) in x₂ || (e.j, e.i) in x₂]))")

      for i in 1:5
        S = generate_scenario(Uniform(), A) # Generating actual scenario
        println("i=", i)
        _, c₁, _ = RRSTExperiments.solve_inc_st_with_model(n, S, A, x₁, k) # Recovery action for RRST
        _, c₂ , _ = RRSTExperiments.solve_inc_st_with_model(n, S, A, x₂, k) # Recovery action for Hurwicz

        println("    c₁=$c₁, c₂=$c₂")
        c1 = calculate_cost(A, x₁, c₁)
        c2 = calculate_cost(A, x₂, c₂)

        push!(ns, n)
        push!(ks, k)
        push!(cs, result1)
        push!(ms, result2)
        push!(λs, λ)
        push!(nums, i)
        push!(c1s, c1)
        push!(c2s, c2)
        println("$n,$k,$λ,$result1,$result2,$i,$c1,$c2")
      end
    end
  end

  # Write results
  filename = "data/london/london_output_$(randstring(4))_hur.csv"
  open(filename, "w")

  df = DataFrame("num_vertices"=>ns, "rec_param"=>ks, "λ"=>λs, "alg_sol_cost"=>cs, "hurwicz_sol_cost"=>ms, "experiment_num"=>nums, "alg_eval_cost"=>c1s, "hurwicz_eval_cost"=>c2s)
  println(df)

  CSV.write(filename, df)
end

experiment(0x00000001)


text:
using PackageCompiler
PackageCompiler.force_native_image!()
build_executable("./hello.jl","hello")

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using SuperStiff.PeriodizeSC
using SuperStiff.Stiffness

paramsfile = "params.json"
params = JSON.parsefile(paramsfile)
Grid_ = 100
M_tol = -0.01
path_ = params["data_loop"]
data = readdlm(path_, header=true)
relevant_data = hcat(data[1][:,find(x->x=="mu",data[2])],data[1][:,find(x->x=="ave_mu",data[2])])

block_params = Dict{String,Any}(
    "path_to_files" => params["path_to_files"],
    "data_loop" => params["data_loop"],
    "pattern" => params["pattern"], # Pattern of the files to be loaded
    "beta" => params["beta"],
    "print_mu_dop" => params["Print_mu_dop"], # Parameter used to create stiffness.dat file (file used with SE_G)
    "w_discretization" => params["w_discretization"], # Matsubara frequency grid
    "AFM_SC_NOCOEX" => params["AFM_SC_NOCOEX"],
    "cumulant" => params["cumulant"],
    "Periodization" => params["Periodization"],
    "fout_name" => params["fout_name"],
    "inplane_axis" => params["inplane_axis"]
    )

super_data_M, super_datap, list_of_files = Stiffness.read_data_loop(block_params,block_params["print_mu_dop"],block_params["pattern"])

#println(length(super_data_M),"\n",super_datap[1],"\n",list_of_files[1],"\n",datap)

list_modulevec = Stiffness.gen_modulevec_args(block_params,list_of_files,super_datap[:,1],params)
println("super_data_M: ", typeof(super_data_M))

match_data_loop = match(r"^(.*?)(?=/)",block_params["data_loop"]).match
println("match_data_loop: ", match_data_loop)

## Choosing the appropriate bare current vertices
if match_data_loop == "COEX" && block_params["Periodization"] == 0
    println("Tr.")
    if block_params["inplane_axis"] == "xx"
        fct_array = [PeriodizeSC.DxDxomegak, PeriodizeSC.DxDxzetakomegak, PeriodizeSC.DxDxomegakepsilonk, PeriodizeSC.DxDxzetakepsilonk, PeriodizeSC.DxDxzetak, PeriodizeSC.DxDxzetakepsilonk, PeriodizeSC.DxDxomegakepsilonk, PeriodizeSC.DxDxzetakepsilonk, PeriodizeSC.DxDxepsilonk]
    elseif block_params["inplane_axis"] == "yy"
        fct_array = [PeriodizeSC.DyDyomegak, PeriodizeSC.DyDyzetakomegak, PeriodizeSC.DyDyomegakepsilonk, PeriodizeSC.DyDyzetakepsilonk, PeriodizeSC.DyDyzetak, PeriodizeSC.DyDyzetakepsilonk, PeriodizeSC.DyDyomegakepsilonk, PeriodizeSC.DyDyzetakepsilonk, PeriodizeSC.DyDyepsilonk]
    end
elseif match_data_loop == "COEX" && block_params["Periodization"] == 1
    println("Per. or Cum.")
    if block_params["inplane_axis"] == "xx"
        fct_array = [PeriodizeSC.DxDxEpsilonbark, PeriodizeSC.DxDxZetakbar, PeriodizeSC.DxDxZetakbarepsilonk, PeriodizeSC.DxDxepsilonk]
    elseif block_params["inplane_axis"] == "yy"
        fct_array = [PeriodizeSC.DyDyEpsilonbark, PeriodizeSC.DyDyZetakbar, PeriodizeSC.DyDyZetakbarepsilonk, PeriodizeSC.DyDyepsilonk]
    end
elseif match_data_loop == "NOCOEX" && block_params["Periodization"] == 1 && block_params["AFM_SC_NOCOEX"] == 1
    println("AFM_SC_1")
    if block_params["inplane_axis"] == "xx"
        fct_array = [PeriodizeSC.DxDxZetakbar, PeriodizeSC.DxDxZetakbarepsilonk, PeriodizeSC.DxDxepsilonk]
    elseif block_params["inplane_axis"] == "yy"
        fct_array = [PeriodizeSC.DyDyZetakbar, PeriodizeSC.DyDyZetakbarepsilonk, PeriodizeSC.DyDyepsilonk]
    end
elseif match_data_loop == "NOCOEX" && (block_params["Periodization"] == 1 || block_params["Periodization"] == 0)
    if block_params["inplane_axis"] == "xx"
        println("NOCOEX along xx")
        fct_array = [PeriodizeSC.DxDxEpsilonbark]
    elseif block_params["inplane_axis"] == "yy"
        println("NOCOEX along yy")
        fct_array = [PeriodizeSC.DyDyEpsilonbark]
    end
end

println(fct_array)
k_array_prebuild = Array{Array{Float64,2},1}()
stiff_array = Array{Float64,1}()
tktilde_prebuild = Array{Array{Complex{Float64},2},2}(0,0)
@time for (l,modulevec_el) in enumerate(list_modulevec)
    if match_data_loop == "COEX"
        super_data_M_el = super_data_M[l]
    end
    w_discretization = modulevec_el.vals_["w_discretization"]
    beta = modulevec_el.vals_["beta"]
    zvec = 1.0im*[(2*n+1)*pi/beta for n in 1:w_discretization]
    modelvec = PeriodizeSC.ModelVector(Stiffness.ts(modulevec_el)[1],Stiffness.ts(modulevec_el)[2],Stiffness.ts(modulevec_el)[3],Stiffness.ts(modulevec_el)[4],zvec[1:w_discretization], modulevec_el.data_[1,1:w_discretization,:,:])
    if l <= 1
        for f in fct_array
            push!(k_array_prebuild, PeriodizeSC.k_grid(PeriodizeSC.Model(modelvec,l),Grid_,f))
        end
        tktilde_prebuild = PeriodizeSC.k_grid(PeriodizeSC.Model(modelvec,l),Grid_,PeriodizeSC.tktilde)
    end
    println("k_array_prebuild: ", typeof(k_array_prebuild))
    println("Iteration number: ",l)
    arr_w_kIntegrated = Matrix{Complex{Float64}}(0,0)
    if match_data_loop == "COEX" && block_params["Periodization"] == 0
        @time arr_w_kIntegrated = PeriodizeSC.stiffness_cluster_G_ab_k_grid(modelvec, k_array_prebuild, tktilde_prebuild, Grid_)
    elseif match_data_loop == "COEX" && block_params["Periodization"] == 1
        @time arr_w_kIntegrated = PeriodizeSC.stiffness_COEX_Per_Cum_ab_k_grid(modelvec, k_array_prebuild, tktilde_prebuild, Grid_, block_params["cumulant"], super_data_M_el, M_tol)
    elseif match_data_loop == "NOCOEX" && block_params["Periodization"] == 1
        println("NOCOEX Per. or Cum.")
        @time arr_w_kIntegrated = PeriodizeSC.stiffness_NOCOEX_Per_Cum_ab_k_grid(modelvec, k_array_prebuild, tktilde_prebuild, Grid_, block_params["cumulant"], block_params["AFM_SC_NOCOEX"])
    elseif match_data_loop == "NOCOEX" && block_params["Periodization"] == 0
        println("NOCOEX Tr.")
        @time arr_w_kIntegrated = PeriodizeSC.stiffness_NOCOEX_Tr_ab_k_grid(modelvec, k_array_prebuild, tktilde_prebuild, Grid_)
    end
    dop_stiff = 2.0/beta*sum(arr_w_kIntegrated[:,2])
    f = open(params["fout_name"], "a")
    if l == 1
        write(f, "#Grid "*"$(Grid_)"*"\n")
    end
    write(f, "$(relevant_data[:,1][l])"*"\t\t"*"$(relevant_data[:,2][l])"*"\t\t"*"$(real(dop_stiff))"*"\n")
    close(f)
    push!(stiff_array,dop_stiff)
    println(dop_stiff)
end
println("Superfluid Stiffness: ", stiff_array)

max_stars_repo_name: oliviersimard/SuperStiff

max_issues_repo_name: oliviersimard/SuperStiff

max_forks_repo_name: oliviersimard/SuperStiff


text:
module CacheStore
using ..SymbolServer: VarRef, FakeTypeName, FakeTypeofBottom, FakeTypeVar, FakeUnion, FakeUnionAll
using ..SymbolServer: ModuleStore, Package, FunctionStore, MethodStore, DataTypeStore, GenericStore
@static if !(Vararg isa Type)
    using ..SymbolServer: FakeTypeofVararg
end

const NothingHeader = 0x01
const SymbolHeader = 0x02
const CharHeader = 0x03
const IntegerHeader = 0x04
const StringHeader = 0x05
const VarRefHeader = 0x06
const FakeTypeNameHeader = 0x07
const FakeTypeofBottomHeader = 0x08
const FakeTypeVarHeader = 0x09
const FakeUnionHeader = 0x0a
const FakeUnionAllHeader = 0xb
const ModuleStoreHeader = 0x0c
const MethodStoreHeader = 0x0d
const FunctionStoreHeader = 0x0e
const DataTypeStoreHeader = 0x0f
const GenericStoreHeader = 0x10
const PackageHeader = 0x11
const TrueHeader = 0x12
const FalseHeader = 0x13
const TupleHeader = 0x14
const FakeTypeofVarargHeader = 0x15
const UndefHeader = 0x16

function write(io, x::VarRef)
    Base.write(io, VarRefHeader)
    write(io, x.parent)
    write(io, x.name)
end
function write(io, x::Nothing)
    Base.write(io, NothingHeader)
end
function write(io, x::Char)
    Base.write(io, CharHeader)
    Base.write(io, UInt32(x))
end
function write(io, x::Bool)
    x ? Base.write(io, TrueHeader) : Base.write(io, FalseHeader)
end
function write(io, x::Int)
    Base.write(io, IntegerHeader)
    Base.write(io, x)
end
function write(io, x::Symbol)
    Base.write(io, SymbolHeader)
    Base.write(io, sizeof(x))
    Base.write(io, String(x))
end
function write(io, x::NTuple{N,Any}) where N
    Base.write(io, TupleHeader)
    Base.write(io, N)
    for i = 1:N
        write(io, x[i])
    end
end
function write(io, x::String)
    Base.write(io, StringHeader)
    Base.write(io, sizeof(x))
    Base.write(io, x)
end
function write(io, x::FakeTypeName)
    Base.write(io, FakeTypeNameHeader)
    write(io, x.name)
    write_vector(io, x.parameters)
end
write(io, x::FakeTypeofBottom) = Base.write(io, FakeTypeofBottomHeader)
function write(io, x::FakeTypeVar)
    Base.write(io, FakeTypeVarHeader)
    write(io, x.name)
    write(io, x.lb)
    write(io, x.ub)
end
function write(io, x::FakeUnion)
    Base.write(io, FakeUnionHeader)
    write(io, x.a)
    write(io, x.b)
end
function write(io, x::FakeUnionAll)
    Base.write(io, FakeUnionAllHeader)
    write(io, x.var)
    write(io, x.body)
end

@static if !(Vararg isa Type)
    function write(io, x::FakeTypeofVararg)
        Base.write(io, FakeTypeofVarargHeader)
        isdefined(x, :T) ? write(io, x.T) : Base.write(io, UndefHeader)
        isdefined(x, :N) ? write(io, x.N) : Base.write(io, UndefHeader)
    end
end

function write(io, x::MethodStore)
    Base.write(io, MethodStoreHeader)
    write(io, x.name)
    write(io, x.mod)
    write(io, x.file)
    Base.write(io, x.line)
    Base.write(io, length(x.sig))
    for p in x.sig
        write(io, p[1])
        write(io, p[2])
    end
    write_vector(io, x.kws)
    write(io, x.rt)
end

function write(io, x::FunctionStore)
    Base.write(io, FunctionStoreHeader)
    write(io, x.name)
    write_vector(io, x.methods)
    write(io, x.doc)
    write(io, x.extends)
    write(io, x.exported)
end

function write(io, x::DataTypeStore)
    Base.write(io, DataTypeStoreHeader)
    write(io, x.name)
    write(io, x.super)
    write_vector(io, x.parameters)
    write_vector(io, x.types)
    write_vector(io, x.fieldnames)
    write_vector(io, x.methods)
    write(io, x.doc)
    write(io, x.exported)
end

function write(io, x::GenericStore)
    Base.write(io, GenericStoreHeader)
    write(io, x.name)
    write(io, x.typ)
    write(io, x.doc)
    write(io, x.exported)
end

function write(io, x::ModuleStore)
    Base.write(io, ModuleStoreHeader)
    write(io, x.name)
    Base.write(io, length(x.vals))
    for p in x.vals
        write(io, p[1])
        write(io, p[2])
    end
    write(io, x.doc)
    write(io, x.exported)
    write_vector(io, x.exportednames)
    write_vector(io, x.used_modules)
end

function write(io, x::Package)
    Base.write(io, PackageHeader)
    write(io, x.name)
    write(io, x.val)
    Base.write(io, UInt128(x.uuid))
    Base.write(io, x.sha === nothing ? zeros(UInt8, 32) : x.sha)
end

function write_vector(io, x)
    Base.write(io, length(x))
    for p in x
        write(io, p)
    end
end

function read(io, t = Base.read(io, UInt8))
    # There are a bunch of `yield`s in potentially expensive code paths.
    # One top-level `yield` would probably increase responsiveness in the
    # LS, but increases runtime by 3x. This seems like a good compromise.

    if t === VarRefHeader
        VarRef(read(io), read(io))
    elseif t === NothingHeader
        nothing
    elseif t === SymbolHeader
        n = Base.read(io, Int)
        out = Vector{UInt8}(undef, n)
        readbytes!(io, out, n)
        Symbol(String(out))
    elseif t === StringHeader
        yield()
        n = Base.read(io, Int)
        out = Vector{UInt8}(undef, n)
        readbytes!(io, out, n)
        String(out)
    elseif t === CharHeader
        Char(Base.read(io, UInt32))
    elseif t === IntegerHeader
        Base.read(io, Int)
    elseif t === FakeTypeNameHeader
        FakeTypeName(read(io), read_vector(io, Any))
    elseif t === FakeTypeofBottomHeader
        FakeTypeofBottom()
    elseif t === FakeTypeVarHeader
        FakeTypeVar(read(io), read(io), read(io))
    elseif t === FakeUnionHeader
        FakeUnion(read(io), read(io))
    elseif t === FakeUnionAllHeader
        FakeUnionAll(read(io), read(io))
    elseif t === FakeTypeofVarargHeader
        T, N = read(io), read(io)
        if T === nothing
            FakeTypeofVararg()
        elseif N === nothing
            FakeTypeofVararg(T)
        else
            FakeTypeofVararg(T, N)
        end
    elseif t === UndefHeader
        nothing
    elseif t === MethodStoreHeader
        yield()
        name = read(io)
        mod = read(io)
        file = read(io)
        line = Base.read(io, UInt32)
        nsig = Base.read(io, Int)
        sig = Vector{Pair{Any, Any}}(undef, nsig)
        for i in 1:nsig
            sig[i] = read(io) => read(io)
        end
        kws = read_vector(io, Symbol)
        rt = read(io)
        MethodStore(name, mod, file, line, sig, kws, rt)
    elseif t === FunctionStoreHeader
        yield()
        FunctionStore(read(io), read_vector(io, MethodStore), read(io), read(io), read(io))
    elseif t === DataTypeStoreHeader
        yield()
        DataTypeStore(read(io), read(io), read_vector(io, Any), read_vector(io, Any), read_vector(io, Any), read_vector(io, MethodStore), read(io), read(io))
    elseif t === GenericStoreHeader
        yield()
        GenericStore(read(io), read(io), read(io), read(io))
    elseif t === ModuleStoreHeader
        yield()
        name = read(io)
        n = Base.read(io, Int)
        vals = Dict{Symbol,Any}()
        sizehint!(vals, n)
        for _ = 1:n
            k = read(io)
            v = read(io)
            vals[k] = v
        end
        doc = read(io)
        exported = read(io)
        exportednames = read_vector(io, Symbol)
        used_modules = read_vector(io, Symbol)
        ModuleStore(name, vals, doc, exported, exportednames, used_modules)
    elseif t === TrueHeader
        true
    elseif t === FalseHeader
        false
    elseif t === TupleHeader
        N = Base.read(io, Int)
        ntuple(i->read(io), N)
    elseif t === PackageHeader
        yield()
        name = read(io)
        val = read(io)
        uuid = Base.UUID(Base.read(io, UInt128))
        sha = Base.read(io, 32)
        Package(name, val, uuid, all(x == 0x00 for x in sha) ? nothing : sha)
    else
        error("Unknown type: $t")
    end
end

function read_vector(io, T)
    n = Base.read(io, Int)
    v = Vector{T}(undef, n)
    for i in 1:n
        v[i] = read(io)
    end
    v
end

function storeunstore(x)
    io = IOBuffer()
    write(io, x)
    bs = take!(io)
    read(IOBuffer(bs))
end
end


text:
@testset "▶ types/legend                " begin
    # LEGEND
    l = G.Legend()
    @test isnothing(l.position)
    @test l.textstyle isa G.TextStyle
    @test l.offset == (0.0,0.0)
    @test isnothing(l.bgcolor)
    @test isnothing(l.margins)
    @test !l.nobox
    @test !l.off
end

@testset "▶ /legend                     " begin
    f = Figure()
    erase!(f)
    legend(position="top-left")
    @test f.axes[1].legend.position == "tl"
    legend(position="bottom-right")
    @test f.axes[1].legend.position == "br"
end

@testset "▶ set_prop/legend             " begin
    # Legend
    cla()
    line([0,0],[1,1],label="line")
    legend(nobox=true, margins=(1.0,1.0), offset=(.5,.5))
    @test gca().legend.nobox
    @test gca().legend.margins == (1.0,1.0)
    @test gca().legend.offset == (0.5,0.5)
end

@testset "▶ apply_gle/legend            " begin

end


text:
#const BG_COLOR = RGB(248/255, 247/255, 247/255)
# BG: 243 250 255 - LN: 67 75 86
#const BG_COLOR = RGB(243/255, 250/255, 255/255)
const BG_COLOR = RGB(255/255, 255/255, 255/255)
const LN_COLOR = RGB(67/255, 75/255, 86/255)
const MK_COLOR = RGB(67/255, 75/255, 86/255)
const LN01_COLOR = RGB(202/255,0/255,32/255)
const LN02_COLOR = RGB(5/255,113/255,176/255)
const FL01_COLOR = RGB(239/255, 138/255, 98/255)
const FL02_COLOR = RGB(103/255, 169/255, 207/255)

function plot_histogram(df::DataFrame, ycol::Symbol, output_dir::String)
    ENV["GKSwstype"] = "100"
    hist_path = output_dir * String(ycol) * "_total_hist.png"

    # plot histogram
    ht = Plots.histogram(float.(df[!, ycol]), title="Histogram of " * String(ycol),
        size = (2560, 1080),
        xlabel=String(ycol), ylabel="#", bins=200, legend=false,
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, fillcolor = FL02_COLOR)
    png(ht, hist_path)

end

function plot_lineplot_total(df::DataFrame, ycol::Symbol, output_dir::String)
    ENV["GKSwstype"] = "100"
    line_path = output_dir * String(ycol) * "_total_line.png"
    # plot in dates
    dates = DateTime.(df[!, :date], Local)
    pl = Plots.plot(dates, float.(df[!, ycol]),
        size = (2560, 1080),
        title=String(ycol) * " in dates", xlabel="date", ylabel=String(ycol),
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, linecolor = LN01_COLOR, legend=false)
    png(pl, line_path)
end

function plot_pcorr(_df::DataFrame, feas::AbstractArray, label_feas::AbstractArray, output_dir::String)
    ENV["GKSwstype"] = "100"

    crpl_path = output_dir * "pearson_corr"
    crcsv_path = output_dir * "pearson_corr.csv"

    dfM = convert(Matrix, _df[:, feas])

    dfm_cor = Statistics.cor(dfM)

    ann = []
    for i in 1:length(feas)
    for j in 1:length(feas)
        _ann = (i - 0.5, j - 0.5, Plots.text(Formatting.format(dfm_cor[i, j], precision=2), 18, :white))
        push!(ann, _ann)
    end
    end

    crpl = Plots.heatmap(string.(label_feas), string.(label_feas), float.(dfm_cor),
        size = (2560, 1080),
        clim = (-1.0, 1.0), c=:blues, legend=true, annotations = ann,
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        title="CORR", background_color = BG_COLOR)

    Plots.png(crpl, crpl_path * ".csv")
    Plots.svg(crpl, crpl_path * ".svg")

    df_cor = DataFrame(dfm_cor)
    DataFrames.rename!(df_cor, Symbol.(label_feas))
    CSV.write(crcsv_path, df_cor, writeheader = true)
end

function plot_corr(df_corr::DataFrame, output_size::Integer, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    plot_path = output_dir * "$(output_prefix)_corr_hourly"

    pl = Plots.plot(float.(df_corr[!, :hour]), float.(df_corr[!, :corr]),
        size = (2560, 1080),
        title="Correlation on hourly prediction", xlabel="hour", ylabel="corr",
        line=:solid, linewidth=5, label="OBS",
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, linecolor = LN01_COLOR, legend=false)

    png(pl, plot_path)

    nothing
end

function plot_corr_input(corr_arr::Array{F, 1}, title::String, output_dir::String, output_prefix::String) where F<:AbstractFloat
    ENV["GKSwstype"] = "100"

    line_path = output_dir * "$(output_prefix)_corr.png"

    pl = Plots.plot(float.(corr_arr),
        size = (2560, 1080),
        line=:solid, linewidth=5,
        ylim=(min(0.0, minimum(corr_arr)), 1.0),
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, color= :black,
        title = title,
        xlabel = "lag", ylabel = "corr", legend=false)
    plot!([72], color= :black, linewidth=5, line=:dash, seriestype="vline")

    Plots.png(pl, line_path)

    nothing
end

function plot_corr_OU(lag_arr::Array{I, 1}, corr_arr::Array{F, 1},
    title::String, output_dir::String, output_prefix::String) where {F<:AbstractFloat, I<:Integer}
    ENV["GKSwstype"] = "100"

    line_path = output_dir * "$(output_prefix).png"

    pl = Plots.plot(float.(lag_arr), float.(corr_arr),
        size = (2560, 1080),
        line=:solid, linewidth=5,
        ylim=(min(0.0, minimum(corr_arr)), 1.0),
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=30PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, color= :black,
        title = title,
        xlabel = "time", ylabel = "corr", legend=false)
    # plot 72 hour line to check correlation is zero at 72 hour.
    plot!([72], color= :black, linewidth=5, line=:dash, seriestype="vline")

    Plots.png(pl, line_path)

    nothing
end

function plot_corr_OU(corr_arr::Array{F, 1},
    title::String, output_dir::String, output_prefix::String) where {F<:AbstractFloat, I<:Integer}
    ENV["GKSwstype"] = "100"

    line_path = output_dir * "$(output_prefix).png"

    pl = Plots.plot(float.(corr_arr),
        size = (2560, 1080),
        line=:solid, linewidth=5,
        ylim=(min(0.0, minimum(corr_arr)), 1.0),
        guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=30PlotMeasures.px,
        guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
        background_color = BG_COLOR, color= :black,
        title = title,
        xlabel = "time", ylabel = "corr", legend=false)
    # plot 72 hour line to check correlation is zero at 72 hour.
    plot!([72], color= :black, linewidth=5, line=:dash, seriestype="vline")

    Plots.png(pl, line_path)

    nothing
end

function plot_DNN_scatter(dnn_table::Array{IndexedTable, 1}, ycol::Symbol,
    output_size::Integer, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        sc_path = output_dir * "$(i_pad)/" * "$(output_prefix)_scatter_$(i_pad)h.png"
        lim = max(maximum(JuliaDB.select(dnn_table[i], :y)),
            maximum(JuliaDB.select(dnn_table[i], :ŷ)))

        sc = Plots.scatter(float.(JuliaDB.select(dnn_table[i], :y)), float.(JuliaDB.select(dnn_table[i], :ŷ)),
            size = (1080, 1080),
            xlim = (0, lim), ylim = (0, lim), legend=false,
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            title="DNN/OBS ($(i_pad)h)", xlabel="OBS", ylabel="DNN",
            background_color = BG_COLOR, markercolor = MK_COLOR)

        # diagonal line (corr = 1)
        Plots.plot!(collect(0:0.1:lim), collect(0:0.1:lim),
            xlim = (0, lim), ylim = (0, lim), legend=false,
            background_color = BG_COLOR, linecolor = LN_COLOR)
        Plots.png(sc, sc_path)

        sc_csvpath = output_dir * "$(i_pad)/" * "$(output_prefix)_scatter_$(i_pad)h.csv"
        df = DataFrame(obs = float.(JuliaDB.select(dnn_table[i], :y)),
            model = float.(JuliaDB.select(dnn_table[i], :ŷ)))
        CSV.write(sc_csvpath, df, writeheader = true)
    end

    nothing
end

function plot_DNN_scatter(dfs::Array{DataFrame, 1}, ycol::Symbol,
    output_size::Integer, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        sc_path = output_dir * "$(i_pad)/" * "$(output_prefix)_scatter_$(i_pad)h.png"
        df = dfs[i]
        lim = max(maximum(float.(df[!, :y])), maximum(float.(df[!, :yhat])))

        sc = Plots.scatter(float.(df[!, :y]), float.(df[!, :yhat]),
            size = (1080, 1080),
            xlim = (0, lim), ylim = (0, lim), legend=false,
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            title="DNN/OBS ($(i_pad)h)", xlabel="OBS", ylabel="DNN",
            background_color = BG_COLOR, markercolor = MK_COLOR)

        # diagonal line (corr = 1)
        Plots.plot!(collect(0:0.1:lim), collect(0:0.1:lim),
            xlim = (0, lim), ylim = (0, lim), legend=false,
            background_color = BG_COLOR, linecolor = LN_COLOR)
        Plots.png(sc, sc_path)

        sc_csvpath = output_dir * "$(i_pad)/" * "$(output_prefix)_scatter_$(i_pad)h.csv"
        outdf = DataFrame(obs = float.(df[!, :y]), model = float.(df[!, :yhat]))
        CSV.write(sc_csvpath, outdf, writeheader = true)
    end

    nothing
end

function plot_DNN_histogram(dnn_table::Array{IndexedTable, 1}, ycol::Symbol,
    output_size::Integer, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        hs_OBS_path = output_dir * "$(i_pad)/" * "$(output_prefix)_hist(obs)_$(i_pad)h.png"
        hs_DNN_path = output_dir * "$(i_pad)/" * "$(output_prefix)_hist(dnn)_$(i_pad)h.png"

        ht = Plots.histogram(float.(JuliaDB.select(dnn_table[i], :y)), label="OBS",
            size = (2560, 1080),
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            title="Histogram of OBS ($(i_pad)h)", ylabel="#",
            background_color = BG_COLOR, fillcolor = FL02_COLOR)
        Plots.png(ht, hs_OBS_path)

        ht = Plots.histogram(float.(JuliaDB.select(dnn_table[i], :ŷ)), label="DNN",
            size = (2560, 1080),
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            title="Histogram of DNN ($(i_pad)h)", ylabel="#",
            background_color = BG_COLOR, fillcolor = FL02_COLOR)
        Plots.png(ht, hs_DNN_path)
    end

    nothing
end

function plot_DNN_histogram(dfs::Array{DataFrame, 1}, ycol::Symbol,
    output_size::Integer, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        hs_OBS_path = output_dir * "$(i_pad)/" * "$(output_prefix)_hist(obs)_$(i_pad)h.png"
        hs_DNN_path = output_dir * "$(i_pad)/" * "$(output_prefix)_hist(dnn)_$(i_pad)h.png"
        df = dfs[i]

        ht = Plots.histogram(float.(df[!, :y]), label="OBS",
            size = (2560, 1080),
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            title="Histogram of OBS ($(i_pad)h)", ylabel="#",
            background_color = BG_COLOR, fillcolor = FL02_COLOR)
        Plots.png(ht, hs_OBS_path)

        ht = Plots.histogram(float.(df[!, :yhat]), label="DNN",
            size = (2560, 1080),
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            title="Histogram of DNN ($(i_pad)h)", ylabel="#",
            background_color = BG_COLOR, fillcolor = FL02_COLOR)
        Plots.png(ht, hs_DNN_path)
    end

    nothing
end

"""
    plot_DNN_lineplot(dates, dnn_table, ycol, output_size, output_dir, output_prefix)

Plot prediction by dates (full length).
"""
function plot_DNN_lineplot(dates::Array{DateTime, 1}, dnn_table::Array{IndexedTable, 1}, ycol::Symbol,
    output_size::Integer, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        line_plotpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.png"
        dates_h = dates .+ Dates.Hour(i)

        pl = Plots.plot(dates_h, float.(JuliaDB.select(dnn_table[i], :y)),
            size = (2560, 1080),
            ylim = (0.0,
                max(maximum(JuliaDB.select(dnn_table[i], :y)), maximum(JuliaDB.select(dnn_table[i], :ŷ)))),
            line=:solid, linewidth=5, label="OBS",
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            background_color = BG_COLOR, color=LN01_COLOR,
            title=String(ycol) * " by dates ($(i_pad)h)",
            xlabel="date", ylabel=String(ycol), legend=:best)

        pl = Plots.plot!(dates_h, float.(JuliaDB.select(dnn_table[i], :ŷ)),
            line=:solid, linewidth=5, color=LN02_COLOR, label="DNN")
        Plots.png(pl, line_plotpath)

        line_csvpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.csv"
        df = DataFrame(dates = dates_h,
            obs = float.(JuliaDB.select(dnn_table[i], :y)),
            model = float.(JuliaDB.select(dnn_table[i], :ŷ)))
        CSV.write(line_csvpath, df, writeheader = true)
    end

    nothing
end

"""
    plot_DNN_lineplot(dates, dnn_table, s_date, f_date, ycol, output_size, output_dir, output_prefix)

Plot prediction by dates (given date range).
"""
function plot_DNN_lineplot(dates::Array{DateTime, 1}, dnn_table::Array{IndexedTable, 1},
    s_date::DateTime, f_date::DateTime, ycol::Symbol,
    output_size::Integer, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        line_plotpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.png"
        dates_h = dates .+ Dates.Hour(i)

        pl = Plots.plot(dates_h, float.(JuliaDB.select(dnn_table[i], :y)),
            size = (2560, 1080),
            ylim = (0.0,
                max(maximum(JuliaDB.select(dnn_table[i], :y)), maximum(JuliaDB.select(dnn_table[i], :ŷ)))),
            line=:solid, linewidth=5, label="OBS",
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            background_color = BG_COLOR, color=LN01_COLOR,
            title=String(ycol) * " by dates ($(i_pad)h)",
            xlabel="date", ylabel=String(ycol), legend=:best)

        pl = Plots.plot!(dates_h, float.(JuliaDB.select(dnn_table[i], :ŷ)),
            line=:solid, linewidth=5, color=LN02_COLOR, label="DNN")
        Plots.png(pl, line_plotpath)

        line_csvpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.csv"
        df = DataFrame(dates = dates_h,
            obs = float.(JuliaDB.select(dnn_table[i], :y)),
            model = float.(JuliaDB.select(dnn_table[i], :ŷ)))
        CSV.write(line_csvpath, df, writeheader = true)
    end

    nothing
end

"""
    plot_DNN_lineplot(dates, dnn_table, ycol, output_size, output_dir, output_prefix)

Plot prediction by dates (full length).
"""
function plot_DNN_lineplot(dfs::Array{DataFrame, 1}, ycol::Symbol,
    output_size::Integer, output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        line_plotpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.png"
        df = dfs[i]
        dates_h = df[!, :date] .+ Dates.Hour(i)

        @show first(df, 5)

        pl = Plots.plot(dates_h, float.(df[!, :y]),
            size = (2560, 1080),
            ylim = (0.0,
                max(maximum(float.(df[!, :y])), maximum(float.(df[!, :yhat])))),
            line=:solid, linewidth=5, label="OBS",
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            background_color = BG_COLOR, color=LN01_COLOR,
            title=String(ycol) * " by dates ($(i_pad)h)",
            xlabel="date", ylabel=String(ycol), legend=:best)

        pl = Plots.plot!(dates_h, float.(df[!, :yhat]),
            line=:solid, linewidth=5, color=LN02_COLOR, label="DNN")
        Plots.png(pl, line_plotpath)

        line_csvpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.csv"
        outdf = DataFrame(dates = dates_h, obs = float.(df[!, :y]), model = float.(df[!, :yhat]))
        CSV.write(line_csvpath, outdf, writeheader = true)
    end

    nothing
end

"""
    plot_RNN_lineplot(dates, dnn_table, ycol, test_fdate, test_tdate, output_size, output_dir, output_prefix)

Plot prediction by dates (full length).
"""
function plot_RNN_lineplot(dfs::Array{DataFrame, 1}, ycol::Symbol,
    output_size::Integer, test_fdate::DateTime, test_tdate::DateTime,
    output_dir::String, output_prefix::String)

    ENV["GKSwstype"] = "100"

    for i = 1:output_size
        i_pad = lpad(i, 2, '0')
        line_plotpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.png"
        # filter by date because batched input makes `undef` DateTime values and throw Overflow
        df = filter(row -> test_fdate <= row[:date] <= test_tdate, dfs[i])
        dates_h = df[!, :date] .+ Dates.Hour(i)

        pl = Plots.plot(dates_h, float.(df[!, :y]),
            size = (2560, 1080),
            ylim = (0.0,
                max(maximum(float.(df[!, :y])), maximum(float.(df[!, :yhat])))),
            line=:solid, linewidth=5, label="OBS",
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            background_color = BG_COLOR, color=LN01_COLOR,
            title=String(ycol) * " by dates ($(i_pad)h)",
            xlabel="date", ylabel=String(ycol), legend=:best)

        pl = Plots.plot!(dates_h, float.(df[!, :yhat]),
            line=:solid, linewidth=5, color=LN02_COLOR, label="DNN")
        Plots.png(pl, line_plotpath)

        line_csvpath = output_dir * "$(i_pad)/" * "$(output_prefix)_line_$(i_pad)h.csv"
        outdf = DataFrame(dates = dates_h, obs = float.(df[!, :y]), model = float.(df[!, :yhat]))
        CSV.write(line_csvpath, outdf, writeheader = true)
    end

    nothing
end

function plot_evaluation(df::DataFrame, ycol::Symbol, eval_syms, output_dir::String)
    ENV["GKSwstype"] = "100"

    last_epoch = df[end, :epoch]

    for eval_sym in eval_syms
        gr(size = (1920, 1080))

        plot_path = output_dir * String(ycol) * "_eval_$(String(eval_sym)).png"

        pl = Plots.plot(float.(df[:, :epoch]), float.(df[:, eval_sym]),
            size = (2560, 1080),
            color=:black, linewidth=6,
            guidefontsize = 18, titlefontsize = 24, tickfontsize = 18, legendfontsize = 18, margin=15PlotMeasures.px,
            guidefontcolor = LN_COLOR, titlefontcolor = LN_COLOR, tickfontcolor = LN_COLOR, legendfontcolor = LN_COLOR,
            background_color = BG_COLOR,
            title="$(String(eval_sym)) " * String(ycol),
            xlabel="epoch", ylabel="$(String(eval_sym))", legend=false)
        annotate!([(last_epoch, df[last_epoch, eval_sym], text("Value: " *  string(df[last_epoch, eval_sym]), 18, :black, :right))])
        Plots.png(pl, plot_path)
    end

    eval_path = output_dir * String(ycol) * "_evals.csv"
    CSV.write(eval_path, df, writeheader = true)
    nothing
end


text:
using SqState
using Flux

dir = "sqth"
device = SqState.get_device()
model_name = "q2args_sqth_fno"
batch_size = 25
epochs = 30
η₀ = 1e-3

m = SqState.q2args_sqth_fno() |> device
loss(𝐱, 𝐲) = Flux.mse(m(𝐱), 𝐲)
opt = Flux.Optimiser(WeightDecay(1e-4), Flux.ADAM(η₀))

# prepare data
data_files = readdir(joinpath(SqState.training_data_path(), dir), join=true)
loader4test = SqState.preprocess_q2args(data_files[1], batch_size=batch_size)
@info "numbers of data fragments: $(length(data_files)-1)"
data_loaders = Channel(5, spawn=true) do ch
    for e in 1:epochs, (i, data_file) in enumerate(data_files[2:end])
        put!(ch, SqState.preprocess_q2args(data_file, batch_size=batch_size))
        @info "Load epoch $e, file $i into buffer"
        flush(stdout)
    end
end

t = Ref{Int64}(1)
losses = Float32[]
data4validation = [(𝐱, 𝐲) for (𝐱, 𝐲) in loader4test] |> device
for loader4train in data_loaders
    @time begin
        data = [(𝐱, 𝐲) for (𝐱, 𝐲) in loader4train] |> device

        # descent η
        (t[] > 300) && (opt.os[2].eta = η₀ / 2^ceil((t[]-300)/300))

        # training
        Flux.train!(loss, params(m), data, opt)

        # collect loss
        validation_loss = sum(loss(𝐱, 𝐲) for (𝐱, 𝐲) in data4validation)/length(data4validation)
        in_data_loss = sum(loss(𝐱, 𝐲) for (𝐱, 𝐲) in data)/length(data)
        @info "$(t[])0k data\n η: $(opt.os[2].eta)\n in  loss: $in_data_loss\n out loss: $validation_loss"

        # update saved model
        push!(losses, validation_loss)
        (losses[end] == minimum(losses)) && SqState.update_model!(SqState.model_path(), model_name, m)

        # update indicator
        t[] += 1
    end
end


text:
################################################################################
#
#  final computation of the maximal order and automorphisms
#
################################################################################

function _simplify_components(L::Hecke.NfRelNS{nf_elem}, autL::Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem})
  if ngens(L) == 1
    return L, autL
  end
  to_simplify = Int[]
  for i = 1:length(L.pol)
    if total_degree(L.pol[i]) > 2
      push!(to_simplify, i)
    end
  end
  if isempty(to_simplify)
    return L, autL
  end
  pols = Vector{Generic.Poly{nf_elem}}(undef, ngens(L))
  maps = Vector{NfRelNSElem{nf_elem}}(undef, ngens(L))
  for i = 1:length(pols)
    Li, mLi = component(L, i)
    if !(i in to_simplify)
      pols[i] = Li.pol
      maps[i] = L[i]
      continue
    end
    Linew, mLinew = simplify(Li)
    pols[i] = Linew.pol
    maps[i] = mLi(mLinew.prim_img)
  end
  Lnew, gLnew = number_field(pols, cached = false, check = false)
  iso = hom(Lnew, L, maps)
  inv_iso = inv(iso)
  autsLnew = Vector{Hecke.NfRelNSToNfRelNSMor_nf_elem}(undef, length(autL))
  for i = 1:length(autL)
    autsLnew[i] = (iso*autL[i])*inv_iso
  end
  return Lnew, autsLnew
end

function _from_relative_to_abs_with_embedding(L1::Hecke.NfRelNS{nf_elem}, autL1::Array{Hecke.NfRelNSToNfRelNSMor_nf_elem, 1})
  #@time res = _from_relative_to_abs_with_embedding1(L1, autL1)
  res = _relative_to_absolute(L1, autL1)
  return res
end

function _from_relative_to_abs_with_embedding1(L1::Hecke.NfRelNS{nf_elem}, autL1::Array{Hecke.NfRelNSToNfRelNSMor_nf_elem, 1})
  @vtime :Fields 3 L, autL = _simplify_components(L1, autL1)
  S, mS = simple_extension(L, cached = false)
  K, mK, MK = absolute_field(S, cached = false)
  #First, we compute the maximal order of the absolute field.
  #We start from the maximal orders of the relative extension and of the base field.
  #FALSE: Since the computation of the relative maximal order is slow, I prefer to bring to the absolute field the elements
  # generating the equation order.
  pols = L.pol
  gL = gens(L)
  B = Array{nf_elem, 1}(undef, degree(K))
  B[1] = K(1)
  ind = total_degree(pols[1])
  genjj = mK\(mS\gL[1])
  for i = 2:ind
    B[i] = genjj*B[i-1]
  end
  for jj = 2:length(pols)
    genjj = mK\(mS\gL[jj])
    el = deepcopy(genjj)
    new_deg = total_degree(pols[jj])
    for i = 2:new_deg
      for j = 1:ind
        B[(i-1)* ind + j] = B[j]* el
      end
      mul!(el, el, genjj)
    end
    ind *= new_deg
  end

  #Now, I add the elements of the maximal order
  OB = maximal_order(base_field(S))
  for i = 1:degree(OB)
    el = MK(OB.basis_nf[i])
    for j = 1:ind
      B[(i-1)* ind + j] = B[j] * el
    end
  end
  @vprint :Fields 2 "Product basis computed\n"
  #Now, we compute the maximal order. Then we simplify.
  #We simplify only if the degree of the field is lower than 30

  BasisMat = basis_matrix(B, FakeFmpqMat)
  @vtime :Fields 3 Hecke.hnf_modular_eldiv!(BasisMat.num, BasisMat.den, :lowerleft)
  NewBMat = FakeFmpqMat(BasisMat.num, BasisMat.den)
  @vtime :Fields 3 Ostart = NfAbsOrd(K, NewBMat)
  Ostart.index = divexact(NewBMat.den^degree(K), prod_diagonal(NewBMat.num))
  Ostart.gen_index = fmpq(Ostart.index)
  Ostart.disc = divexact(numerator(discriminant(K)), Ostart.index^2)
  ram_primes_rel = numerator(norm(discriminant(L)))
  ram_primes_down = Hecke.ramified_primes(OB)
  for p in ram_primes_down
    if isone(gcd(p, ram_primes_rel))
      push!(Ostart.primesofmaximality, p)
    end
  end
  @vtime :Fields 3 O1 = MaximalOrder(Ostart)
  O1.ismaximal = 1
  Hecke._set_maximal_order_of_nf(K, O1)
  #simplify takes care of translating the order.
  @vtime :Fields 3 Ks, mKs = Hecke.simplify(K, cached = false)
  #I want also the embedding of the old field into the new one.
  #It is enough to find the image of the primitive element.
  k = base_field(S)
  a = MK(gen(k))
  embed = hom(k, Ks, mKs\a, check = false)
  #@assert iszero(k.pol(img_a))
  @vprint :Fields 3 "MaximalOrder Computed. Now Automorphisms\n"
  #Now, the automorphisms.
  # I need both generators and the whole group.
  autos = Array{NfToNfMor, 1}(undef, length(autL))
  el = mS(mK(image_primitive_element(mKs)))
  el1 = mS(mK(gen(K)))
  for i=1:length(autL)
    #@assert iszero(K.pol(mK(mS\(autL[i](el1)))))
    x = mKs\(mK\(mS\(autL[i](el))))
    #@assert Ks.pol(y) == 0
    autos[i] = hom(Ks, Ks, x, check = false)
  end
  @vprint :Fields 2 "Finished\n"
  #@assert codomain(embed) == Ks
  return Ks, autos, embed
end

###############################################################################
#
#  From relative to absolute
#
###############################################################################

function _relative_to_absoluteQQ(L::NfRelNS{nf_elem}, auts::Vector{NfRelNSToNfRelNSMor_nf_elem})
  K, gK = number_field(NfAbsNS, L)
  Ks, mKs = simplified_simple_extension1(K, isabelian = true, cached = false)
  #Now, I have to translate the automorphisms.
  #First, to automorphisms of K.
  autsK = Vector{NfAbsNSToNfAbsNS}(undef, length(auts))
  Qxy = parent(K.pol[1])
  for i = 1:length(auts)
    embs = Vector{NfAbsNSElem}(undef, ngens(K))
    imgs = image_generators(auts[i])
    for j = 1:length(imgs)
      embs[j] = K(map_coeffs(FlintQQ, imgs[j].data, parent = Qxy))
    end
    autsK[i] = hom(K, K, embs, check = false)
  end
  #Final step: translate the automorphisms to Ks
  Hecke._assure_has_inverse_data(mKs)
  autsKs = Vector{NfToNfMor}(undef, length(autsK))
  for i = 1:length(autsK)
    autsKs[i] = hom(Ks, Ks, mKs\(image_primitive_element(mKs*autsK[i])), check = false)
  end
  return Ks, autsKs
end

function _relative_to_absolute(L::NfRelNS{nf_elem}, auts::Vector{NfRelNSToNfRelNSMor_nf_elem})
  Ks, mKs = simplified_absolute_field(L)
  Hecke._assure_has_inverse_data(mKs)
  #Now, I have to translate the automorphisms.
  #First, to automorphisms of K.
  autsKs = Vector{NfToNfMor}(undef, length(auts))
  imggen = mKs(gen(Ks))
  for i = 1:length(auts)
    autsKs[i] = hom(Ks, Ks, mKs\(auts[i](imggen)))
  end
  #finally, the embedding of the base field of L into Ks
  k = base_field(L)
  img = mKs\(L(gen(k)))
  embed = hom(k, Ks, img, check = false)
  return Ks, autsKs, embed
end


text:
######################### Boundary Conditions #########################
# The job of the BC functions is to calculate/enforce face fluxes at boundary faces.
# The JST method only calculates fluxes at internal faces, then these conditions are applied to calculate them at the boundaries

# InletConditions: [ Static Pressure, Static Temperture, Ux, Uy, Uz ]
function supersonicInletBoundary(mesh::Mesh, sln::SolutionState, boundaryNumber, inletConditions, fluid::Fluid)
    P, T, Ux, Uy, Uz = inletConditions

    rho = idealGasRho(T, P, fluid.R)
    xMom, yMom, zMom = [Ux, Uy, Uz] .* rho
    e = calPerfectEnergy(T, fluid)
    eV2 = rho*(e + (mag([Ux, Uy, Uz])^2)/2)

    boundaryFluxes = Vector{Float64}(undef, 15)
    calculateFluxes3D!(boundaryFluxes, [P, T, Ux, Uy, Uz],  [rho, xMom, yMom, zMom, eV2])
    currentBoundary = mesh.boundaryFaces[boundaryNumber]
    @inbounds for face in currentBoundary
        sln.faceFluxes[face,:] = boundaryFluxes
    end
end

# InletConditions: [ totalPressure, totalTemp, nx, ny, nz ]
# Where n is the unit vector representing the direction of inlet velocity
# Using method from FUN3D solver
function subsonicInletBoundary(mesh, sln::SolutionState, boundaryNumber, inletConditions, fluid::Fluid)
    Pt, Tt, nx, ny, nz = inletConditions
    Cp, R, gamma = (fluid.Cp, fluid.R, fluid.gamma)

    boundaryFluxes = Vector{Float64}(undef, 15)
    primitives = Vector{Float64}(undef, 5)
    state = Vector{Float64}(undef, 5)
    velUnitVector = [ nx, ny, nz ]
    currentBoundary = mesh.boundaryFaces[boundaryNumber]

    @inbounds @fastmath for face in currentBoundary
        ownerCell = mesh.faces[face][1]
        @views adjustedVelocity = dot(velUnitVector, sln.cellPrimitives[ownerCell, 3:5])
        primitives[2] = Tt - (gamma-1)/2 * (adjustedVelocity^2)/(gamma*R)
        machNum = abs(adjustedVelocity) / sqrt(gamma * R * primitives[2])
        primitives[1] = Pt*(1 + (gamma-1)/2 * machNum^2)^(-gamma/(gamma-1))
        primitives[3:5] = adjustedVelocity .* velUnitVector

        # Calculate state variables
        state[1] = idealGasRho(primitives[2], primitives[1], R)
        state[2:4] .= primitives[3:5] .* state[1]
        e = calPerfectEnergy(primitives[2], fluid)
        state[5] = state[1]*(e + (mag(primitives[3:5])^2)/2 )

        calculateFluxes3D!(boundaryFluxes, primitives, state)
        sln.faceFluxes[face, :] = boundaryFluxes
    end
end

function pressureOutletBoundary(mesh, sln::SolutionState, boundaryNumber, outletPressure, fluid::Fluid)
    nFluxes = size(sln.cellFluxes, 2)

    # Directly extrapolate cell center flux to boundary (zero gradient between the cell center and the boundary)
    currentBoundary = mesh.boundaryFaces[boundaryNumber]
    @inbounds @fastmath for face in currentBoundary
        ownerCell = mesh.faces[face][1]
        # Mass fluxes are unaffected by pressure boundary
        for flux in 1:nFluxes
            sln.faceFluxes[face, flux] = sln.cellFluxes[ownerCell, flux]
        end

        origP = sln.cellPrimitives[ownerCell, 1]
        # Adjust the momentum fluxes containing pressure
        for flux in [4, 8, 12]
            sln.faceFluxes[face, flux] += outletPressure - origP
        end
        #Adjust the energy fluxes
        for d in 1:3
            sln.faceFluxes[face, 12+d] += sln.cellPrimitives[ownerCell, 2+d]*(outletPressure - origP)
        end
    end
end

function zeroGradientBoundary(mesh::Mesh, sln::SolutionState, boundaryNumber, _, __)
    nFluxes = size(sln.cellFluxes, 2)

    # Directly extrapolate cell center flux to boundary (zero gradient between the cell center and the boundary)
    currentBoundary = mesh.boundaryFaces[boundaryNumber]
    @inbounds for face in currentBoundary
        ownerCell = mesh.faces[face][1] #One of these will be -1 (no cell), the other is the boundary cell we want
        for flux in 1:nFluxes
            sln.faceFluxes[face, flux] = sln.cellFluxes[ownerCell, flux]
        end
    end
end

function wallBoundary(mesh::Mesh, sln::SolutionState, boundaryNumber, _, __)
    currentBoundary = mesh.boundaryFaces[boundaryNumber]
    @inbounds for f in currentBoundary
        ownerCell = max(mesh.faces[f][1], mesh.faces[f][2]) #One of these will be -1 (no cell), the other is the boundary cell we want

        faceP = sln.cellPrimitives[ownerCell, 1]
        # Momentum flux is Pressure in each of the normal directions (dot product)
        sln.faceFluxes[f, 4] = faceP
        sln.faceFluxes[f, 8] = faceP
        sln.faceFluxes[f, 12] = faceP

        # Mass Flux is zero
        sln.faceFluxes[f, 1:3] .= 0.0
        # Energy Flux is zero
        sln.faceFluxes[f, 13:15] .= 0.0
    end
end

# Symmetry is identical to a wall in Euler simulations
symmetryBoundary = wallBoundary

function emptyBoundary(mesh::Mesh, sln::SolutionState, boundaryNumber, _, __)
    return
end

function applyBoundaryConditions(mesh::Mesh, sln::SolutionState, boundaryConditions, nBoundaries, fluid::Fluid)
    for boundaryNumber in 1:nBoundaries
        bFunctionIndex = 2*boundaryNumber-1
        boundaryConditionFunction = boundaryConditions[bFunctionIndex]
        boundaryParameters = boundaryConditions[bFunctionIndex+1]

        boundaryConditionFunction(mesh, sln, boundaryNumber, boundaryParameters, fluid)
    end
end

max_stars_repo_name: henrystoldt/fvCFD


text:
# Grid search algorithm.
# Reference:  https://github.com/cstjean/ScikitLearn.jl/blob/master/src/grid_search.jl

export ParameterGrid, recursive_index, UniformGrid, LogGrid

"""Grid search iterator for parameters. The sequence of ranges defines the grids to search.

Parameters
----------
- `ParameterGrid`: Contain list of ranges. The range is in `(start, end, points)` order.

Examples
--------
```julia-repl
julia> ranges = [ (1.,10.,10.), (1.,10.,10.) ] # list of ranges (start_num, stop_num, number of grids`{int}`)
julia> param_range = ParameterRange(ranges)
```

```julia-repl

julia> ParameterGrid([ [1,1000,5], [1,3,1]]; method=LogGrid())
5-element ParameterGrid:
 [1.0999, 3.0]
 [1.999, 3.0]
 [10.99, 3.0]
 [100.9, 3.0]
 [1000.0, 3.0]
```

"""
@with_kw struct ParameterGrid <: ParameterGenerator
    param_ranges
    len
    indexes
    method =UniformGrid()
end

"""
Construct parameter grid from list of ranges (`[start, end, grid_num]`). The grid distribution is default to be uniform.
"""
function ParameterGrid(param_ranges; method= UniformGrid() :: GridSampler, SHUFFLE=false )
    len = mul(Int,[ i[end] for i in param_ranges ])
    indexes = Array(1:1:len)
    SHUFFLE ? shuffle!(indexes) : nothing # indexes to be sampled, shuffled
    return ParameterGrid(param_ranges, len, indexes, method);
end

function ParameterGrid(param_ranges::AbstractArray{T,1}; method=UniformGrid()::GridSampler, SHUFFLE=false) where T<:AbstractRange

    param_ranges_ = ranges2tuples(param_ranges)

    return ParameterGrid(param_ranges_; method=method, SHUFFLE=SHUFFLE)
end

"""
Grid object with ranges and distribution function

"""
function ParameterGrid(param_ranges, method)
    len = mul(Int,[ i[end] for i in param_ranges ])
    indexes = Array(1:1:len)
    shuffle!(indexes) # indexes to be sampled, shuffledP
    return ParameterGrid(param_ranges, len, indexes, method);
end

function Base.length(self::ParameterGrid)
    return self.len
end

Base.size(self::ParameterGrid) = (Base.length(self),)

## Iterator
"""Get the parameters that would be ```ind``` th in iteration

The total sample number is # to do

Paremeters
----------
- `ind` {int}: The iteration index

Returns
-------
- `params`: list of numbers in the range of `ParameterGrid.param_ranges`
"""
function Base.getindex(self::ParameterGrid, ind::Int)
    ind_ = self.indexes[ind]
    vec_i = recursive_index(self, ind_) # number in grid sample space

    values = zeros(length(vec_i))
    for i in eachindex(values)
        rang = self.param_ranges[i]
        str,ed, grid_num = rang
        values[i] = self.method(str,ed, grid_num,vec_i[i])
    end
    return values
end

"""
get the number of a given range from index with uniform distribution.

```julia-repl
julia> FindSteadyStates.uniformGrid(1,10,3, 2)
5
```
"""
struct UniformGrid <: GridSampler end

function (self::UniformGrid)(str_num, end_num, grid_num, ind)
    @assert end_num >= str_num  # end number is equal or bigger than start
    @assert ind <= grid_num # grid size should greater or equal to the specified index
    move = (end_num - str_num) / (grid_num-1)  # Step size. Start and end is included. Therefore: grid_num -1
    return str_num + move * (ind - 1.)
end

@with_kw struct LogGrid <: GridSampler
    decay = 10.
end

function (self::LogGrid)(str_num, end_num, grid_num, ind)
    @assert end_num > str_num  # end number is bigger than start
    @assert ind <= grid_num # grid size should greater or equal to the specified index
    step = (end_num - str_num)
    move = step/(self.decay)^(grid_num-ind)
    return str_num + move
end

"""
Recursive indexing.
This function uses list of ranges and index to provide a vector of indexes.

Example
-------
```julia-repl
julia> ranges = [(1,2,2), (4,45,2)]
julia> param_range = ParameterGrid(ranges)
julia> vec_i = recursive_index(param_range, 4)
julia> vec_i
[2,2]
```
"""
function recursive_index(self::ParameterGrid, ind) :: Array{Integer,1}

    lens = [ i[end] for i in self.param_ranges]

    index = self.indexes[ind]

    vec_param = zeros(length(self.param_ranges)) #output vector

    for i in eachindex(vec_param)
        if i == 1 # first index
            id = ind % lens[i] # mod
            if id == 0 # end of Array
                id = lens[i]
            end
            vec_param[i] = id
        else # otherwise
            l = mul(lens[1:i-1]) # multiply from first to the previous neighbor digit
            mov = (ind-1) ÷ l  # move 1 digit when the prevous layer exceed length
            id = mov % lens[i] + 1
            vec_param[i] = id
        end
    end

    return vec_param

end