text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля

#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки поля промакированы
function mark!(r)
    back!(r)
    putmarker!(r)
    while (!isborder(r,Sud) || !isborder(r,Nord)) && !isborder(r,Ost)
        for side in (Nord, Ost, Sud)
            if side == Ost
                putmarker!(r)
                move!(r,Ost)
            else
                while !isborder(r,side)
                    putmarker!(r)
                    move!(r,side)
                end
                if side == Sud && !isborder(r,Ost)
                    putmarker!(r)
                    move!(r,Ost)
                end
            end
        end
    end
    putmarker!(r)
end
function back!(r)
    for side in (West,Sud)
        while !isborder(r,side)
            move!(r,side)
        end
    end
end

ext: jl


text:
using Test
using Carrot

@testset "post_forex_prediction" begin
    @test create_request_text("chat_id", "the message") == "chat_id=chat_id&text=the%20message"
    @test create_request_text("chat_id", "message", message_id="message_id") == "chat_id=chat_id&text=message&reply_to_message_id=message_id"
end


text:
function wasserstein_distance(u_values::AbstractVector{T1}, v_values::AbstractVector{T2}, p::Integer) where { T1<:Real, T2<:Real}
  # code adapted from SciPy.stats._cdf_distance
  u_sorter = issorted(u_values) ? (1:length(u_values)) : sortperm(u_values)
  v_sorter = issorted(v_values) ? (1:length(v_values)) : sortperm(v_values)
  all_values = sort(vcat(u_values,v_values),alg=MergeSort)
  deltas = all_values[2:end] .- all_values[1:end-1]
  u_cdf = map(x->(searchsortedlast(u_values[u_sorter],x))/ length(u_values),all_values[1:end-1])
  v_cdf = map(x->(searchsortedlast(v_values[v_sorter],x))/ length(v_values),all_values[1:end-1])
  #println("deltas: ",deltas)
  #println("u_cdf:  ",u_cdf)
  #println("v_cdf:  ",u_cdf)
  #return nothing
  if p==1
     return sum(abs.(u_cdf.-v_cdf).*deltas)
  elseif p==2
      return sqrt(sum((u_cdf.-v_cdf).^2 .* deltas))
  else
      return (sum(abs.(u_cdf.-v_cdf).^p .* deltas))^(1/p)
  end
end

function wasserstein_distance(u_values::AbstractVector{T1}, v_values::AbstractVector{T2},
                              u_weights::AbstractVector{T3}, v_weights::AbstractVector{T4}, p::Integer
            ) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real }
  # code adapted from SciPy.stats._cdf_distance
  @assert length(u_values) == length(u_weights)
  @assert length(v_values) == length(v_weights)
  u_sorter = issorted(u_values) ? (1:length(u_values)) : sortperm(u_values)
  v_sorter = issorted(v_values) ? (1:length(v_values)) : sortperm(v_values)
  all_values = sort(vcat(u_values,v_values),alg=MergeSort)
  deltas = all_values[2:end] .- all_values[1:end-1]
  u_cdf_indices = map(x->searchsortedlast(u_values[u_sorter],x),all_values[1:end-1]) .+1
  v_cdf_indices = map(x->searchsortedlast(v_values[v_sorter],x),all_values[1:end-1]) .+1
  u_sorted_weights = vcat(0.,cumsum(u_weights[u_sorter]))
  v_sorted_weights = vcat(0.,cumsum(v_weights[v_sorter]))
  u_cdf_max = u_sorted_weights[end]
  v_cdf_max = v_sorted_weights[end]
  u_cdf = u_sorted_weights[u_cdf_indices] / u_cdf_max
  v_cdf = v_sorted_weights[v_cdf_indices] / v_cdf_max
  if p==1
     return sum(abs.(u_cdf.-v_cdf).*deltas)
    elseif p==2
      return sqrt(sum((u_cdf.-v_cdf).^2 .* deltas))
    else
      return (sum(abs.(u_cdf.-v_cdf).^p .* deltas))^(1/p)
    end
end

function wasserstein_distance_presorted(u_values::AbstractVector{T1}, v_values::AbstractVector{T2},
                              u_weights::AbstractVector{T3}, v_weights::AbstractVector{T4}, p::Integer
            ) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real }
  # code adapted from SciPy.stats._cdf_distance
  @assert length(u_values) == length(u_weights)
  @assert length(v_values) == length(v_weights)
  @assert issorted(u_values) && issorted(v_values)
  all_values = sort(vcat(u_values,v_values),alg=MergeSort)
  deltas = all_values[2:end] .- all_values[1:end-1]
  u_cdf_indices = map(x->searchsortedlast(u_values,x),all_values[1:end-1]) .+1
  v_cdf_indices = map(x->searchsortedlast(v_values,x),all_values[1:end-1]) .+1
  #return u_cdf_indices, v_cdf_indices, deltas, all_values
  u_sorted_weights = vcat(0.,cumsum(u_weights))
  v_sorted_weights = vcat(0.,cumsum(v_weights))
  u_cdf_max = u_sorted_weights[end]
  v_cdf_max = v_sorted_weights[end]
  u_cdf = u_sorted_weights[u_cdf_indices] / u_cdf_max
  v_cdf = v_sorted_weights[v_cdf_indices] / v_cdf_max
  if p==1
     return sum(abs.(u_cdf.-v_cdf).*deltas)
    elseif p==2
      return sqrt(sum((u_cdf.-v_cdf).^2 .* deltas))
    else
      return (sum(abs.(u_cdf.-v_cdf).^p .* deltas))^(1/p)
    end
end

function wasserstein_distance_presorted(u_values::AbstractVector{T1}, v_values::AbstractVector{T2},
                              u_weights::AbstractVector{T3}, v_weights::AbstractVector{T4}, p_list::AbstractVector{T5}
            ) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real }
  # code adapted from SciPy.stats._cdf_distance
  @assert length(u_values) == length(u_weights)
  @assert length(v_values) == length(v_weights)
  @assert issorted(u_values) && issorted(v_values)
  all_values = sort(vcat(u_values,v_values),alg=MergeSort)
  deltas = all_values[2:end] .- all_values[1:end-1]
  u_cdf_indices = map(x->searchsortedlast(u_values,x),all_values[1:end-1]) .+1
  v_cdf_indices = map(x->searchsortedlast(v_values,x),all_values[1:end-1]) .+1
  #return u_cdf_indices, v_cdf_indices, deltas, all_values
  u_sorted_weights = vcat(0.,cumsum(u_weights))
  v_sorted_weights = vcat(0.,cumsum(v_weights))
  u_cdf_max = u_sorted_weights[end]
  v_cdf_max = v_sorted_weights[end]
  u_cdf = u_sorted_weights[u_cdf_indices] / u_cdf_max
  v_cdf = v_sorted_weights[v_cdf_indices] / v_cdf_max
  map(p->(sum(abs.(u_cdf.-v_cdf).^p .* deltas))^(1/p), p_list)
end

function wasserstein_distance_presorted_common_x(u_values::AbstractVector{T1}, #v_values::AbstractVector{T2},
                              u_weights::AbstractVector{T3}, v_weights::AbstractVector{T4}, p::Integer
            ) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real }
  # code adapted from SciPy.stats._cdf_distance
  @assert length(u_values) == length(u_weights)
  @assert length(u_values) == length(v_weights)
  @assert issorted(u_values) #&& issorted(v_values)
  all_values = zeros(T1,2*length(u_values))
  all_values[1:2:end-1] .= u_values
  all_values[2:2:end] .= u_values
  #deltas = all_values[2:end] .- all_values[1:end-1]
  deltas = zeros(T1,2*length(u_values)-1)
  deltas[1:2:end] .= zero(T1)
  deltas[2:2:end-1] .= u_values[2:end].-u_values[1:end-1]
  u_cdf_indices = map(x->searchsortedlast(u_values,x),all_values[1:end-1]) .+1
  #v_cdf_indices = map(x->searchsortedlast(v_values,x),all_values[1:end-1]) .+1
  u_sorted_weights = vcat(0.,cumsum(u_weights))
  v_sorted_weights = vcat(0.,cumsum(v_weights))
  u_cdf_max = u_sorted_weights[end]
  v_cdf_max = v_sorted_weights[end]
  u_cdf = u_sorted_weights[u_cdf_indices] / u_cdf_max
  v_cdf = v_sorted_weights[u_cdf_indices] / v_cdf_max
  if p==1
     return sum(abs.(u_cdf.-v_cdf).*deltas)
    elseif p==2
      return sqrt(sum((u_cdf.-v_cdf).^2 .* deltas))
    else
      return sum(abs.(u_cdf.-v_cdf).^p .* deltas)^(1/p)
    end
end

function wasserstein_distance_presorted_common_x(u_values::AbstractVector{T1}, #v_values::AbstractVector{T2},
                              u_weights::AbstractVector{T3}, v_weights::AbstractVector{T4}, p_list::AbstractVector{T5}
            ) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real }
  # code adapted from SciPy.stats._cdf_distance
  @assert length(u_values) == length(u_weights)
  @assert length(u_values) == length(v_weights)
  @assert issorted(u_values) #&& issorted(v_values)
  all_values = zeros(T1,2*length(u_values))
  all_values[1:2:end-1] .= u_values
  all_values[2:2:end] .= u_values
  #deltas = all_values[2:end] .- all_values[1:end-1]
  deltas = zeros(T1,2*length(u_values)-1)
  deltas[1:2:end] .= zero(T1)
  deltas[2:2:end-1] .= u_values[2:end].-u_values[1:end-1]
  u_cdf_indices = map(x->searchsortedlast(u_values,x),all_values[1:end-1]) .+1
  #v_cdf_indices = map(x->searchsortedlast(v_values,x),all_values[1:end-1]) .+1
  u_sorted_weights = vcat(0.,cumsum(u_weights))
  v_sorted_weights = vcat(0.,cumsum(v_weights))
  u_cdf_max = u_sorted_weights[end]
  v_cdf_max = v_sorted_weights[end]
  u_cdf = u_sorted_weights[u_cdf_indices] / u_cdf_max
  v_cdf = v_sorted_weights[u_cdf_indices] / v_cdf_max
  map(p->(sum(abs.(u_cdf.-v_cdf).^p .* deltas))^(1/p), p_list)
end

#earth_mover_distance(u_values::AbstractVector{T1}, v_values::AbstractVector{T2}) where { T1<:Real, T2<:Real} = wasserstein_distance(u_values,v_values,1)
earth_mover_distance(u_values::AbstractVector{T1}, v_values::AbstractVector{T2}, u_weights::AbstractVector{T3}, v_weights::AbstractVector{T4}) where { T1<:Real, T2<:Real, T3<:Real, T4<:Real } = wasserstein_distance(u_values,v_values,u_weights,v_weights,1)
earth_mover_distance(u_values::AbstractVector{T1}, u_weights::AbstractVector{T3}, v_weights::AbstractVector{T4} ) where { T1<:Real, T3<:Real, T4<:Real } = wasserstein_distance_presorted_common_x(u_values,u_weights,v_weights,1)


text:
using Turing
using Turing.RandomMeasures
using DynamicPPL
using Random
using AutoGibbs
using ProgressMeter
using Dates
using InteractiveUtils

include("models.jl")

# parameters for chain replication
const CHAIN_LENGTH = 10_000   # sampling steps
const HMC_LF_SIZE = 0.05     # parameter 1 for HMC
const HMC_N_STEP = 10        # parameter 2 for HMC
const DATA_RNG = MersenneTwister(424242)

# experimental parameters
const DATA_SIZES = (10, 25, 50, 100)
const N_PARTICLES = (100,)
const BENCHMARK_COMPILATIONS = 5
const BENCHMARK_CHAINS = 10

function run_experiments(
    modelname,
    max_data_size,
    setup,
    compilation_times_channel,
    chains_channel
)
    generate, example, tarray_example, p_discrete, p_continuous = setup

    old_prog = Turing.PROGRESS[]
    Turing.turnprogress(false)

    GC.gc()

    for data_size in DATA_SIZES
        if data_size > max_data_size
            continue
            @info "[$modelname] Skipped data size $data_size > $max_data_size"
        end

        dataname, data = generate(DATA_RNG, data_size)
        model_ag = example(; dataname => data)
        model_pg = tarray_example(; dataname => data)

        # This is done outside the combination loop, so that we have several compile time
        # samples .The last one is reused for the actual sampler.
        local static_conditional
        @info "[$modelname] Compiling conditionals for data size $data_size"
        @showprogress for c in 1:BENCHMARK_COMPILATIONS
            start_time = time_ns()
            static_conditional = StaticConditional(model_ag, p_discrete)
            compilation_time = (time_ns() - start_time) / 1e9

            put!(compilation_times_channel,
                 (model = modelname,
                  data_size = data_size,
                  repetition = c,
                  compilation_time = compilation_time))
        end

        for (i, particles) in enumerate(N_PARTICLES)
            # mh = MH(p_continuous...)
            hmc = HMC(HMC_LF_SIZE, HMC_N_STEP, p_continuous...)
            pg = PG(particles, p_discrete)

            combinations = Dict([
                # ("AG", "MH") => (model_ag, Gibbs(static_conditional, mh)),
                ("AG", "HMC") => (model_ag, Gibbs(static_conditional, hmc)),
                # ("PG", "MH") => (model_pg, Gibbs(pg, mh)),
                ("PG", "HMC") => (model_pg, Gibbs(pg, hmc))
            ])

            GC.gc()

            for ((d_alg, c_alg), (model, sampler)) in combinations
                C = BENCHMARK_CHAINS

                @info "[$modelname] Sampling $C chains using $d_alg+$c_alg with data size $data_size and $particles particles"
                @showprogress for repetition in 1:C
                    start_time = time_ns()
                    chain = sample(model, sampler, CHAIN_LENGTH)
                    sampling_time = (time_ns() - start_time) / 1e9

                    put!(chains_channel,
                         (model = modelname,
                          discrete_algorithm = d_alg,
                          continuous_algorithm = c_alg,
                          particles = particles,
                          data_size = data_size,
                          repetition = repetition,
                          sampling_time = sampling_time,
                          chain = chain,
                          ))

                    # GC.gc()
                end
            end
        end
    end

    Turing.turnprogress(old_prog)
end

VARIABLE_NAMES = (
    "model",
    "discrete_algorithm",
    "continuous_algorithm",
    "particles",
    "data_size",
    "repetition"
)

function print_csv_line(io, values...)
    join(io, ("\"$(escape_string(string(v)))\"" for v in values), ",")
    # withoug flushing, this might not be written when the buffer is small and the process is killed
    flush(io)
end

function serialize_chains(filename_times, filename_diagnostics, filename_chains, observations)
    open(filename_times, "w") do f_times
        open(filename_diagnostics, "w") do f_diag
            open(filename_chains, "w") do f_chains
                # write header lines
                print_csv_line(f_times, VARIABLE_NAMES..., "sampling_time")
                print_csv_line(f_diag, VARIABLE_NAMES..., "parameter", "diagnostic", "value")
                print_csv_line(f_chains, VARIABLE_NAMES..., "parameter", "step", "value")

                for obs in observations
                    condition_data = Tuple(obs)[1:end-2]
                    param_names = names(obs.chain, :parameters)

                    # write line for the sampling time to `f_times`
                    print(f_times, '\n')
                    print_csv_line(f_times, condition_data..., obs.sampling_time)

                    # write line for parameter and diagnostic per chain to `f_diag`
                    diagnostics = ess(obs.chain)
                    for p in param_names, d in (:ess, :r_hat)
                        print(f_diag, '\n')
                        print_csv_line(f_diag, condition_data..., p, d, diagnostics[p, d])
                    end

                    # write line for each parameter and value in the chain to `f_chains`
                    for (step, sampling) in enumerate(eachrow(Array(obs.chain)))
                        for (p, value) in zip(param_names, sampling)
                            print(f_chains, '\n')
                            print_csv_line(f_chains, condition_data..., p, step, value)
                        end
                    end
                end
            end
        end
    end
end

function serialize_compilation_times(filename, observations)
    open(filename, "w") do f
        print_csv_line(f, "model", "data_size", "repetition", "compilation_time")

        for obs in observations
            print(f, '\n')
            print_csv_line(f, obs...)
        end
    end
end

const MODEL_SETUPS = Dict([
    "GMM" => (gmm_generate, gmm_example, gmm_tarray_example, :z, (:w, :μ)),
    "HMM" => (hmm_generate, hmm_example, hmm_tarray_example, :s, (:T, :m)),
    "IMM" => (imm_stick_generate, imm_stick_example, imm_stick_tarray_example, :z, (:μ, :v))
])

function main(
    modelname,
    max_data_size=maximum(DATA_SIZES),
    results_path=nothing,
    symlinks=false
)
    @info "[$modelname] Executed on:" version=sprint(InteractiveUtils.versioninfo)

    modelname = uppercase(modelname)
    max_data_size = parse(Int, max_data_size)

    if isnothing(results_path)
        results_path = "."
    end

    timestamp = round(Dates.now(), Dates.Minute)
    samplingtimes_fn = abspath(results_path, "$modelname-sampling_times-$timestamp.csv")
    diagnostics_fn = abspath(results_path, "$modelname-diagnostics-$timestamp.csv")
    chains_fn = abspath(results_path, "$modelname-chains-$timestamp.csv")
    compiletimes_fn = abspath(results_path, "$modelname-compile_times-$timestamp.csv")

    if haskey(MODEL_SETUPS, modelname)
        compilation_times_channel = Channel()
        chains_channel = Channel()

        @async try
            run_experiments(modelname,
                            max_data_size,
                            MODEL_SETUPS[modelname],
                            compilation_times_channel,
                            chains_channel)
        catch e
            @error "[$modelname] An error occured:" exception=(e, catch_backtrace())
        finally
            close(compilation_times_channel)
            close(chains_channel)
        end

        @sync begin
            @async serialize_chains(samplingtimes_fn, diagnostics_fn, chains_fn, chains_channel)
            @async serialize_compilation_times(compiletimes_fn, compilation_times_channel)
        end

        if symlinks
            # overwrite a symlink to the latest version, without time stamp
            samplingtimes_fn_nodate = abspath(results_path, "$modelname-sampling_times.csv")
            diagnostics_fn_nodate = abspath(results_path, "$modelname-diagnostics.csv")
            chains_fn_nodate = abspath(results_path, "$modelname-chains.csv")
            compiletimes_fn_nodate = abspath(results_path, "$modelname-compile_times.csv")

            for (file, link) in zip(
                (samplingtimes_fn, diagnostics_fn, chains_fn, compiletimes_fn),
                (samplingtimes_fn_nodate, diagnostics_fn_nodate, chains_fn_nodate, compiletimes_fn_nodate)
            )
                rm(link, force=true)
                symlink(file, link)
            end
        end

        @info "[$modelname] Finished successfully"
    else
        @error "Unknown model: $modelname"
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS...)
end


text:
### A Pluto.jl notebook ###
# v0.12.15

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 2fcbf600-19e9-11eb-3f4e-812a4a557c20
using Revise, PlutoUI, Aggregator, Aggregator.Api, Plots, StatsBase

# ╔═╡ eab20564-19fc-11eb-039f-dde6b6ef710a
begin
	using Aggregator.MangroveBase: ReviewsSummary
	using Aggregator.Model: MeanBeta
	using Turing, MCMCChains, StatsPlots
end

# ╔═╡ 9a68f5a6-33b8-11eb-3ca1-6f6ae024c66a
using Aggregator.Model: neutralities

# ╔═╡ 5e865152-19f3-11eb-0434-7b4fe9bde74c
md"""
# Looking around
## Fetch some data via the open API

See [documentation](https://docs.mangrove.reviews), to try it out directly.
"""

# ╔═╡ e6055cee-19ea-11eb-0a18-dd5847a9b78b
all_data = Api.all_data();

# ╔═╡ 6e1573f8-1e00-11eb-2c32-37d99db3a81d
all_data

# ╔═╡ 52f702b4-19f3-11eb-12db-bbf15e67ac47
md"## Define useful filters"

# ╔═╡ 4f4d918e-19ec-11eb-3b85-53cdeccf6c6c
begin
	reviews = all_data[:reviews]
	issuers = all_data[:issuers]
	maresi_subjects = all_data[:maresi_subjects]
	# Extract payload for each review.
	payloads = [r[:payload] for r in reviews]
end;

# ╔═╡ 4aff9e48-19ff-11eb-2e55-5f742ab32834
md"## View individual data points to understand the data"

# ╔═╡ 09950bc8-19f0-11eb-25cd-cbdad6720d00
md"We have **$(length(reviews))** reviews available, to browse through their payloads you can use the slider:"

# ╔═╡ 3a5fa0c4-19eb-11eb-0d57-8b3ae66cc2cb
@bind i Slider(1:length(payloads); show_value=true)

# ╔═╡ fd4eb326-19ea-11eb-064e-9765cb805381
payloads[i]

# ╔═╡ 0cc6d9a2-19ff-11eb-3151-edf0fd984a33
md"Images from selected review:"

# ╔═╡ ee090a7c-19f8-11eb-3b24-4bad4154fb47
begin
	image_urls = [i[:src] for i in get(payloads[i], :images, [])]
	Markdown.parse(join(["![]($image_url)" for image_url in image_urls]))
end

# ╔═╡ 623c0214-19f4-11eb-19f8-25bd8b4759c3
md"## Count properties to get global view"

# ╔═╡ b07ed06c-19f1-11eb-1e0d-f9cc921936f9
bar(countmap([r[:scheme] for r in reviews]); label="Number of reviews with given scheme")

# ╔═╡ afaa3110-19ef-11eb-1300-4f0997354a02
histogram([p[:rating]/25+1 for p in payloads if haskey(p, :rating)]; label="Number of reviews with given rating", legend=:left)

# ╔═╡ d88e9f14-19f2-11eb-25cb-515474ef1d7e
pie(countmap([haskey(p, :opinion) ? "Includes text" : "Only rating" for p in payloads]))

# ╔═╡ e37de39a-19fb-11eb-3dc7-519496db6281
md"# Aggregation algorithm"

# ╔═╡ aa351f6c-1a01-11eb-112a-91098dfc0d9a
md"""
Takes the reviews and gives us information about:
- quality of each subject
- neutrality of each reviewer
"""

# ╔═╡ 7bd88b90-19fc-11eb-284c-3bbaac7bfbd4
data = convert(ReviewsSummary, Api.extract(all_data));

# ╔═╡ 4d525f4c-1a01-11eb-1bf5-93d9f75f60c1
md"## Rating hypothesis/prior"

# ╔═╡ 13784e4c-1a00-11eb-2cde-75e886262ca1
@bind typical_rating Slider(0.1:0.1:0.9; default=0.5, show_value=true)

# ╔═╡ ceb59220-19fc-11eb-0f3c-13401192e6e5
rating_prior = MeanBeta(typical_rating);

# ╔═╡ bde4519e-19fd-11eb-24e2-af7f11145287
rating_sample = rand(rating_prior)

# ╔═╡ 0560381c-19fe-11eb-2687-03de0bd3f1f7
stars(raw_rating) = round(Int, raw_rating * 4) + 1;

# ╔═╡ a310a412-19fd-11eb-15ea-9f628cb434ea
stars(rating_sample)

# ╔═╡ 6d56a51c-1a00-11eb-11e4-a7c3f78de140
md"### Prior distribution of ratings"

# ╔═╡ 67f9ac00-1a84-11eb-1c54-7f71af6356e1
scatter(rand(rating_prior, 100); label="Sequence of draws from rating prior", legend=:right)

# ╔═╡ 7a9aa156-19fd-11eb-089b-1f1e9dc68932
histogram(map(stars, rand(rating_prior, 1000)); label="Number of reviews with given rating", legend=:left)

# ╔═╡ 80cba09a-1a01-11eb-0eca-474971cbd4d1
md"## Data driven prior"

# ╔═╡ 7bfb2c32-1a02-11eb-0a8c-7f0b9ff5675f
collect(data)[i]

# ╔═╡ 15852576-19fd-11eb-227c-97d39d2f8d0e
begin
	typical_rating_ = round(mean(values(data)), digits=2)
	md"Average rating value: **$(typical_rating_)**"
end

# ╔═╡ 4be0bac4-1a07-11eb-31c1-3d175b6fe365
md"## Monte Carlo inference"

# ╔═╡ 05475720-1a03-11eb-3c18-69b65f2ad768
chains = sample(mangrove_model(data, typical_rating_), PG(100), 100);

# ╔═╡ 932f8bfc-1a03-11eb-031b-3944a2157475
plot(chains; colordim=:parameter)

# ╔═╡ 90c60698-1a06-11eb-19ac-77b8fcd23c46
begin
	mean_qualities = mean(group(chains, :qualities))
	subjects = Dict(zip(subs(mean_qualities), qualities(mean_qualities)))
	mean_biases = mean(group(chains, :biases))
	reviewers = Dict(zip(Aggregator.Model.kids(mean_biases), Aggregator.Model.neutralities(mean_biases)))
end;

# ╔═╡ a12cb3c4-1a10-11eb-1494-199c763f04c1
subjects |> values |> collect |> histogram

# ╔═╡ 6ce62da8-1a11-11eb-22f6-573137e77293
worst_reviewer = collect(keys(reviewers))[argmin(collect(values(reviewers)))]

# ╔═╡ 8b3ad6c8-1a11-11eb-1f6a-077185c89857
[r[:payload] for r in reviews if r[:kid] == worst_reviewer]

# ╔═╡ 0609d5fa-19f1-11eb-13d9-25c772bc2a10
histogram(collect(values(reviewers)); label="Neutralities histogram", legend=:left)

# ╔═╡ c3f32ac6-1e06-11eb-0263-05927fbc8215
reviewers

# ╔═╡ d2ec9cf6-1e06-11eb-105a-edc6d5fbf2cc
subjects

# ╔═╡ d50c4ca0-1e08-11eb-29ba-cd06e4559d16
rand(Beta(0.04, 0.08), 10000) |> histogram

# ╔═╡ 8bf807d4-1e0a-11eb-2050-27bf41d7aa3e
rand(Beta(20, 1), 10000) |> histogram

# ╔═╡ 0adba460-1e0a-11eb-213d-09d963634bcd
rand(MeanBeta(0.8, 10), 10000) |> histogram

# ╔═╡ 81f061f4-33b9-11eb-179d-e78d8cb7e41e
plot(x -> pdf(MeanBeta(0.5, 2.), x), 0, 1)

# ╔═╡ 76768ae4-33ba-11eb-065f-2561fb2ea0fb
[rand(MeanBeta(0.73, 3))*rand(MeanBeta(0.5 / 101, 3)) for _ in 1:10000] |> mean

# ╔═╡ 98160576-33ba-11eb-3434-937ff3a831b6
issuers[Symbol("-----BEGIN PUBLIC KEY-----MFkwEwYHKoZIzj0CAQYIKoZIzj0DAQcDQgAEQs7EX2DHYSpif48XQBp3vjH4Df39kbqef7P5RkBRmDSQUffTW80BSuPM8rHjzSgehCy2VJHUPfG719PDKpUmMA==-----END PUBLIC KEY-----")]

# ╔═╡ cbbe7ec6-33d3-11eb-1e56-ed7b81097851
subjects["geo:47.557018,8.5930537?q=Burgstelle Hebelstein&u=30"]

# ╔═╡ e12fb6d8-33d3-11eb-1b74-6f0a7f437699
[p[:payload][:rating] for p in reviews if p[:kid] == "-----BEGIN PUBLIC KEY-----MFkwEwYHKoZIzj0CAQYIKoZIzj0DAQcDQgAEQs7EX2DHYSpif48XQBp3vjH4Df39kbqef7P5RkBRmDSQUffTW80BSuPM8rHjzSgehCy2VJHUPfG719PDKpUmMA==-----END PUBLIC KEY-----"]

# ╔═╡ 23a6f398-33d4-11eb-2659-7f2673bb9ab4
MeanBeta(0+0.5 / 101)

# ╔═╡ 8c958bc2-33d8-11eb-304c-71995153657b
reviews[1]

# ╔═╡ dcaee4ea-33d9-11eb-2bdb-056764479e43
Uniform() |> rand

# ╔═╡ Cell order:
# ╠═2fcbf600-19e9-11eb-3f4e-812a4a557c20
# ╟─5e865152-19f3-11eb-0434-7b4fe9bde74c
# ╠═e6055cee-19ea-11eb-0a18-dd5847a9b78b
# ╠═6e1573f8-1e00-11eb-2c32-37d99db3a81d
# ╟─52f702b4-19f3-11eb-12db-bbf15e67ac47
# ╠═4f4d918e-19ec-11eb-3b85-53cdeccf6c6c
# ╟─4aff9e48-19ff-11eb-2e55-5f742ab32834
# ╟─09950bc8-19f0-11eb-25cd-cbdad6720d00
# ╟─3a5fa0c4-19eb-11eb-0d57-8b3ae66cc2cb
# ╟─fd4eb326-19ea-11eb-064e-9765cb805381
# ╟─0cc6d9a2-19ff-11eb-3151-edf0fd984a33
# ╠═ee090a7c-19f8-11eb-3b24-4bad4154fb47
# ╟─623c0214-19f4-11eb-19f8-25bd8b4759c3
# ╠═b07ed06c-19f1-11eb-1e0d-f9cc921936f9
# ╠═afaa3110-19ef-11eb-1300-4f0997354a02
# ╠═d88e9f14-19f2-11eb-25cb-515474ef1d7e
# ╟─e37de39a-19fb-11eb-3dc7-519496db6281
# ╟─aa351f6c-1a01-11eb-112a-91098dfc0d9a
# ╠═eab20564-19fc-11eb-039f-dde6b6ef710a
# ╠═7bd88b90-19fc-11eb-284c-3bbaac7bfbd4
# ╟─4d525f4c-1a01-11eb-1bf5-93d9f75f60c1
# ╟─13784e4c-1a00-11eb-2cde-75e886262ca1
# ╠═ceb59220-19fc-11eb-0f3c-13401192e6e5
# ╠═bde4519e-19fd-11eb-24e2-af7f11145287
# ╠═0560381c-19fe-11eb-2687-03de0bd3f1f7
# ╠═a310a412-19fd-11eb-15ea-9f628cb434ea
# ╟─6d56a51c-1a00-11eb-11e4-a7c3f78de140
# ╠═67f9ac00-1a84-11eb-1c54-7f71af6356e1
# ╠═7a9aa156-19fd-11eb-089b-1f1e9dc68932
# ╟─80cba09a-1a01-11eb-0eca-474971cbd4d1
# ╠═7bfb2c32-1a02-11eb-0a8c-7f0b9ff5675f
# ╟─15852576-19fd-11eb-227c-97d39d2f8d0e
# ╟─4be0bac4-1a07-11eb-31c1-3d175b6fe365
# ╠═05475720-1a03-11eb-3c18-69b65f2ad768
# ╠═932f8bfc-1a03-11eb-031b-3944a2157475
# ╟─90c60698-1a06-11eb-19ac-77b8fcd23c46
# ╠═a12cb3c4-1a10-11eb-1494-199c763f04c1
# ╠═6ce62da8-1a11-11eb-22f6-573137e77293
# ╠═8b3ad6c8-1a11-11eb-1f6a-077185c89857
# ╠═9a68f5a6-33b8-11eb-3ca1-6f6ae024c66a
# ╠═0609d5fa-19f1-11eb-13d9-25c772bc2a10
# ╠═c3f32ac6-1e06-11eb-0263-05927fbc8215
# ╠═d2ec9cf6-1e06-11eb-105a-edc6d5fbf2cc
# ╠═d50c4ca0-1e08-11eb-29ba-cd06e4559d16
# ╠═8bf807d4-1e0a-11eb-2050-27bf41d7aa3e
# ╠═0adba460-1e0a-11eb-213d-09d963634bcd
# ╠═81f061f4-33b9-11eb-179d-e78d8cb7e41e
# ╠═76768ae4-33ba-11eb-065f-2561fb2ea0fb
# ╠═98160576-33ba-11eb-3434-937ff3a831b6
# ╠═cbbe7ec6-33d3-11eb-1e56-ed7b81097851
# ╠═e12fb6d8-33d3-11eb-1b74-6f0a7f437699
# ╠═23a6f398-33d4-11eb-2659-7f2673bb9ab4
# ╠═8c958bc2-33d8-11eb-304c-71995153657b
# ╠═dcaee4ea-33d9-11eb-2bdb-056764479e43

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove


text:
function opengene_homedir()
    dir = joinpath(homedir(), ".OpenGene")
    if !isdir(dir)
        mkpath(dir)
    end
    return dir
end

function parse_config()
    homedir = opengene_homedir()
    file = joinpath(homedir,"config.ini")
    config = Dict()
    if !isfile(file)
        return config
    end
    io = open(file)
    # TODO: work around for readall missing in master
    if isdefined(Base, :readstring)
        text = readstring(io)
    else
        text = readall(io)
    end
    lines = split(text, '\n')
    for line in lines
        line = strip(line)
        if length(line) < 3
            continue
        end
        items = split(line, '=')
        if length(items) < 2
            warn("error format of $line")
            warn("each line in the config file should have format key = value")
            continue
        end
        key = strip(items[1])
        value = strip(items[2])
        config[key] = value
    end
    return config
end

function opengene_datadir()
    config = parse_config()
    dir = joinpath(opengene_homedir(), "data")
    if haskey(config, "data_folder")
        dir = config["data_folder"]
    end
    if !isdir(dir)
        mkpath(dir)
    end
    return dir
end


text:
"""
    asnamedtuple(x) :: NamedTuple

# Examples
```jldoctest
julia> using DataUtils: asnamedtuple

julia> asnamedtuple(Dict(:a => 1, :b => 2))
(a = 1, b = 2)

julia> asnamedtuple(Dict("a" => 1, "b" => 2))
(b = 2, a = 1)

julia> asnamedtuple(1 + 2im)
(re = 1, im = 2)
```
"""
asnamedtuple(d::AbstractDict) =
    NamedTuple{Tuple(Symbol.(keys(d)))}(Tuple(values(d)))

asnamedtuple(obj) = NamedTuple{fieldnames(typeof(obj))}(fieldvalues(obj))

fieldvalues(obj) = ntuple(i -> getfield(obj, i), nfields(obj))


text:
"""
To implement a sub-type of this you need to implement the methods below.
"""
abstract type Results end
function get_base_power(r::T) where {T <: Results}
    error("get_base_power must be implemented for $T")
end

function get_variables(r::T) where {T <: Results}
    error("get_variables must be implemented for $T")
end

function get_parameters(r::T) where {T <: Results}
    error("get_parameters must be implemented for $T")
end

function get_total_cost(r::T) where {T <: Results}
    error("get_total_cost must be implemented for $T")
end

function get_optimizer_stats(r::T) where {T <: Results}
    error("get_optimizer_stats must be implemented for $T")
end

function get_timestamp(r::T) where {T <: Results}
    error("get_timestamp must be implemented for $T")
end

function write_results(r::T) where {T <: Results}
    error("write_results must be implemented for $T")
end


text:
const DCLR_EX_DAYS = 9
const DIST_DAYS = 12
const DIST_DAYS_EXTENDED = 40
const DIST_RET = :net
const DIST_RET_BOUND = 0.2

const DIST_MIN_ALLOWED = 1#DIST_DAYS*2
const DIST_DAYS_TO_STALE_PRICE = 22
const DIST_PRIMARY_YEAR_RANGE = YEAR_RANGE[]
const DIST_YEAR_RANGE = 1926:2018

function oldreaddist(;fname::String = CRSP_DISTRIBUTIONS_NAME,
  path::String=EVENT_PATH,
  incsvstream::Function = IN_CSV_STREAM,
  crspdateformat::DateFormat = CRSP_DATE_FORMAT)
  local dist::DataFrame

  dist = incsvstream("$path\\$fname.csv.lz4") |> CSV.read |> DataFrame
  names!(dist, (s::Symbol->Symbol(lowercase(string(s)))).(names(dist)))
  rename!(dist, [:dclrdt=>:dclrdate, :exdt=>:exdate])

  @warn "Using crsp event file from crsp, but without the checks against main crsp.
    (There are a bunch of checks that were used against the main crsp.
    Check the 11/20 backup for details.)"

  return dist
end

function preprocessdist(dist::DataFrame;
      distdays::Int = DIST_DAYS)::DataFrame
  for f ∈ [:dclrdate, :exdate]
    (Date <: eltype(dist[!,f])) && continue
    dist[!,f] = (i::MInt->ismissing(i) ? missing : Date("$i", CRSP_DATE_FORMAT)).(dist[!,f])
  end

  #drop some bad data
  dist = dist[completecases(dist[!, [:permno, :exdate, :distcd]]), :]

  #get the last dividend for table 10
  sort!(dist, [:permno, :exdate])
  dist.Lexdate = Vector{MDate}(missing, size(dist,1))
  sdist::SubDataFrame = view(dist, (cd::MInt-> (cd ≥ 1000) && (cd < 2000)).(dist.distcd),:)
  ssdists::GroupedDataFrame = groupby(sdist, :permno)
  @mpar for i ∈ 1:length(ssdists)
    ssdist::SubDataFrame = ssdists[i]
    for (i::Int, r::DataFrameRow) ∈ enumerate(eachrow(ssdist))
      (i==1) && continue
      r.Lexdate = ssdist[i-1,:exdate]
    end
  end

  #get the interval since the last dividend
  #dist.days2lastdiv = ((Ldt::MDate, dt::Date)->
  #  (ismissing(Ldt) ? missing : (dt-Ldt).value)).(dist.Lexdate, dist.exdate)

  #drop some more bad data
  dist = dist[completecases(dist[!, [:dclrdate, :divamt]]), :]

  #time between exdate and dclrdate is important
  dist.dclrexdays = (d::Day->d.value).(dist.exdate .- dist.dclrdate)
  dist.exyear = (year).(dist.exdate)

  #filter!(r::DataFrameRow->r.exyear ∈ yearrange, dist) #restrict to desired years

  #ensures integrity
  dist = dist[dist.dclrexdays .≥ 0, :]

  return dist
end

function primaryreaddist(;  prefix = error("prefix is required"),
  path::String=EVENT_PATH,
  incsvstream::Function = IN_CSV_STREAM,
  crspdateformat::DateFormat = CRSP_DATE_FORMAT,
  refreshdist::Bool=true,
  distyearrange::UnitRange = DIST_YEAR_RANGE
  )

  if REPLICATION_TYPE[] == :monthly
    @warn "REPLICATION TYPE set to monthly but monthly is not supported for event.
      Event processing will run using daily data."
  end

  local dist::DataFrame = prepcrsp(
    crspcolumns = [:date, :permno, :divamt, :distcd, :dclrdt],
    crsptype = :dailyexpanded,
    refreshcrsp=refreshdist,
    validpermnos = nothing,
    prefix = "altdist",
    yearrange=distyearrange) do df::DataFrame

    filter!(r::DataFrameRow->(!ismissing(r.distcd)), df)
    #df.Tprice = Vector{MFloat64}(missing, size(df,1))
    #(!isnothing(daystostaleprice)) && lagprice!(df)
  end

  rename!(dist, [:date=>:exdate, :dclrdt=>:dclrdate])

  println(describe(dist))

  return dist
end

function processdist!(dist::DataFrame,
    dclrexdays::Int = DCLR_EX_DAYS,
    Fret::Symbol = DIST_RET,
    retbound::Float64=DIST_RET_BOUND,
    minallowed = DIST_MIN_ALLOWED,
    daystostaleprice::Int=DIST_DAYS_TO_STALE_PRICE,
    days::Int = DIST_DAYS,
    daysextended::Int = DIST_DAYS_EXTENDED,
    primaryyearrange::UnitRange = DIST_PRIMARY_YEAR_RANGE,
    distyieldtol::Float64=10^-5)

  local sdists::GroupedDataFrame

  dist.primaryyearrange = (y::Int->y ∈ primaryyearrange).(dist.exyear)

  #now compute the dividend yield from the declaration date
  dist.yield = Vector{MFloat64}(undef, size(dist,1))
  dist.rawyield = Vector{MFloat64}(undef, size(dist,1))
  dist.yield2 = Vector{MFloat64}(undef, size(dist,1))
  dist.winsorizedyield = falses(size(dist,1))
  dist.dclrextradingdays = Vector{MInt}(undef, size(dist,1))
  dist.lowerdays = Vector{MInt}(undef, size(dist,1))
  dist.upperdays = Vector{MInt}(undef, size(dist,1))

  @assert issorted(dist, [:eid, :date])

  #lowerdays::Int = max(daystostaleprice,days,daysextended)
  #upperdays::Int = max(days, daysextended)
  sdists = groupby(dist, :eid)
  @mpar for i ∈ 1:length(sdists)
    sdist::SubDataFrame = sdists[i]
    #some basic integrity checks

    #find the upper and lower limits of the window
    local lowerdays::Int = -minimum(sdist.day)
    local upperdays::Int = maximum(sdist.day)
    sdist.lowerdays .= lowerdays
    sdist.upperdays .= upperdays

    @assert maximum(sdist.divamt) == minimum(sdist.divamt)
    @assert maximum(sdist.dclrexdays) == minimum(sdist.dclrexdays)
    @assert sdist.dclrdate[1] == sdist.date[lowerdays + 1] - Day(sdist.dclrexdays[1])

    #find the row for the declaration date
    #per the paper, we take the later of price at declaration or 22 days prior
    dclridx::Int = searchsortedfirst(sdist.date, sdist.dclrdate[1])
    (dclridx>size(sdist,1)) && continue
    ((sdist.day[dclridx]) < -daystostaleprice) && (dclridx = lowerdays+1-daystostaleprice)

    #compute the yield and winsorize if necessary
    yield::MFloat64 = sdist.divamt[1]/sdist.price[dclridx]
    sdist.rawyield .= ismissing(yield) || (!isfinite(yield)) ? missing : yield
    yield2::MFloat64 = sdist.divamt2[lowerdays+1]/sdist.price[dclridx]

    #no negative yields
    (!ismissing(yield)) && (yield < 0.0) && (yield=missing)
    (!ismissing(yield2)) && (yield2 < 0.0) && (yield2=missing)

    #if we can't get the yield from the event file, try the proper crsp file (yield2)
    ((ismissing(yield)) && (!ismissing(yield2)) && (!ismissing(dsdist.distcd2[lowerdays+1])) &&
      (dsdist.distcd[lowerdays+1] == sdist.distcd2[lowerdays+1]) && (yield=yield2))
    (ismissing(yield)) && continue

    sdist.dclrextradingdays .=  abs(sdist.day[dclridx])
    sdist.winsorizedyield .= yield > retbound
    sdist.yield .= min(yield,retbound)
    (!ismissing(yield2)) && (sdist.yield2 .= min(yield2,retbound))

  end
  #now determine the dividend yields from the exiration

  filter!(r::DataFrameRow->((r.day ≥ -daysextended) && (r.day ≤ daysextended)), dist)

  dist.primary = Vector{Bool}(undef, size(dist,1))
  dist.inrange = Vector{Bool}(undef, size(dist,1))
  dist.match = Vector{Bool}(undef, size(dist,1))
  dist.ordinary = Vector{Bool}(undef, size(dist,1))
  dist.dclrexok = Vector{Bool}(undef, size(dist,1))

  sdists = groupby(dist, :eid)
  @mpar for i ∈ 1:length(sdists)
    sdist = sdists[i]
    r::DataFrameRow = sdist[1, :] #examine fields common to each eid

    sdist.ordinary .= r.distcd == 1232
    sdist.dclrexok .= r.dclrexdays ≥ dclrexdays
    sdist.inrange .= (sdist.day .≥ -days) .* (sdist.day .≤ days)
    notmissing::Int = sum((!ismissing).(sdist[sdist.inrange,Fret]))
    sdist.match .= (notmissing ≥ minallowed) && (!ismissing(r.yield))

  end

  dist.primary .= (dist.ordinary .*
    dist.match .*
    dist.oneex .*
    dist.dclrexok .*
    dist.inrange .*
    dist.primaryyearrange)

  #println("inrange: $(sum(dist.inrange))")

  #some more integrity checks
  yieldmismatches::Int = 0
  avgyieldmismatch::Float64 = 0.0
  distcdmismatches::Int = 0
  dist.mismatches = falses(size(dist,1))
  for sdist ∈ groupby(dist, :eid)
    r::DataFrameRow = sdist[sdist.lowerdays[1]+1,:]
    @assert r.day == 0 #true by the prior filter

    if ((!ismissing(r.yield2)) && (!ismissing(r.yield)) && (r.distcd == r.distcd2) &&
      (abs(r.yield2 - r.yield)>distyieldtol))

      avgyieldmismatch += abs(r.yield2 - r.yield)
      yieldmismatches += 1
      sdist.mismatches .= true
    end

    if   ((!ismissing(r.distcd2)) && (!ismissing(r.distcd)) &&
      (r.distcd2 ≠ r.distcd))
      distcdmismatches += 1
    #  sdist.mismatches .= true
    end
  end
  avgyieldmismatch /= yieldmismatches

  #volume data for figure 4
  dist.dolvol = dist.vol .* dist.price
  dist.lnumtrd = (log).(dist.numtrd)
  dist.lvol = (log).(dist.vol)
  dist.ldolvol = (log).(dist.dolvol)
  (s::Symbol->finiteormissing!(dist, s)).([:lvol, :lnumtrd, :ldolvol])

  println("yieldmismatches: $yieldmismatches, avgyieldmismatch=$avgyieldmismatch
    distcdmismatches: = $distcdmismatches")

  #dist[(max(end-20_000,1)):end,:] |> CSV.write(
  #  "output\\distest-$(REPLICATION_TYPE[])-$(DIST_TYPE[]).csv")

  if occursin("short", string(DIST_TYPE[]))
    dist |> CSV.write(
      "output\\disttest-$(REPLICATION_TYPE[])-$(DIST_TYPE[]).csv")
  end
  dist[dist.mismatches,:] |> CSV.write(
    "output\\distmismatches-$(REPLICATION_TYPE[])-$(DIST_TYPE[]).csv")
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Add `Primitive` type to BalanceLaws, AtmosModel

# Vars wrapper
function prognostic_to_primitive!(
    bl::BalanceLaw,
    prim::AbstractArray,
    prog::AbstractArray,
    aux::AbstractArray,
)
    FT = eltype(prim)
    prognostic_to_primitive!(
        bl,
        Vars{vars_state(bl, Primitive(), FT)}(prim),
        Vars{vars_state(bl, Prognostic(), FT)}(prog),
        Vars{vars_state(bl, Auxiliary(), FT)}(aux),
    )
end
function primitive_to_prognostic!(
    bl::BalanceLaw,
    prog::AbstractArray,
    prim::AbstractArray,
    aux::AbstractArray,
)
    FT = eltype(prog)
    primitive_to_prognostic!(
        bl,
        Vars{vars_state(bl, Prognostic(), FT)}(prog),
        Vars{vars_state(bl, Primitive(), FT)}(prim),
        Vars{vars_state(bl, Auxiliary(), FT)}(aux),
    )
end

# By default the primitive is the prognostic
vars_state(bl::BalanceLaw, ::Primitive, FT) = vars_state(bl, Prognostic(), FT)

function prognostic_to_primitive!(bl, prim::Vars, prog::Vars, aux)
    prim_arr = parent(prim)
    prog_arr = parent(prog)
    prim_arr .= prog_arr
end
function primitive_to_prognostic!(bl, prog::Vars, prim::Vars, aux)
    prim_arr = parent(prim)
    prog_arr = parent(prog)
    prog_arr .= prim_arr
end

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA


text:
@testset "Simplification" begin
  c = Chain(P2[(0,0),(1,0),(1.5,0.5),(1,1),(0,1),(0,0)])
  s1 = simplify(c, DouglasPeucker(T(0.1)))
  s2 = simplify(c, DouglasPeucker(T(0.5)))
  @test s1 == Chain(P2[(0,0),(1,0),(1.5,0.5),(1,1),(0,1),(0,0)])
  @test s2 == Chain(P2[(0,0),(1.5,0.5),(1,1),(0,1),(0,0)])
end


text:
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using StatsFuns: logistic

function time_normal_response(
    n :: Int64,  # number of cases
    p :: Int64,  # number of predictors (SNPs)
    k :: Int64,  # number of true predictors per group
    debias :: Bool # whether to debias
)
	Random.seed!(1111)

	#construct snpmatrix, covariate files, and true model b
	x, maf = simulate_random_snparray(n, p)
	xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
	z = ones(n, 1) # non-genetic covariates, just the intercept
	true_b = zeros(p)
	true_b[1:k] = randn(k)
	shuffle!(true_b)
	correct_position = findall(x -> x != 0, true_b)
	noise = rand(Normal(0, 0.1), n) # noise vectors from N(0, s)

	#simulate phenotypes (e.g. vector y) via: y = Xb + noise
	y = xbm * true_b + noise

    #free memory from xbm.... does this work like I think it works?
	xbm = nothing
	GC.gc()

	#time the result and return
    result = @timed L0_normal_reg(x, z, y, 1, k, debias=debias)
    iter = result[1].iter
    runtime = result[2]      # seconds
    memory = result[3] / 1e6 # MB

    return iter, runtime, memory
end

function main(n :: Int64, i :: Int64)
	p = 1000000
	k = 10
	debias = false
	response = "normal"

	#run result: first run on super small data to force compile, then run result
	iter, runtime, memory = time_normal_response(100, 100, k, debias)
	iter, runtime, memory = time_normal_response(n, p, k, debias)

	#write result to corresponding file
	if debias
		dir = response * "_results/"
	else
		dir = response * "_results_nodebias/"
	end
	file = string(n) * "_by_" * string(p) * "_run$i"
	open(dir * file, "w") do f
		write(f, "time (seconds), memory (MB), iteration\n")
		write(f, string(runtime) * ", " * string(memory) * ", " * string(iter))
	end
end

# arguments passed from command line
n = parse(Int64, ARGS[1])
i = parse(Int64, ARGS[2])

main(n, i)


text:
using Test
using BasisFunctions, DomainIntegrals

using BasisFunctions: support

@testset "Normalization" begin
    F = Fourier(10)
    f = exp
    for g in (PeriodicEquispacedGrid(10,support(F)), interpolation_grid(F), MidpointEquispacedGrid(10,support(F)))
        s = sampling_normalization(GridBasis(g), FourierWeight())
        @test abs(sum(s*s*f)-applymeasure(FourierWeight(),f)) < .1
        s = sampling_normalization(GridBasis(g), lebesguemeasure(support(F)))
        @test abs(sum(s*s*f)-applymeasure(lebesguemeasure(support(F)),f)) < .1
    end

    C = ChebyshevT(9)
    f = exp
    g = ChebyshevNodes(9)
    s = sampling_normalization(GridBasis(g), lebesguemeasure(support(C)))
    @test abs(sum(s*s*f)-applymeasure(lebesguemeasure(support(C)),f)) < .1
    s = sampling_normalization(GridBasis(g), LegendreWeight{Float64}())
    @test abs(sum(s*s*f)-applymeasure(lebesguemeasure(support(C)),f)) < .1

    g = ChebyshevExtremae(9)
    s = sampling_normalization(GridBasis(g), lebesguemeasure(support(C)))
    @test abs(sum(s*s*f)-applymeasure(lebesguemeasure(support(C)),f)) < .1
    s = sampling_normalization(GridBasis(g), LegendreWeight{Float64}())
    @test abs(sum(s*s*f)-applymeasure(lebesguemeasure(support(C)),f)) < .1

    C = ChebyshevT(1000)
    f = exp
    g = MidpointEquispacedGrid(1000,-1,1)
    s = sampling_normalization(GridBasis(g), ChebyshevTWeight{Float64}())
    @test abs(sum(s*s*f)-applymeasure(ChebyshevTWeight{Float64}(),f)) < .1

    x, w = rectangular_rule(10)
    μ = gauss_rule(Fourier(10))
    @test points(μ) ≈ x
    @test weights(μ) ≈ w
end


text:
push!(LOAD_PATH,"../src/")

using Documenter
using DocumenterTools
using QuantumLattices

makedocs(
    format=     Documenter.HTML(
                    prettyurls= !("local" in ARGS),
                    canonical=  "https://quantum-many-body.github.io/QuantumLattices.jl/latest/",
                    assets=     ["assets/favicon.ico"],
                    analytics=  "UA-89508993-1",
                    ),
    clean=      false,
    sitename=   "QuantumLattices.jl",
    pages=      [
                "Home"      =>  "index.md",
                "Tutorials" =>  [
                    "Unitcell Description" => [
                        "tutorials/UnitcellDescription/Introduction.md",
                        "tutorials/UnitcellDescription/SpatialInfoOfAUnitcell.md",
                        "tutorials/UnitcellDescription/InternalDegreesOfFreedom.md",
                        "tutorials/UnitcellDescription/CouplingsAmongDifferentDegreesOfFreedom.md",
                        "tutorials/UnitcellDescription/GeneratorOfOperators.md",
                    ],
                    "Engine App Interface" => [
                        "tutorials/EngineAppInterface/Introduction.md",
                    ]
                ],
                "Manual"    =>  [
                    "man/Interfaces.md",
                    "Prerequisites" => [
                        "man/Prerequisites/Introduction.md",
                        "man/Prerequisites/TypeTraits.md",
                        "man/Prerequisites/Factories.md",
                        "man/Prerequisites/CompositeStructures.md",
                        "man/Prerequisites/SimpleTrees.md",
                        "man/Prerequisites/NamedVectors.md",
                        ],
                    "Mathematics" => [
                        "man/Mathematics/Introduction.md",
                        "man/Mathematics/Combinatorics.md",
                        "man/Mathematics/VectorSpaces.md",
                        "man/Mathematics/AlgebraOverFields.md",
                        "man/Mathematics/QuantumNumbers.md",
                        ],
                    "Essentials" => [
                        "man/Essentials/Introduction.md",
                        "man/Essentials/Spatials.md",
                        "man/Essentials/DegreesOfFreedom.md",
                        "man/Essentials/Terms.md",
                        "man/Essentials/FockPackage.md",
                        "man/Essentials/SpinPackage.md",
                        "man/Essentials/Frameworks.md",
                        ],
                    ],
                "Developer Documentation"   =>  [
                    "developer/CodeStyle.md",
                    ]
                ]
)

deploydocs(
    repo=       "github.com/Quantum-Many-Body/QuantumLattices.jl.git",
    target=     "build",
    deps=       nothing,
    make=       nothing,
)


text:
# Schrodinger.jl Tests
using Schrodinger
using Test, LinearAlgebra, SparseArrays

@testset "All Tests" begin
include("quobj.jl")
include("math.jl")
include("special.jl")
include("dynamics.jl")
include("superops.jl")
include("control.jl")
include("tomography.jl")
include("misc.jl")
end


text:
module TestBarnesHut
using LinearAlgebra
using CovarianceFunctions
using CovarianceFunctions: BarnesHutFactorization, barneshut!, vector_of_static_vectors,
        node_sums, euclidean, GradientKernel, taylor!
using NearestNeighbors
using NearestNeighbors: isleaf, getleft, getright, get_leaf_range
using Test

function barneshut_no_far_field!(b::AbstractVector, F::BarnesHutFactorization, w::AbstractVector,
                    α::Number = 1, β::Number = 0, θ::Real = F.θ)
    D = length(eltype(F.x))
    sums_a = node_sums(w, F.Tree)
    # centers_of_mass = compute_centers_of_mass(F.y, w, F.Tree)
    centers_of_mass = [hyp.center for hyp in F.Tree.hyper_spheres]
    for i in eachindex(F.x)
        Fa_i = bh_recursion_no_far_field(1, F.k, F.x[i], F.y, w, sums_a, θ, F.Tree, centers_of_mass)
        b[i] = α * Fa_i + β * b[i]
    end
    return b
end

function bh_recursion_no_far_field(index, k, x, y::AbstractVector,
            w::AbstractVector, sums_w::AbstractVector{<:Number}, θ::Real, T::BallTree, centers_of_mass)
    h = T.hyper_spheres[index]
    if isleaf(T.tree_data.n_internal_nodes, index) # do direct computation
        y_ind = get_leaf_range(T.tree_data, index)
        y_i, w_i = @views y[T.indices[y_ind]], w[T.indices[y_ind]]

        elty = promote_type(eltype(x), eltype(eltype(y_i)), eltype(w_i))
        val = zero(elty)
        @inbounds @simd for j in eachindex(y_i)
          val += k(x, y_i[j]) * w_i[j]
        end
        return val

    elseif h.r < θ * euclidean(x, centers_of_mass[index]) # compress

        return 0

    else # recurse

        l = bh_recursion_no_far_field(getleft(index), k, x, y, w, sums_w, θ, T, centers_of_mass) # IDEA: these two are spawnable
        r = bh_recursion_no_far_field(getright(index), k, x, y, w, sums_w, θ, T, centers_of_mass)
        return l + r
    end
end

verbose = false
# IDEA: loop over element types, loop over kernels
@testset "barneshut" begin
    n = 1024
    d = 2
    rand_sign = rand(n) .* rand((-1, 1), n)
    weight_vectors = [ones(n), rand(n), rand_sign, randn(n)]
    weight_names = ["ones", "rand", "rand + sign", "randn"]

    X = randn(d, n)
    x = vector_of_static_vectors(X)
    # k = (x, y) -> inv(1 + sum(abs2, x-y))
    k = CovarianceFunctions.Cauchy()
    # k = (x, y) -> CovarianceFunctions.Exp()(x, y) # making it anonymous to test independence to kernel types
    K = gramian(k, x)
    # k = CovarianceFunctions.EQ()
    # k = function (x, y)
    #         r = sum(abs2, x-y)
    #         iszero(r) ? 0 : inv(r^2)
    #     end

    # test factorization syntax
    # leafsize = 16 # could also loop over this
    D = 1e-2
    F = BarnesHutFactorization(k, x, x, D; θ = 1/8, leafsize = 16)

    @testset "weight vector $name" for (w, name) in zip(weight_vectors, weight_names)
        # display(sum(eigvals(K + D*I) .> 1e-6))

        b = K*w + D*w
        b_bh = F*w
        @test isapprox(b, b_bh, rtol = 1e-3)
        w_bh = \(F, b_bh, maxiter = 128, verbose = false)
        @test isapprox(b_bh, F*w_bh, rtol = 1e-3)
        # @test isapprox(w, w_bh, rtol = 1e-3) # usually 1e-2 and less robust to smaller diagonal
        if verbose
            println("relative errors")
            println(norm(b-b_bh) / norm(b))
            println(norm(b_bh - F*w_bh) / norm(b_bh))
            println(norm(w - w_bh) / norm(w))
        end
        θ = 0
        b_bh = zero(b)
        barneshut!(b_bh, F, w, 1, 0, θ)
        @test b ≈ b_bh # testing that barnes hut is exact with θ = 0

        nexp = 16
        err = zeros(nexp)
        err_no_split = zeros(nexp)
        err_hyp = zeros(nexp)
        err_hyp_no_split = zeros(nexp)
        err_nff = zeros(nexp)
        err_taylor = zeros(nexp)
        theta_array = range(1e-1, 1, length = nexp)
        for (i, θ) in enumerate(theta_array)
            barneshut!(b_bh, F, w, 1, 0, θ, split = true)
            err[i] = norm(b - b_bh)

            barneshut!(b_bh, F, w, 1, 0, θ, split = false)
            err_no_split[i] = norm(b - b_bh)

            barneshut_no_far_field!(b_bh, F, w, 1, 0, θ) # compare against pseudo barnes hut where far field = 0
            err_nff[i] = norm(b - b_bh)

            taylor!(b_bh, F, w, 1, 0, θ) # compare against pseudo barnes hut where far field = 0
            err_taylor[i] = norm(b - b_bh)
        end

        rel_err = err / norm(b)
        @test all(<(1e-1), rel_err[1:findlast(<(0.5), theta_array)])
        @test all(<(1e-2), rel_err[1:findlast(<(0.2), theta_array)])

        # plotting code
        # using Plots
        # plotly()
        #
        # rel_err_nff = err_nff / norm(b)
        # rel_err_no_split = err_no_split / norm(b)
        # rel_err_taylor = err_taylor / norm(b)
        #
        # plot(theta_array, rel_err, yscale = :log10, label = "barneshut", ylabel = "relative error", xlabel = "θ")
        # plot!(theta_array, rel_err_no_split, yscale = :log10, label = "no split")
        # plot!(theta_array, rel_err_nff, yscale = :log10, label = "sparse")
        # plot!(theta_array, rel_err_taylor, yscale = :log10, label = "taylor")
        # gui()

    end # testset weight vectors

    # @testset "gradient kernels" begin
    #     k = CovarianceFunctions.Cauchy()
    #     g = CovarianceFunctions.GradientKernel(k)
    #
    #     F = BarnesHutFactorization(g, x)
    #     @test F isa BarnesHutFactorization
    #     @test eltype(F) <: Diagonal
    #     @test size(F) == (n, n)
    #     @test size(F[1, 1]) == (d_out, d_out)
    # end # testset matrix valued bh

end # testset

end # module


text:
using MyPackage
using Test

@testset "MyPackage.jl" begin
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
end

max_stars_repo_name: Damien-Black/MyPackage

max_issues_repo_name: Damien-Black/MyPackage

max_forks_repo_name: Damien-Black/MyPackage


text:
# Tested and developed in Windows with Julia 1.5.1, October 2020.
# 2-D Topoplogy Optimization in Julia
# IMPORT PACKAGES
using Printf
using SparseArrays
using LinearAlgebra
using Plots: heatmap
# FUNCTIONS
function optimize!(F::Array{Float64,1}, U::Array{Float64,1},
                  freedofs::Array{Int64,1},
                  penal::Float64, E0::Float64, Emin::Float64, ft::Int64,
                  nelx::Int64, nely::Int64,
                  iK::Array{Int64,1},
                  jK::Array{Int64,1},
                  edofMat::Array{Int64,2},
                  KE::Array{Float64,2},
                  H::SparseMatrixCSC{Float64,Int64},
                  Hs::Array{Float64,2}, volfrac::Float64)::Array{Float64,2}
    x = repeat([volfrac],nely,nelx)
    xPhys = copy(x)
    xnew = copy(x)
    loop = 0
    change = 1.0
    sK = zeros(64*nelx*nely)
    K = spzeros(2*(nely+1)*(nelx+1),2*(nely+1)*(nelx+1))
    dc = zeros(nely, nelx)
    dv = zeros(nely, nelx)
    ce = zeros(nely, nelx)
    while change > 0.01
      loop = loop +1
      # FE-ANALYSIS
      sK = reshape(vec(KE)*(Emin.+vec(xPhys)'.^penal*(E0-Emin)),64*nelx*nely)
      K = sparse(iK,jK,sK)
      U[freedofs] .= Symmetric(K[freedofs,freedofs])\F
      # OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
      ce = reshape(sum((U[edofMat]*KE).*U[edofMat], dims=2),nely,nelx)
      c = sum((Emin.+xPhys.^penal.*(E0-Emin)).*ce)
      dc = -penal*(E0-Emin).*xPhys.^(penal-1).*ce
      dv = ones(nely,nelx)
      # FILTERING/MODIFICATION OF SENSITIVITIES
      if ft == 1 # sensitivity filtering
        dc = reshape(H*(vec(x).*vec(dc))./Hs./max.(1e-3,vec(x)),nely,nelx)
      elseif ft == 2 # density filtering
        dc = reshape(H*(vec(dc)./Hs),nely,nelx)
        dv = reshape(H*(vec(dv)./Hs),nely,nelx)
      end
      # OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES AND PHYSICAL DENSITIES
      xPhys, xnew = optimalityCriteria!(x, xPhys, xnew, dc, dv, H, Hs, nelx, nely, ft, volfrac)
      change = maximum(abs.(vec(xnew) - vec(x)))
      x = copy(xnew)
      s = @sprintf "It: %i  Obj.: %.2f  Ch: %.3f" loop c change;
      println(s)
    end
    return xPhys
end
function optimalityCriteria!(x::Array{Float64,2}, xPhys::Array{Float64,2}, xnew::Array{Float64,2},
                            dc::Array{Float64,2}, dv::Array{Float64,2}, H::SparseMatrixCSC{Float64,Int64},
                            Hs::Array{Float64,2}, nelx::Int64, nely::Int64, ft::Int64,
                            volfrac::Float64)::Tuple{Array{Float64,2},Array{Float64,2}}
    l1 = 0.0; l2 = 1.0e9; move = 0.2
    while (l2-l1)/(l1+l2) > 1.0e-3
      lmid = 0.5*(l2+l1)
      xnew = max.(0.0,max.(x.-move,min.(1.0,min.(x.+move,x.*sqrt.(-dc./dv./lmid)))))
      if ft == 1
        xPhys = copy(xnew)
      elseif ft == 2
        xPhys = reshape((H*vec(xnew))./Hs,nely,nelx)
      end
      if (sum(xPhys) > volfrac*nelx*nely) l1 = lmid else l2 = lmid end
    end
    return xPhys, xnew
end
function prepFilter(nelx::Int64, nely::Int64, rmin::Int64)::Tuple{SparseMatrixCSC{Float64,Int64}, Array{Float64,2}}
  iH = ones(nelx*nely*(2*(ceil(rmin)-1)+1)^2)
  jH = ones(size(iH))
  sH = zeros(size(iH))
  k = 0
  for i1 in 1:nelx
    for j1 in 1:nely
      e1 = (i1-1)*nely+j1
      for i2 in max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
        for j2 in max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
          e2 = (i2-1)*nely+j2
          k = k+1
          iH[k] = e1
          jH[k] = e2
          sH[k] = max(0.0,rmin-sqrt((i1-i2)^2+(j1-j2)^2))
        end
      end
    end
  end
  iH = Int.(iH)
  jH = Int.(jH)
  H = sparse(iH,jH,sH)
  Hs = sum(H, dims=2)
  return (H,Hs)
end
################################################################################
# MAIN CODE
function TopOptOC(nelx = 150::Int64, nely = 50::Int64,
                  rmin = 6::Int64, volfrac = 0.5::Float64,
                  penal = 3.0::Float64, ft = 1::Int64)
  # TOPOPT SETTINGS
  # nelx, elements in x
  # nely, elements in y
  # rmin, filter radius
  # volfrac, volume fraction
  # penal, SIMP penalty
  # ft, filtering option: 1 = sensitivity filter, 2 density filter
  nele = nelx * nely
  # MATERIAL PROPERTIES
  E0 = 1.0
  Emin = 1.0e-9
  nu = 0.3
  # ELEMENT STIFFNESS MATRIX
  A11 = Float64[12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12.0]
  A12 = Float64[-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6.0]
  B11 = Float64[-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4.0]
  B12 = Float64[2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2.0]
  KE = 1/(1-nu^2)/24 * ([A11 A12; A12' A11] + nu*[B11 B12; B12' B11])
  # CONNECTIVITY
  nodenrs = reshape(1:(1+nelx)*(1+nely), 1+nely, 1+nelx)
  edofVec = reshape(2*nodenrs[1:end-1,1:end-1]+ones(nely,nelx), nelx*nely, 1)
  edofMat = Int.(repeat(edofVec, 1, 8) + repeat([0 1 2*nely .+ [2 3 0 1] -2 -1], nelx*nely, 1))
  iK = reshape(kron(edofMat, ones(8,1))', 64*nelx*nely)
  jK = reshape(kron(edofMat, ones(1,8))', 64*nelx*nely)
  iK = convert(Array{Int64}, iK)
  jK = convert(Array{Int64}, jK)
  # DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
  F = zeros(2*(nely+1)*(nelx+1))
  F[2,1] = -1.0
  U = zeros(2*(nely+1)*(nelx+1))
  fixeddofs = union(collect(1:2:2*(nely+1)), [2*(nelx+1)*(nely+1)])
  alldofs = collect(1:2*(nely+1)*(nelx+1))
  freedofs = setdiff(alldofs,fixeddofs)
  # PREPARE FILTER
  H,Hs = prepFilter(nelx,nely,rmin)
  # INITIALIZE ITERATION
  x = Array{Float64}(undef, nely, nelx)
  x = optimize!(F[freedofs], U, freedofs, penal, E0, Emin, ft, nelx, nely, iK, jK, edofMat, KE, H, Hs, volfrac)
  # PLOT FINAL DESIGN
  heatmap(1.0.-x[end:-1:1,:],  yaxis=false, xaxis=false, legend = :none,
  color = :greys, grid=false, border=nothing, aspect_ratio=:equal)
end
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# This code was written by Nicolò Pollini,                                %
# Department of Wind Energy,                                              %
# Technical University of Denmark.                                        %
#                                                                         %
# Contact: nicolo@alumni.technion.ac.il                                   %
#                                                                         %
# Code repository: https://github.com/pollinico/topopt_jl                 %
#                                                                         %
# Disclaimer:                                                             %
# The author reserves all rights but does not guarantee that the code is  %
# free from errors. Furthermore, the author shall not be liable in any    %
# event caused by the use of the program.                                 %
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
include(joinpath(@__DIR__, "pedagogy_graph.jl"))
using Random

dims = [5, 5]
g = LG.grid(dims)

pp = PRS.SinkFreeGraph(g)
rng = Random.MersenneTwister(7)
sink_free_graph_pedagogy(pp, dims, rng)


text:
# # Calibration of `RiBasedVerticalDiffusivity` to five LESbrary simulations

# ## Install dependencies

# ```julia
# using Pkg
# pkg"add ParameterEstimocean, Oceananigans, CairoMakie"
# ```

using Oceananigans
using Oceananigans.Units
using ParameterEstimocean
using LinearAlgebra, CairoMakie, DataDeps, Distributions

using Oceananigans.TurbulenceClosures: RiBasedVerticalDiffusivity

# # Using LESbrary data
#
# `ParameterEstimocean.jl` provides paths to synthetic observations
# derived from high-fidelity large eddy simulations. In this example, we illustrate
# calibration of a turbulence parameterization to one of these simulations:

cases = ["free_convection",
         "strong_wind_weak_cooling",
         "weak_wind_strong_cooling",
         "strong_wind",
         "strong_wind_no_rotation"]

datapaths = [@datadep_str("two_day_suite_1m/$(case)_instantaneous_statistics.jld2") for case in cases]

times = [2hours, 12hours, 24hours]
field_names = (:b, :u, :v)
transformation = ZScore()
regrid = (1, 1, 32)

observations = [SyntheticObservations(path; field_names, times, transformation, regrid)
                for path in datapaths]

# Let's take a look at the observations. We define a few
# plotting utilities along the way to use later in the example:

colorcycle = [:black, :red, :darkblue, :orange, :pink1, :seagreen, :magenta2]
markercycle = [:rect, :utriangle, :star5, :circle, :cross, :+, :pentagon]

function make_figure_axes(n=1)
    fig = Figure(resolution=(1200, n*400))
    axs = []
    for i = 1:n
        ax_b = Axis(fig[i, 1], xlabel = "Buoyancy \n[cm s⁻²]", ylabel = "z [m]")
        ax_u = Axis(fig[i, 2], xlabel = "x-velocity, u \n[cm s⁻¹]")
        ax_v = Axis(fig[i, 3], xlabel = "y-velocity, v \n[cm s⁻¹]")
        push!(axs, (ax_b, ax_u, ax_v))
    end
    n == 1 && (axs = first(axs))
    return fig, axs
end

function plot_fields!(axs, b, u, v, label, color, grid=first(observations).grid)
    z = znodes(Center, grid)
    ## Note unit conversions below, eg m s⁻² -> cm s⁻²:
    lines!(axs[1], 1e2 * b, z; color, label)
    lines!(axs[2], 1e2 * u, z; color, label)
    lines!(axs[3], 1e2 * v, z; color, label)
    return nothing
end

# And then plot the evolution of the observed fields,

fig, axs = make_figure_axes()

for (i, obs) in enumerate(observations)
    Nt = length(obs.times)
    t = obs.times[end]
    fields = map(name -> interior(obs.field_time_serieses[name][Nt], 1, 1, :), field_names)
    plot_fields!(axs, fields..., "t = " * prettytime(t), colorcycle[i])
end

[axislegend(ax, position=:rb, merge=true, labelsize=10) for ax in axs]

save("multi_case_lesbrary_synthetic_observations.svg", fig); nothing # hide

# ![](multi_case_lesbrary_synthetic_observations.svg)

# # Calibration

ri_based_closure = RiBasedVerticalDiffusivity()

simulation = ensemble_column_model_simulation(observations;
                                              Nensemble = 60,
                                              architecture = CPU(),
                                              tracers = (:b, :e),
                                              closure = ri_based_closure)

# The simulation is initialized with neutral boundary conditions
# and a default time-step, which we modify for our particular problem:

Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
Qᵇ = simulation.model.tracers.b.boundary_conditions.top.condition
N² = simulation.model.tracers.b.boundary_conditions.bottom.condition

simulation.Δt = 20minutes

for (i, obs) in enumerate(observations)
    view(Qᵘ, :, i) .= obs.metadata.parameters.momentum_flux
    view(Qᵇ, :, i) .= obs.metadata.parameters.buoyancy_flux
    view(N², :, i) .= obs.metadata.parameters.N²_deep
end

# We identify a subset of the closure parameters to calibrate by specifying
# parameter names and prior distributions:

priors = (ν₀   = lognormal(mean=0.01, std=0.005),
          κ₀   = lognormal(mean=0.1,  std=0.05),
          Ri₀ν = Normal(-0.5, 1.0),
          Ri₀κ = Normal(-0.5, 1.0),
          Riᵟν = lognormal(mean=1.0,  std=0.5),
          Riᵟκ = lognormal(mean=1.0,  std=0.5))

free_parameters = FreeParameters(priors)

# TODO: explain the meaning of each parameter
# The prior information comes from experience, prior calibration runs,
# and educated guesses.

calibration = InverseProblem(observations, simulation, free_parameters)

# Next, we calibrate, using a relatively large noise to reflect our
# uncertainty about how close the observations and model can really get,

eki = EnsembleKalmanInversion(calibration; pseudo_stepping=ConstantConvergence(0.8))
iterate!(eki; iterations = 10)
@show eki.iteration_summaries[end]

# # Results
#
# To analyze the results, we build a new simulation with just one ensemble member
# to evaluate some utilities for analyzing the results:

Nt = length(first(observations).times)
Niter = length(eki.iteration_summaries) - 1
modeled_time_serieses = calibration.time_series_collector.field_time_serieses

observed, modeled = [], []
for (c, obs) in enumerate(observations)
    push!(observed, map(name -> interior(obs.field_time_serieses[name][Nt], 1, 1, :), field_names))
    push!(modeled,  map(name -> interior(  modeled_time_serieses[name][Nt], 1, c, :), field_names))
end

function compare_model_observations(model_label="modeled")
    fig, axs = make_figure_axes(length(observations))
    for (c, obs) in enumerate(observations)
        plot_fields!(axs[c], observed[c]..., "observed at t = " * prettytime(times[end]), :black)
        plot_fields!(axs[c], modeled[c]..., model_label, :blue)
        [axislegend(ax, position=:rb, merge=true, labelsize=10) for ax in axs[c]]
    end
    return fig
end

# Now we execute forward runs for the initial ensemble mean,

initial_parameters = eki.iteration_summaries[0].ensemble_mean
forward_run!(calibration, initial_parameters)
fig = compare_model_observations("modeled after 0 iterations")

save("multi_case_model_observation_comparison_iteration_0.svg", fig); nothing # hide
# ![](multi_case_model_observation_comparison_iteration_0.svg)

# and the final ensemble mean, representing our "best" parameter set,

best_parameters = eki.iteration_summaries[end].ensemble_mean
forward_run!(calibration, best_parameters)
fig = compare_model_observations("modeled after $Niter iterations")

save("multi_case_model_observation_comparison_final_iteration.svg", fig); nothing # hide
# ![](multi_case_model_observation_comparison_final_iteration.svg)

# ## Parameter evolution
#
# To understand how results changed over the EKI iterations,
# we look at the evolution of the ensemble means,

ensemble_means = NamedTuple(n => map(summary -> summary.ensemble_mean[n], eki.iteration_summaries)
                            for n in calibration.free_parameters.names)

fig = Figure()
ax = Axis(fig[1, 1], xlabel = "Ensemble Kalman iteration", ylabel = "Parameter value")

for (i, name) in enumerate(calibration.free_parameters.names)
    label = string(name)
    marker = markercycle[i]
    color = colorcycle[i]
    scatterlines!(ax, 0:Niter, parent(ensemble_means[name]); marker, color, label)
end

axislegend(ax, position=:rb)

save("lesbrary_catke_parameter_evolution.svg", fig); nothing # hide

# ![](lesbrary_catke_parameter_evolution.svg)


text:
T = 0.0

"""
    LanguageServerInstance(pipe_in, pipe_out, debug=false, env="", depot="")

Construct an instance of the language server.

Once the instance is `run`, it will read JSON-RPC from `pipe_out` and
write JSON-RPC from `pipe_in` according to the [language server
specification](https://microsoft.github.io/language-server-protocol/specifications/specification-3-14/).
For normal usage, the language server can be instantiated with
`LanguageServerInstance(stdin, stdout, false, "/path/to/environment")`.

# Arguments
- `pipe_in::IO`: Pipe to read JSON-RPC from.
- `pipe_out::IO`: Pipe to write JSON-RPC to.
- `debug::Bool`: Whether to log debugging information with `Base.CoreLogging`.
- `env::String`: Path to the
  [environment](https://docs.julialang.org/en/v1.2/manual/code-loading/#Environments-1)
  for which the language server is running. An empty string uses julia's
  default environment.
- `depot::String`: Sets the
  [`JULIA_DEPOT_PATH`](https://docs.julialang.org/en/v1.2/manual/environment-variables/#JULIA_DEPOT_PATH-1)
  where the language server looks for packages required in `env`.
"""
mutable struct LanguageServerInstance
    jr_endpoint::JSONRPCEndpoints.JSONRPCEndpoint
    workspaceFolders::Set{String}
    _documents::Dict{URI2,Document}

    debug_mode::Bool
    runlinter::Bool
    ignorelist::Set{String}
    isrunning::Bool

    env_path::String
    depot_path::String
    symbol_server::SymbolServer.SymbolServerInstance
    symbol_results_channel::Channel{Any}
    symbol_store::Dict{String,SymbolServer.ModuleStore}
    symbol_store_ready::Bool
    # ss_task::Union{Nothing,Future}
    format_options::DocumentFormat.FormatOptions
    lint_options::StaticLint.LintOptions

    combined_msg_queue::Channel{Any}

    err_handler::Union{Nothing,Function}

    status::Symbol

    number_of_outstanding_symserver_requests::Int

    current_symserver_progress_token::Union{Nothing,String}

    clientcapability_window_workdoneprogress::Bool

    function LanguageServerInstance(pipe_in, pipe_out, debug_mode::Bool = false, env_path = "", depot_path = "", err_handler=nothing)
        new(
            JSONRPCEndpoints.JSONRPCEndpoint(pipe_in, pipe_out, err_handler),
            Set{String}(),
            Dict{URI2,Document}(),
            debug_mode,
            true,
            Set{String}(),
            false,
            env_path,
            depot_path,
            SymbolServer.SymbolServerInstance(depot_path),
            Channel(Inf),
            deepcopy(SymbolServer.stdlibs),
            false,
            DocumentFormat.FormatOptions(),
            StaticLint.LintOptions(),
            Channel{Any}(Inf),
            err_handler,
            :created,
            0,
            nothing,
            false
        )
    end
end
function Base.display(server::LanguageServerInstance)
    println("Root: ", server.workspaceFolders)
    for d in getdocuments_value(server)
        display(d)
    end
end

function hasdocument(server::LanguageServerInstance, uri::URI2)
    return haskey(server._documents, uri)
end

function getdocument(server::LanguageServerInstance, uri::URI2)
    return server._documents[uri]
end

function getdocuments_key(server::LanguageServerInstance)
    return keys(server._documents)
end

function getdocuments_pair(server::LanguageServerInstance)
    return pairs(server._documents)
end

function getdocuments_value(server::LanguageServerInstance)
    return values(server._documents)
end

function setdocument!(server::LanguageServerInstance, uri::URI2, doc::Document)
    server._documents[uri] = doc
end

function deletedocument!(server::LanguageServerInstance, uri::URI2)
    delete!(server._documents, uri)
end

function create_symserver_progress_ui(server)
    if server.clientcapability_window_workdoneprogress
        server.current_symserver_progress_token = string(uuid4())
        response = JSONRPCEndpoints.send_request(server.jr_endpoint, "window/workDoneProgress/create", Dict("token" => server.current_symserver_progress_token))

        JSONRPCEndpoints.send_notification(server.jr_endpoint, "\$/progress", Dict("token" => server.current_symserver_progress_token, "value" => Dict("kind"=>"begin", "title"=>"Julia Language Server", "message"=>"Indexing packages...")))
    end
end

function destroy_symserver_progress_ui(server)
    if server.clientcapability_window_workdoneprogress
        progress_token = server.current_symserver_progress_token
        server.current_symserver_progress_token = nothing
        JSONRPCEndpoints.send_notification(server.jr_endpoint, "\$/progress", Dict("token" => progress_token, "value" => Dict("kind"=>"end")))
    end
end

function trigger_symbolstore_reload(server::LanguageServerInstance)
    server.symbol_store_ready = false
    if server.number_of_outstanding_symserver_requests==0 && server.status==:running
        create_symserver_progress_ui(server)
    end
    server.number_of_outstanding_symserver_requests += 1

    @async try
        # TODO Add try catch handler that links into crash reporting
        ssi_ret, payload = SymbolServer.getstore(server.symbol_server, server.env_path)

        server.number_of_outstanding_symserver_requests -= 1

        if server.number_of_outstanding_symserver_requests==0
            destroy_symserver_progress_ui(server)
        end

        if ssi_ret==:success
            push!(server.symbol_results_channel, payload)
        elseif ssi_ret==:failure
            error("The symbol server failed with '$(String(take!(payload)))'")
        end
        server.symbol_store_ready = true
    catch err
        bt = catch_backtrace()
        if server.err_handler!==nothing
            server.err_handler(err, bt)
        else
            Base.display_error(stderr, err, bt)
        end
    end
end

"""
    run(server::LanguageServerInstance)

Run the language `server`.
"""
function Base.run(server::LanguageServerInstance)
    server.status=:started

    run(server.jr_endpoint)

    trigger_symbolstore_reload(server)

    @async try
        while true
            msg = JSONRPCEndpoints.get_next_message(server.jr_endpoint)
            put!(server.combined_msg_queue, (type=:clientmsg, msg=msg))
        end
    catch err
        bt = catch_backtrace()
        if server.err_handler!==nothing
            server.err_handler(err, bt)
        else
            Base.display_error(stderr, err, bt)
        end
    end

    @async try
        while true
            msg = take!(server.symbol_results_channel)
            put!(server.combined_msg_queue, (type=:symservmsg, msg=msg))
        end
    catch err
        bt = catch_backtrace()
        if server.err_handler!==nothing
            server.err_handler(err, bt)
        else
            Base.display_error(stderr, err, bt)
        end
    end

    while true
        message = take!(server.combined_msg_queue)

        if message.type==:clientmsg
            msg = message.msg

            request = parse(JSONRPC.Request, msg)

            res = process(request, server)

            if request.id!=nothing
                JSONRPCEndpoints.send_success_response(server.jr_endpoint, msg, res)
            end
        elseif message.type==:symservmsg
            @info "Received new data from Julia Symbol Server."
            msg = message.msg

            server.symbol_store = msg
            roots = Document[]
            for doc in getdocuments_value(server)
                # only do a pass on documents once
                root = getroot(doc)
                if !(root in roots)
                    push!(roots, root)
                    scopepass(root, doc)
                end

                StaticLint.check_all(getcst(doc), server.lint_options, server)
                empty!(doc.diagnostics)
                mark_errors(doc, doc.diagnostics)
                publish_diagnostics(doc, server)
            end
        end
    end
end


text:
"""
Struct for WilsonFermion
"""
struct WilsonFermion_4D_wing{NC} <: AbstractFermionfields_4D{NC}
    NC::Int64
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    NG::Int64
    NDW::Int64
    f::Array{ComplexF64,6}

    γ::Array{ComplexF64,3}
    rplusγ::Array{ComplexF64,3}
    rminusγ::Array{ComplexF64,3}
    hop::Float64 #Hopping parameter
    r::Float64 #Wilson term
    hopp::Array{ComplexF64,1}
    hopm::Array{ComplexF64,1}
    eps::Float64
    Dirac_operator::String
    MaxCGstep::Int64
    BoundaryCondition::Array{Int8,1}

    function WilsonFermion_4D_wing(NC::T,NX::T,NY::T,NZ::T,NT::T,r,hop,eps,MaxCGstep,BoundaryCondition) where T<: Integer
        γ,rplusγ,rminusγ = mk_gamma(r)
        NG = 4
        NDW = 1
        hopp = zeros(ComplexF64,4)
        hopm = zeros(ComplexF64,4)
        hopp .= hop
        hopm .= hop
        #@assert NDW == 1 "only NDW = 1 is supported. Now NDW = $NDW"
        f = zeros(ComplexF64,NC,NX+2NDW,NY+2NDW,NZ+2NDW,NT+2NDW,NG)
        Dirac_operator = "Wilson"
        return new{NC}(NC,NX,NY,NZ,NT,NG,NDW,f,
                γ,rplusγ,rminusγ,hop,r,hopp,hopm,eps,Dirac_operator,MaxCGstep,BoundaryCondition)
    end

    function WilsonFermion_4D_wing(params,NC::T,NX::T,NY::T,NZ::T,NT::T) where T<: Integer
        #eps = params["eps_CG"]
        eps= check_parameters(params,"eps_CG",default_eps_CG)
        #println("eps_CG = ",eps_CG)
        MaxCGstep = check_parameters(params,"MaxCGstep",default_MaxCGstep)
        BoundaryCondition = check_parameters(params,"boundarycondition",[1,1,1,-1])
        r = check_parameters(params,"r",1.0)

        #MaxCGstep = params["MaxCGstep"]
        #BoundaryCondition = params["BoundaryCondition"]
        #BoundaryCondition = params["boundarycondition"]
        #r = params["r"]
        hop = params["κ"]
        return WilsonFermion_4D_wing(NC,NX,NY,NZ,NT,r,hop,eps,MaxCGstep,BoundaryCondition)
    end

end

function Base.similar(x::T) where T <: WilsonFermion_4D_wing
    return WilsonFermion_4D_wing(x.NC,x.NX,x.NY,x.NZ,x.NT,
                x.r,x.hop,x.eps,x.MaxCGstep,x.BoundaryCondition)
end

using InteractiveUtils

function Wx!(xout::T,U::Array{G,1},
    x::T,temps::Array{T,1}) where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields}
    temp = temps[4]
    temp1 = temps[1]
    temp2 = temps[2]

    clear_fermion!(temp)
    set_wing_fermion!(x)
    for ν=1:4

        xplus = shift_fermion(x,ν)

        println("enrer")
        mul!(temp1,U[ν],x)
        error("not")
        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication

        mul!(temp1,view(x.rminusγ,:,:,ν))

        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)

        #
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rplusγ,:,:,ν),temp2)
        mul!(temp2,view(x.rplusγ,:,:,ν))

        add_fermion!(temp,x.hopp[ν],temp1,x.hopm[ν],temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1,x,-1,temp)

    #display(xout)
    #    exit()
    return
end

function Wdagx!(xout::T,U::Array{G,1},
    x::T,temps::Array{T,1}) where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields}
    temp = temps[4]
    temp1 = temps[1]
    temp2 = temps[2]

    clear_fermion!(temp)
    set_wing_fermion!(x)
    for ν=1:4
        xplus = shift_fermion(x,ν)
        mul!(temp1,U[ν],x)
        error("wdag")

        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication
        #mul!(temp1,view(x.rminusγ,:,:,ν),temp1)
        mul!(temp1,view(x.rplusγ,:,:,ν))

        #
        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rminusγ,:,:,ν),temp2)
        mul!(temp2,view(x.rminusγ,:,:,ν))

        add_fermion!(temp,x.hopp[ν],temp1,x.hopm[ν],temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1,x,-1,temp)

    #display(xout)
    #    exit()
    return
end

function Wx!(xout::T,U::Array{G,1},
    x::T,temps::Array{T,1},CloverFμν::TC) where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields,TC <: AbstractArray}
    temp = temps[4]
    temp1 = temps[1]
    temp2 = temps[2]

    clear_fermion!(temp)
    set_wing_fermion!(x)
    for ν=1:4
        xplus = shift_fermion(x,ν)
        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication

        mul!(temp1,view(x.rminusγ,:,:,ν))

        xminus = shift_fermion(x,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)

        #
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rplusγ,:,:,ν),temp2)
        mul!(temp2,view(x.rplusγ,:,:,ν))

        add_fermion!(temp,x.hopp[ν],temp1,x.hopm[ν],temp2)

    end

    clear_fermion!(xout)
    add_fermion!(xout,1,x,-1,temp)

    #cloverterm!(xout,CloverFμν,x)

    #display(xout)
    #    exit()
    return
end

function Wdagx!(xout::T,U::Array{G,1},
    x::T,temps::Array{T,1},CloverFμν::TC) where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields,TC <: AbstractArray}
    temp = temps[4]
    temp1 = temps[1]
    temp2 = temps[2]

    clear_fermion!(temp)
    x5 = temps[3]

    mul_γ5x!(x5,x)
    set_wing_fermion!(x5)
    #set_wing_fermion!(x)

    for ν=1:4
        xplus = shift_fermion(x5,ν)
        mul!(temp1,U[ν],xplus)

        #fermion_shift!(temp1,U,ν,x)

        #... Dirac multiplication
        #mul!(temp1,view(x.rminusγ,:,:,ν),temp1)
        mul!(temp1,view(x.rminusγ,:,:,ν))

        #
        xminus = shift_fermion(x5,-ν)
        Uminus = shift_U(U[ν],-ν)

        mul!(temp2,Uminus',xminus)
        #fermion_shift!(temp2,U,-ν,x)
        #mul!(temp2,view(x.rminusγ,:,:,ν),temp2)
        mul!(temp2,view(x.rplusγ,:,:,ν))

        add_fermion!(temp,x.hopp[ν],temp1,x.hopm[ν],temp2)

    end

    clear_fermion!(temp1)
    add_fermion!(temp1,1,x5,-1,temp)
    #clear_fermion!(xout)
    #add_fermion!(xout,1,x,-1,temp)

    #cloverterm!(temp1,CloverFμν,x5)
    mul_γ5x!(xout,temp1)

    #display(xout)
    #    exit()
    return
end

function LinearAlgebra.mul!(x::WilsonFermion_4D_wing{NC},A::TA) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            x[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            x[ic,ix,iy,iz,it,2] = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            x[ic,ix,iy,iz,it,3] = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            x[ic,ix,iy,iz,it,4] = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_wing{NC},A::TA,x::WilsonFermion_4D_wing{NC}) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            xout[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            xout[ic,ix,iy,iz,it,2] = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            xout[ic,ix,iy,iz,it,3] = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            xout[ic,ix,iy,iz,it,4] = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_wing{NC},x::WilsonFermion_4D_wing{NC},A::TA) where {TA <: AbstractMatrix, NC}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            xout[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[2,1]*e2+A[3,1]*e3+A[4,1]*e4
                            xout[ic,ix,iy,iz,it,2] = A[1,2]*e1+A[2,2]*e2+A[3,2]*e3+A[4,2]*e4
                            xout[ic,ix,iy,iz,it,3] = A[1,3]*e1+A[2,3]*e2+A[3,3]*e3+A[4,3]*e4
                            xout[ic,ix,iy,iz,it,4] = A[1,4]*e1+A[2,4]*e2+A[3,4]*e3+A[4,4]*e4

                    end
                end
            end
        end
    end

end

function set_wing_fermion!(a::WilsonFermion_4D_wing{NC}) where NC
    NT = a.NT
    NZ = a.NZ
    NY = a.NY
    NX = a.NX

    #!  X-direction
    for ialpha=1:4
        for it=1:NT
            for iz = 1:NZ
                for iy=1:NY
                    @simd for k=1:NC
                        a[k,0,iy,iz,it,ialpha] = a.BoundaryCondition[1]*a[k,NX,iy,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for k=1:NC
                        a[k,NX+1,iy,iz,it,ialpha] = a.BoundaryCondition[1]*a[k,1,iy,iz,it,ialpha]
                    end
                end
            end
        end
    end

    #Y-direction
    for ialpha = 1:4
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,0,iz,it,ialpha] = a.BoundaryCondition[2]*a[k,ix,NY,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,NY+1,iz,it,ialpha] = a.BoundaryCondition[2]*a[k,ix,1,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        # Z-direction
        for it=1:NT
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,iy,0,it,ialpha] = a.BoundaryCondition[3]*a[k,ix,iy,NZ,it,ialpha]
                        a[k,ix,iy,NZ+1,it,ialpha] = a.BoundaryCondition[3]*a[k,ix,iy,1,it,ialpha]

                    end
                end
            end
        end

        #T-direction
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,iy,iz,0,ialpha] = a.BoundaryCondition[4]*a[k,ix,iy,iz,NT,ialpha]
                        a[k,ix,iy,iz,NT+1,ialpha] = a.BoundaryCondition[4]*a[k,ix,iy,iz,1,ialpha]
                    end
                end
            end
        end

    end

end

"""
c--------------------------------------------------------------------------c
c     y = gamma_5 * x
c     here
c                  ( -1       )
c        GAMMA5 =  (   -1     )
c                  (     +1   )
c                  (       +1 )
c--------------------------------------------------------------------------c
    """
    function mul_γ5x!(y::WilsonFermion_4D_wing{NC},x::WilsonFermion_4D_wing{NC}) where NC
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        for ig=1:4
            for ic=1:NC
                for it=1:NT
                    for iz=1:NZ
                        for iy=1:NY
                            for ix=1:NX
                                @simd for ic=1:NC
                                    y[ic,ix,iy,iz,it,ig] =x[ic,ix,iy,iz,it,ig]*ifelse(ig <= 2,-1,1)
                                end
                            end
                        end
                    end
                end
            end
        end
    end

function cloverterm!(vec,CloverFμν,x::WilsonFermion_4D_wing{NC}) where NC
    NT = x.NT
    NZ = x.NZ
    NY = x.NY
    NX = x.NX

    i  =0
    for it=1:NT
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    i += 1
                    for k1=1:NC
                        for k2=1:NC

                            c1 = x[k2,ix,iy,iz,it,1]
                            c2 = x[k2,ix,iy,iz,it,2]
                            c3 = x[k2,ix,iy,iz,it,3]
                            c4 = x[k2,ix,iy,iz,it,4]

                            vec[k1,ix,iy,iz,it,1] += CloverFμν[k1,k2,i,1]*(-   c1) +
                                                        + CloverFμν[k1,k2,i,2]*(-im*c2) +
                                                        + CloverFμν[k1,k2,i,3]*(-   c2) +
                                                        + CloverFμν[k1,k2,i,4]*(-   c2) +
                                                        + CloverFμν[k1,k2,i,5]*( im*c2) +
                                                        + CloverFμν[k1,k2,i,6]*(-   c1)

                            vec[k1,ix,iy,iz,it,2] += CloverFμν[k1,k2,i,1]*(   c2) +
                                                        + CloverFμν[k1,k2,i,2]*(im*c1) +
                                                        + CloverFμν[k1,k2,i,3]*(-   c1) +
                                                        + CloverFμν[k1,k2,i,4]*(-   c1) +
                                                        + CloverFμν[k1,k2,i,5]*(-im*c1) +
                                                        + CloverFμν[k1,k2,i,6]*(   c2)

                            vec[k1,ix,iy,iz,it,3] += CloverFμν[k1,k2,i,1]*(   -c3) +
                                                        + CloverFμν[k1,k2,i,2]*(-im*c4) +
                                                        + CloverFμν[k1,k2,i,3]*(   c4) +
                                                        + CloverFμν[k1,k2,i,4]*(-   c4) +
                                                        + CloverFμν[k1,k2,i,5]*(-im*c4) +
                                                        + CloverFμν[k1,k2,i,6]*(   c3)

                            vec[k1,ix,iy,iz,it,4] += CloverFμν[k1,k2,i,1]*(   c4) +
                                                        + CloverFμν[k1,k2,i,2]*(im*c3) +
                                                        + CloverFμν[k1,k2,i,3]*(   c3) +
                                                        + CloverFμν[k1,k2,i,4]*(-   c3) +
                                                        + CloverFμν[k1,k2,i,5]*(im*c3) +
                                                        + CloverFμν[k1,k2,i,6]*( -  c4)

                        end
                    end
                end
            end

        end
    end

    #println("vec = ",vec*vec)

end


text:
include("default_params.jl")
include("reflecting_model/model_def.jl")

ks_errors = DataFrame(DG_limit=[])
l1_errors = DataFrame(DG_limit=[])
l1_cell_probs_errors = DataFrame(DG_limit=[])
ks_error_row = zeros(1)
l1_error_row = zeros(1)
l1_cell_probs_row = zeros(1)

sim_cdf = CSV.read(
    (@__DIR__)*"/reflecting_model/data/"*ic*"/data/sims_cdf_evaluated.csv",
    DataFrame,
)

for o in orders
        approx_data = CSV.read(
            (@__DIR__)*"/reflecting_model/data/"*ic*"/data/cdfs/order_"*string(o)*"_model_limiter.csv",
            DataFrame,
        )
        ks_error = max(
            maximum(abs.(approx_data.phase_1-sim_cdf.phase_1)),
            maximum(abs.(approx_data.phase_2-sim_cdf.phase_2)),
        )
        ks_error_row[1] = ks_error
        l1_cdf = sum(abs.(approx_data.phase_1[1:end-1]+approx_data.phase_1[2:end]-
                sim_cdf.phase_1[1:end-1]-sim_cdf.phase_1[2:end]).*diff(approx_data.x))/2 +
            sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                sim_cdf.phase_2[1:end-1]-sim_cdf.phase_2[2:end]).*diff(approx_data.x))/2
        l1_error_row[1] = l1_cdf
        cell_edges_idx = findall(mod.(approx_data.x,1.0).==0.0)
        sim_cell_probs = DataFrame(
            x=[approx_data.x[1]; approx_data.x[cell_edges_idx[2:end]] - approx_data.x[cell_edges_idx[1:end-1]]],
            phase_1=[sim_cdf.phase_1[1]; sim_cdf.phase_1[cell_edges_idx[2:end]] - sim_cdf.phase_1[cell_edges_idx[1:end-1]]],
            phase_2=[sim_cdf.phase_2[1]; sim_cdf.phase_2[cell_edges_idx[2:end]] - sim_cdf.phase_2[cell_edges_idx[1:end-1]]],
        )
        approx_cell_probs = DataFrame(
            x=[approx_data.x[1]; approx_data.x[cell_edges_idx[2:end]] - approx_data.x[cell_edges_idx[1:end-1]]],
            phase_1=[approx_data.phase_1[1]; approx_data.phase_1[cell_edges_idx[2:end]] - approx_data.phase_1[cell_edges_idx[1:end-1]]],
            phase_2=[approx_data.phase_2[1]; approx_data.phase_2[cell_edges_idx[2:end]] - approx_data.phase_2[cell_edges_idx[1:end-1]]],
        )
        l1_cell_probs_row[1] = sum(abs.(approx_cell_probs.phase_1-sim_cell_probs.phase_1)) +
            sum(abs.(approx_cell_probs.phase_2-sim_cell_probs.phase_2))
    push!(ks_errors,ks_error_row)
    push!(l1_errors,l1_error_row)
    push!(l1_cell_probs_errors,l1_cell_probs_row)
end

CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/ks_error.csv",ks_errors)
CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/l1_error.csv",l1_errors)
CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error.csv",l1_cell_probs_errors)

ks_errors_lwr = DataFrame(DG_limit=[])
l1_errors_lwr = DataFrame(DG_limit=[])
l1_cell_probs_errors_lwr = DataFrame(DG_limit=[])
ks_error_row_lwr = zeros(1)
l1_error_row_lwr = zeros(1)
l1_cell_probs_error_row_lwr = zeros(1)

ks_errors_upr = DataFrame(DG_limit=[])
l1_errors_upr = DataFrame(DG_limit=[])
l1_cell_probs_errors_upr = DataFrame(DG_limit=[])
ks_error_row_upr = zeros(1)
l1_error_row_upr = zeros(1)
l1_cell_probs_error_row_upr = zeros(1)

ks_error_boot = zeros(n_boot)
l1_error_boot = zeros(n_boot)
l1_cell_probs_error_boot = zeros(n_boot)
for o in orders
    # for (c,m) in enumerate(string.(approx_types[[1;3;4]]).*string.([1;3;4]))
        approx_data = CSV.read(
            (@__DIR__)*"/reflecting_model/data/"*ic*"/data/cdfs/order_"*string(o)*"_model_limiter.csv",
            DataFrame,
        )
        cell_edges_idx = findall(mod.(approx_data.x,1.0).==0.0)
        approx_cell_probs = DataFrame(
                x=[approx_data.x[1]; approx_data.x[cell_edges_idx[2:end]] - approx_data.x[cell_edges_idx[1:end-1]]],
                phase_2=[approx_data.phase_2[1]; approx_data.phase_2[cell_edges_idx[2:end]] - approx_data.phase_2[cell_edges_idx[1:end-1]]],
                phase_1=[approx_data.phase_1[1]; approx_data.phase_1[cell_edges_idx[2:end]] - approx_data.phase_1[cell_edges_idx[1:end-1]]],
            )
        for n in 1:n_boot
            sim_cdf = CSV.read(
                (@__DIR__)*"/reflecting_model/data/"*ic*"/data/bootstrap/sample_"*string(n)*".csv",
                DataFrame,
            )
            ks_error = max(
                maximum(abs.(approx_data.phase_2-sim_cdf.phase_2)),
                maximum(abs.(approx_data.phase_1-sim_cdf.phase_1)),
            )
            ks_error_boot[n] = ks_error
            l1_cdf = sum(abs.(approx_data.phase_2[1:end-1]+approx_data.phase_2[2:end]-
                    sim_cdf.phase_2[1:end-1]-sim_cdf.phase_2[2:end]).*diff(approx_data.x))/2 +
                sum(abs.(approx_data.phase_1[1:end-1]+approx_data.phase_1[2:end]-
                    sim_cdf.phase_1[1:end-1]-sim_cdf.phase_1[2:end]).*diff(approx_data.x))/2
            l1_error_boot[n] = l1_cdf
            sim_cell_probs = DataFrame(
                x=[approx_data.x[1]; approx_data.x[cell_edges_idx[2:end]] - approx_data.x[cell_edges_idx[1:end-1]]],
                phase_2=[sim_cdf.phase_2[1]; sim_cdf.phase_2[cell_edges_idx[2:end]] - sim_cdf.phase_2[cell_edges_idx[1:end-1]]],
                phase_1=[sim_cdf.phase_1[1]; sim_cdf.phase_1[cell_edges_idx[2:end]] - sim_cdf.phase_1[cell_edges_idx[1:end-1]]],
            )
            l1_cell_probs_error_boot[n] = sum(abs.(approx_cell_probs.phase_2-sim_cell_probs.phase_2)) +
                sum(abs.(approx_cell_probs.phase_1-sim_cell_probs.phase_1))
        end
        ks_error_row_lwr[1] = quantile(ks_error_boot,qtiles[1])
        l1_error_row_lwr[1] = quantile(l1_error_boot,qtiles[1])
        l1_cell_probs_error_row_lwr[1] = quantile(l1_cell_probs_error_boot,qtiles[1])
        ks_error_row_upr[1] = quantile(ks_error_boot,qtiles[2])
        l1_error_row_upr[1] = quantile(l1_error_boot,qtiles[2])
        l1_cell_probs_error_row_upr[1] = quantile(l1_cell_probs_error_boot,qtiles[2])
    push!(ks_errors_lwr,ks_error_row_lwr)
    push!(l1_errors_lwr,l1_error_row_lwr)
    push!(l1_cell_probs_errors_lwr,l1_cell_probs_error_row_lwr)
    push!(ks_errors_upr,ks_error_row_upr)
    push!(l1_errors_upr,l1_error_row_upr)
    push!(l1_cell_probs_errors_upr,l1_cell_probs_error_row_upr)
end

CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/ks_error_lwr.csv",ks_errors_lwr)
CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/l1_error_lwr.csv",l1_errors_lwr)
CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error_lwr.csv",l1_cell_probs_errors_lwr)
CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/ks_error_upr.csv",ks_errors_upr)
CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/l1_error_upr.csv",l1_errors_upr)
CSV.write((@__DIR__)*"/reflecting_model/data/"*ic*"/data/errors/l1_cell_probs_error_upr.csv",l1_cell_probs_errors_upr)