text:
# Load and save results, statistics and run parameters
function save_results_csv(results::Dict, file_name; save_order=DEFAULT_SAVE_ORDER, dir="")
	AC_inputs = results["inputs"]
	AC_outputs = results["outputs"]
	duals = results["duals"]
	return save_results_csv(AC_inputs, AC_outputs, duals, file_name, save_order=save_order, dir=dir)
end

"""
Saves the AC\\_output, AC\\_input, duals to a csv file
"""
function save_results_csv(AC_inputs, AC_outputs, duals, file_name; save_order=DEFAULT_SAVE_ORDER, dir="")

	data = results_to_array(AC_inputs, AC_outputs, duals, save_order=save_order, header=false, imag_j=true)

	save_path = joinpath(dir, file_name*".csv")

	#TASK: Create file directory automatically to save files?
	if isfile(save_path)  # File already exists append data to it
		open(save_path, "a") do io
		writedlm(io, data, ',')
		end
	else  # Build header for network and create new file
		header = build_header(AC_inputs, AC_outputs, duals, save_order)
		data = vcat(header, data)
		writedlm(save_path, data, ',')
	end
end

"""
Saves the given sample to the given csv file, appending it to the end of the file.
"""
function save_sample_csv(x, results, file_name, save_order, dual_vars; dir="")
	# Gen. & Ext. and Dual Results
	primal_vals, dual_vals = results

    # Load results
	new_AC_inputs = extract_load(x)

    # Duals results
	dual_vals["unique_active_sets"] = extract_active_set(dual_vals, dual_vars)

	save_results_csv(new_AC_inputs, primal_vals, dual_vals, file_name, save_order=save_order, dir=dir)
end

"""
Builds the header as a new row to add to a data array with input, output, then dual data.
The order of the header is dictated by the variables in the save_order argument with each
element of a variable becoming a single column. Uses element labels to determine the type
of element the result variable refers to.
"""
function build_header(AC_inputs, AC_outputs, duals, save_order)
	header = []

	all_data = merge(AC_inputs, AC_outputs, duals)

	# Remove any keys in save_order that are not in the given results
	all_data_keys = keys(all_data)
	matched_label_indxs = [key in all_data_keys for key in save_order]
	save_order = save_order[matched_label_indxs]

	for key in save_order
		element = ELEMENT_LABELS[key]
		append!(header, ["$(element)$(i):"*string(key) for i in 1:size(all_data[key],2)])
	end

	return reshape(header,(1,length(header)))
end

"""
Saves the A and b matrices defining the polytope Ax < b to two seperate comma
delimited csv files.
"""
function save_polytope(A, b, file_name; dir="")
    save_path_A = joinpath(dir, file_name * "_A" * ".csv")
    save_path_b = joinpath(dir, file_name * "_b" * ".csv")

    writedlm(save_path_A, A, ',')
    writedlm(save_path_b, b, ',')
end

"""
Saves an infeasibility certificate by appending new rows for the A and b matrices
to csv files. If these csv files do not already exist then calls save_polytope
to save the entire A and b matrices in new csv files.
"""
function save_cert(A, b, file_name; dir="")
	A_last = A[end, :]'
	b_last = b[end, :]

    save_path_A = joinpath(dir, file_name * "_A" * ".csv")
    save_path_b = joinpath(dir, file_name * "_b" * ".csv")

	if isfile(save_path_A) & isfile(save_path_b)  # File already exists append data to it
		for (x, save_path) in zip((A_last, b_last), (save_path_A, save_path_b))
			open(save_path, "a") do io
				writedlm(io, x, ',')
			end
		end
	else  # Build header for network and create new file
		save_polytope(A, b, file_name, dir=dir)
	end
end

"""
Saves the last iteration of level 1 stats to a csv file.
If the csv file does not already exist then calls save_stats
to save the stats data with a header to a new csv file.
"""
function save_stats(stats, file_name; dir="")
    save_path = joinpath(dir, file_name * ".csv")

	stat_labels = add_dim_h(sort(collect(keys(stats))))

	if stats[stat_labels[1]] isa Array && size(stats[stat_labels[1]], 1) > 1  # More than one iteration of data
		data = stats[stat_labels[1]]
	else
		data = [stats[stat_labels[1]]]
	end
	saved_stat_labels = [stat_labels[1]]
	for stat_label in stat_labels[2:end]
		if !(stats[stat_label] isa Dict)
			push!(saved_stat_labels, stat_label)
			data = hcat(data, stats[stat_label])
		end
	end
	saved_stat_labels = add_dim_h(saved_stat_labels)

	if isfile(save_path) # File already exists append data to it
		open(save_path, "a") do io
		writedlm(io, data, ',')
		end
	else  # Create new file with header
		stat_data = vcat(saved_stat_labels, data)
		writedlm(save_path, stat_data, ',')
	end
end

"Makes an array, x, have a second dimension and be vertically oriented"
add_dim_v(x::Array) = reshape(x, (size(x)...,1))

"Makes an array, x, have a second dimension and be horizontally oriented"
add_dim_h(x::Array) = reshape(x, (1, size(x)...))


text:
@sde_model CoxIngersollRoss dr = κ*(θ-r)*dt + σ*sqrt(r)*dW

function _cir_transition_params(m, s, x0)
  @unpack κ, θ, σ = m
  c = 2 * κ / (σ^2 * (1 - exp(-κ * s.Δt)))
  df = 4 * κ * θ / σ^2
  nc = x0 * 2c * exp(-κ * s.Δt)
  c, df, nc
end

function sample(model::CoxIngersollRoss, scheme::Exact, t0, x0)
  c, df, nc = _cir_transition_params(model, scheme, x0)
  rand(NoncentralChisq(df, nc)) / 2c
end

function transition(model::CoxIngersollRoss, scheme::Exact, t0, x0, x1)
  c, df, nc = _cir_transition_params(model, scheme, x0)
  2c * pdf(NoncentralChisq(df, nc), x1 * 2c)
end

function logtransition(model::CoxIngersollRoss, scheme::Exact, t0, x0, x1)
  c, df, nc = _cir_transition_params(model, scheme, x0)
  logpdf(NoncentralChisq(df, nc), x1 * 2c) + log(2c)
end

ext: jl


text:
"""
    estimate_rule_counts(grammar, terminals; n_samples = length(terminals)^2)

Parse `terminals` using `grammar` and calculate a Monte Carlo estimate of
of the expected rule usage in the parse tree for all dependent category
components and rules, return a tuple of the expected counts as a matrix and
the likelihood as a `LogProb`

# Assumptions
- the number of used rules is constant for each parse tree
- `grammar` has a score named `:prob` to calculate the likelihood
- `grammar` has a score named `:forest` from which trees can be sampled as
  iterable objects of tuples of category indices and rule indices
"""
function estimate_rule_counts(
        grammar::Grammar, terminals; n_samples=length(terminals)^2)

    m, n = length(grammar.depcomps), length(grammar.all_rules)

    p = try
        parse(grammar, terminals)
    catch e
        println(e)
        println("not parsable: ", terminals)
        return fill(0.0, m, n), LogProb(1)
    end

    likelihood           = sum(values(p)).prob
    parse_forest         = sum(values(p)).forest
    random_tree_iterator = rand(parse_forest)
    n_rules              = length(collect(random_tree_iterator))

    trees = fill((0,0), n_rules, n_samples)
    for j in 1:n_samples
        trees[:, j] = collect(random_tree_iterator)
    end

    counts = foldr(vec(trees), init = fill(0, m, n)) do (c, r), M
        M[grammar.catidx2dcompidx[c], r] += 1
        M
    end

    counts / n_samples, likelihood
end

"""
    variational_probs(pseudocounts::Union{AbstractVector, AbstractMatrix})

Calculate the updated variational probabilities according to the `pseudocounts`.

Rowwise calculation if pseudocounts is a matric
"""
function variational_probs(pseudocounts::AbstractVector)
    s = sum(pseudocounts)
    normalize(map(a ->exp(digamma(a)), pseudocounts), 1)
end

function variational_probs(pseudocounts::AbstractMatrix)
    mapslices(variational_probs, pseudocounts, dims=2)
end

"""
    train_grammar(grammar, dataset)

Train `grammar` iterating through the `dataset` once.

The result is an updated grammar and the likelihood of the dataset according to
the initial grammar.

# Assumptions
- the number of used rules is constant for each parse tree
- `grammar` has a score named `:prior` that returns real numbers
- `grammar` has a score named `:prob` to calculate the likelihood
- `grammar` has a score named `:forest` from which trees can be sampled as
  iterable objects of tuples of category indices and rule indices
"""
function train_grammar(grammar::Grammar, dataset)
    # add counts, multiply likelihoods
    reduce((c1, p1), (c2, p2)) = (c1 + c2, p1 * p2)

    rule_counts, likelihood = @distributed (reduce) for terminals in dataset
        estimate_rule_counts(grammar, terminals)
    end

    pseudocounts = rule_counts + tabulate_score(grammar, :prior)

    vprobs = variational_probs(pseudocounts)
    probs  = mapslices(row -> normalize(row, 1), pseudocounts, dims=2)

    vprob_score = @closure (catidx, ruleidx) ->
        LogProb(vprobs[grammar.catidx2dcompidx[catidx], ruleidx])

    prob_score  = @closure (catidx, ruleidx) ->
        LogProb(probs[grammar.catidx2dcompidx[catidx], ruleidx])

    trained_grammar = set_scores(grammar,
        (prior = grammar.scores.prior,
         vprob = vprob_score,
         prob  = prob_score)
    )

    add_forest_score(trained_grammar, :vprob), likelihood
end


text:
module POMDPGallery

export gen_readme, run_scripts

function gen_readme(output=joinpath(dirname(@__FILE__()), "..", "README.md"))
    readme = IOBuffer()
    print(readme, """
        # POMDPGallery

        [![Build Status](https://travis-ci.org/JuliaPOMDP/POMDPGallery.jl.svg?branch=master)](https://travis-ci.org/JuliaPOMDP/POMDPGallery.jl)

        A gallery of models written for [POMDPs.jl](https://github.com/JuliaPOMDP/POMDPs.jl) with visualizations. You should be able to copy and paste the code below each visualization to run it on your local machine.

        For instructions on how to add new models, see [INSTRUCTIONS.md](INSTRUCTIONS.md).

        For the older version of this package with julia-0.6 models, see [this branch](https://github.com/JuliaPOMDP/POMDPGallery.jl/tree/julia-0.6).
        """)

    problemsdir = joinpath(dirname(@__FILE__()), "..", "problems")
    for problem in readdir(problemsdir)
        problemdir = joinpath(problemsdir, problem)
        url = strip(read(joinpath(problemdir, "url.txt"), String))
        desc = read(joinpath(problemdir, "description.txt"), String)

        println(readme, "## [$problem]($url)\n")
        println(readme, desc*"\n")

        println(readme, "![$problem](problems/$problem/out.gif)\n")
        println(readme, """
            ```julia
            $(strip(read(joinpath(problemdir, "script.jl"), String)))
            ```

            """)
    end
    file = open(output, "w")
    println(file, String(take!(readme)))
    close(file)
    return true
end

function run_scripts(;allow_failure=String[])
    pkgdir = joinpath(dirname(@__FILE__()), "..")
    problemsdir = joinpath(pkgdir, "problems")
    problems = readdir(problemsdir)
    results = similar(problems, Any)

    # @sync for (i, problem) in enumerate(problems)
    for (i, problem) in enumerate(problems)
        # @async begin
        begin
            println("Launching $problem...")
            problemdir = joinpath(problemsdir, problem)
            script = joinpath(problemdir, "script.jl")
            runs = """cd("$problemdir"); using Pkg; Pkg.activate("."); Pkg.instantiate(); include("$script")"""
            # outfile = joinpath(problemdir, "stdout.log")
            # errfile = joinpath(problemdir, "stderr.log")
            # pipe = pipeline(`julia --project=$pkgdir -e $runs`, stdout=outfile, stderr=errfile)
            # procs[i] = run(pipe, wait=false)
            # wait(procs[i])
            results[i] = try
                run(`julia --project=$pkgdir -e $runs`)
            catch ex
                ex
            end

            println("Finished $problem: $(passed(results[i]) ? "passed" : "FAILED")")
            # if !success(procs[i])
            #     println("""
            #             ====================
            #             stderr for $problem:
            #             ====================
            #             """)
            #     println(read(outfile, String))
            #     println("""
            #             ====================
            #             stdout for $problem:
            #             ====================
            #             """)
            #     println(read(errfile, String))
            # end
        end
    end

    pass = true
    println("""\n\n\n
            =====================
            POMDPGallery Summary:
            =====================

            """)
    for i in 1:length(problems)
        print(problems[i]*": ")
        if passed(results[i])
            println("passed")
        else
            println("FAILED")
            if !(problems[i] in allow_failure)
                pass = false
            end
        end
    end
    println("\n\n")

    return pass
end

passed(proc) = success(proc)
passed(ex::ProcessFailedException) = false

end # module


text:
using FFTW

function linear(U, temps, T, Ntau)

    dtau   = T / Ntau
    Ubis   = vcat(U, transpose(U[1, :]))
    temps  = temps % T
    repere = temps / dtau + 1
    indice = trunc(Int64,repere)

    @assert (1 <= indice <= Ntau)

    (indice + 1 - repere) * Ubis[indice, :] .+ (repere - indice) * Ubis[indice + 1, :]

end

function trigo(U, temps, T, Ntau)

    W   = zeros(ComplexF64,Ntau)
    W  .= vcat(0:Ntau÷2-1,-Ntau÷2:-1)
    W   = exp.(1im * 2 * pi / T * W * temps)
    V   = fft(U, 1)

    sum(V .* W, dims=1) / Ntau

end

function reconstr(U, temps, T, Ntau, type = 2)

    if type == 1

        res = linear(U, temps, T, Ntau)

    else

        res = trigo(U, temps, T, Ntau)

    end

    res

end


text:
# Taken from https://github.com/piotrmurach/tty-spinner/blob/master/lib/tty/spinner/formats.rb

const FORMATS = Dict(
      :classic => Dict(
        :frequency => 10,
        :frames => ["|", "/", "-", "\\"]
      ),
      :spin => Dict(
        :frequency => 10,
        :frames => ["◴", "◷", "◶", "◵"]
      ),
      :spin_2 => Dict(
        :frequency => 10,
        :frames => ["◐" "◓" "◑" "◒"]
      ),
      :spin_3 => Dict(
        :frequency => 10,
        :frames => ["◰", "◳", "◲", "◱"]
      ),
      :spin_4 => Dict(
        :frequency => 10,
        :frames => ["╫, ╪"]
      ),
      :pulse => Dict(
        :frequency => 10,
        :frames => ["⎺", "⎻", "⎼", "⎽", "⎼", "⎻"]
      ),
      :pulse_2 => Dict(
        :frequency => 15,
        :frames => ["▁", "▃", "▅", "▆", "▇", "█", "▇", "▆", "▅", "▃"]
      ),
      :pulse_3 => Dict(
        :frequency => 20,
        :frames => ["▉", "▊", "▋", "▌", "▍", "▎", "▏", "▎", "▍", "▌", "▋", "▊", "▉"]
      ),
      :dots => Dict(
        :frequency => 10,
        :frames => ["⠋", "⠙", "⠹", "⠸", "⠼", "⠴", "⠦", "⠧", "⠇", "⠏"]
      ),
      :dots_2 => Dict(
        :frequency => 10,
        :frames => ["⣾", "⣽", "⣻", "⢿", "⡿", "⣟", "⣯", "⣷"]
      ),
      :dots_3 => Dict(
        :frequency => 10,
        :frames => ["⠋", "⠙", "⠚", "⠞", "⠖", "⠦", "⠴", "⠲", "⠳", "⠓"]
      ),
      :dots_4 => Dict(
        :frequency => 10,
        :frames => ["⠄", "⠆", "⠇", "⠋", "⠙", "⠸", "⠰", "⠠", "⠰", "⠸", "⠙", "⠋", "⠇", "⠆"]
      ),
      :dots_5 => Dict(
        :frequency => 10,
        :frames => ["⠋", "⠙", "⠚", "⠒", "⠂", "⠂", "⠒", "⠲", "⠴", "⠦", "⠖", "⠒", "⠐", "⠐", "⠒", "⠓", "⠋"]
      ),
      :dots_6 => Dict(
        :frequency => 10,
        :frames => ["⠁", "⠉", "⠙", "⠚", "⠒", "⠂", "⠂", "⠒", "⠲", "⠴", "⠤", "⠄", "⠄", "⠤", "⠴", "⠲", "⠒", "⠂", "⠂", "⠒", "⠚", "⠙", "⠉", "⠁"]
      ),
      :dots_7 => Dict(
        :frequency => 10,
        :frames => ["⠈", "⠉", "⠋", "⠓", "⠒", "⠐", "⠐", "⠒", "⠖", "⠦", "⠤", "⠠", "⠠", "⠤", "⠦", "⠖", "⠒", "⠐", "⠐", "⠒", "⠓", "⠋", "⠉", "⠈"]
      ),
      :dots_8 => Dict(
        :frequency => 10,
        :frames => ["⠁", "⠁", "⠉", "⠙", "⠚", "⠒", "⠂", "⠂", "⠒", "⠲", "⠴", "⠤", "⠄", "⠄", "⠤", "⠠", "⠠", "⠤", "⠦", "⠖", "⠒", "⠐", "⠐", "⠒", "⠓", "⠋", "⠉", "⠈", "⠈"]
      ),
      :dots_9 => Dict(
        :frequency => 10,
        :frames => ["⢹", "⢺", "⢼", "⣸", "⣇", "⡧", "⡗", "⡏"]
      ),
      :dots_10 => Dict(
        :frequency => 10,
        :frames => ["⢄", "⢂", "⢁", "⡁", "⡈", "⡐", "⡠"]
      ),
      :dots_11 => Dict(
        :frequency => 10,
        :frames => ["⠁", "⠂", "⠄", "⡀", "⢀", "⠠", "⠐", "⠈"]
      ),
      :arrow => Dict(
        :frequency => 10,
        :frames => ["←", "↖", "↑", "↗", "→", "↘", "↓", "↙"]
      ),
      :arrow_pulse => Dict(
        :frequency => 10,
        :frames => [
          "▹▹▹▹▹",
          "▸▹▹▹▹",
          "▹▸▹▹▹",
          "▹▹▸▹▹",
          "▹▹▹▸▹",
          "▹▹▹▹▸"
        ]
      ),
      :triangle => Dict(
        :frequency => 10,
        :frames => ["◢", "◣", "◤", "◥"]
      ),
      :arc => Dict(
        :frequency => 10,
        :frames => ["◜", "◠", "◝", "◞", "◡", "◟"]
      ),
      :pipe => Dict(
        :frequency => 10,
        :frames => ["┤", "┘", "┴", "└", "├", "┌", "┬", "┐"]
      ),
      :bouncing => Dict(
        :frequency => 10,
        :frames => [
          "[    ]",
          "[   =]",
          "[  ==]",
          "[ ===]",
          "[====]",
          "[=== ]",
          "[==  ]",
          "[=   ]"
        ]
      ),
      :bouncing_ball => Dict(
        :frequency => 10,
        :frames => [
          "( ●    )",
          "(  ●   )",
          "(   ●  )",
          "(    ● )",
          "(     ●)",
          "(    ● )",
          "(   ●  )",
          "(  ●   )",
          "( ●    )",
          "(●     )"
        ]
      ),
      :bounce => Dict(
        :frequency => 10,
        :frames => ["⠁", "⠂", "⠄", "⠂"]
      ),
      :box_bounce => Dict(
        :frequency => 10,
        :frames => ["▌", "▀", "▐", "▄", ""]
      ),
      :box_bounce_2 => Dict(
        :frequency => 10,
        :frames => ["▖", "▘", "▝", "▗"]
      ),
      :star => Dict(
        :frequency => 10,
        :frames => ["✶", "✸", "✹", "✺", "✹", "✷"]
      ),
      :toggle => Dict(
        :frequency => 10,
        :frames => ["■", "□", "▪", "▫"]
      ),
      :balloon => Dict(
        :frequency => 10,
        :frames => ["", ".", "o", "O", "@", "*"]
      ),
      :balloon_2 => Dict(
        :frequency => 10,
        :frames => [".", "o", "O", "°", "O", "o", "."]
      ),
      :flip => Dict(
        :frequency => 10,
        :frames => ["-", "◡", "⊙", "-", "◠"]
      ),
      :burger => Dict(
        :frequency => 6,
        :frames => ["☱", "☲", "☴"]
      ),
      :dance => Dict(
        :frequency => 10,
        :frames => [">))'>", " >))'>", "  >))'>", "   >))'>", "    >))'>", "   <'((<", "  <'((<", " <'((<"]
      ),
      :shark => Dict(
        :frequency => 10,
        :frames => [
          "▐|\\____________▌",
          "▐_|\\___________▌",
          "▐__|\\__________▌",
          "▐___|\\_________▌",
          "▐____|\\________▌",
          "▐_____|\\_______▌",
          "▐______|\\______▌",
          "▐_______|\\_____▌",
          "▐________|\\____▌",
          "▐_________|\\___▌",
          "▐__________|\\__▌",
          "▐___________|\\_▌",
          "▐____________|\\▌",
          "▐____________/|▌",
          "▐___________/|_▌",
          "▐__________/|__▌",
          "▐_________/|___▌",
          "▐________/|____▌",
          "▐_______/|_____▌",
          "▐______/|______▌",
          "▐_____/|_______▌",
          "▐____/|________▌",
          "▐___/|_________▌",
          "▐__/|__________▌",
          "▐_/|___________▌",
          "▐/|____________▌"
        ]
      ),
      :pong => Dict(
        :frequency => 10,
        :frames => [
          "▐⠂       ▌",
          "▐⠈       ▌",
          "▐ ⠂      ▌",
          "▐ ⠠      ▌",
          "▐  ⡀     ▌",
          "▐  ⠠     ▌",
          "▐   ⠂    ▌",
          "▐   ⠈    ▌",
          "▐    ⠂   ▌",
          "▐    ⠠   ▌",
          "▐     ⡀  ▌",
          "▐     ⠠  ▌",
          "▐      ⠂ ▌",
          "▐      ⠈ ▌",
          "▐       ⠂▌",
          "▐       ⠠▌",
          "▐       ⡀▌",
          "▐      ⠠ ▌",
          "▐      ⠂ ▌",
          "▐     ⠈  ▌",
          "▐     ⠂  ▌",
          "▐    ⠠   ▌",
          "▐    ⡀   ▌",
          "▐   ⠠    ▌",
          "▐   ⠂    ▌",
          "▐  ⠈     ▌",
          "▐  ⠂     ▌",
          "▐ ⠠      ▌",
          "▐ ⡀      ▌",
          "▐⠠       ▌"
        ]
    )
)


text:
include("./imports.jl")

@testset ExtendedTestSet "DataAugmentation.jl" begin
    include("base.jl")
    include("buffered.jl")
    include("projective/base.jl")
    include("projective/bounds.jl")
    include("projective/compose.jl")
    include("projective/crop.jl")
    include("projective/affine.jl")
    include("items/image.jl")
    include("items/keypoints.jl")
    include("items/mask.jl")
    include("preprocessing.jl")
end


text:
# map operations to elements

#################################################
#
#   Generic functor based mapping
#
#################################################

# the macro to generate functor-based mapping functions
# for various number of arguments

macro compose_mapfuns(fname, AN)
    @assert AN != 0
    FN = AN > 0 ? AN : 1
    fname! = symbol(string(fname, '!'))
    _fname! = symbol(string('_', fname!))

    # code-gen preparation
    h = codegen_helper(AN)
    he = codegen_helper_ex(AN)
    aargs = h.args[2:end]
    ti = h.term(:i)
    t1 = h.term(1)
    mapker! = AN > 0 ? :vecmap! : :vecmapdiff!

    quote
        global $(fname)
        $(fname)($(h.dense_aparams...)) =
            $(fname!)(fun, Array($(h.termtype), $(h.inputsize)), $(aargs...))

        global $(fname!)
        function $(fname!)(fun::Functor{$FN}, dst::ContiguousArray, $(h.contiguous_aparams[2:end]...))
            n = length(dst)
            n == $(h.inputlen) || throw(DimensionMismatch("Inconsistent argument lengths."))
            for i = 1:n
                @inbounds dst[i] = $(ti)
            end
            return dst
        end

        function $(fname!)(fun::Functor{$FN}, dst::DenseArray, $(h.dense_aparams[2:end]...))
            siz = $(h.inputsize)
            size(dst) == siz || throw(DimensionMismatch("Inconsistent argument size."))
            $(_fname!)(siz, length(siz), dst, $(h.args...))
            return dst
        end

        global $(_fname!)
        function $(_fname!)(siz::Dims, dim::Int, dst::DenseArray, $(h.dense_aparams...))
            if dim == 1
                n = siz[1]::Int
                $(he.getstrides1)
                sdst1 = stride(dst, 1)::Int
                $(he.getoffsets)
                idst = offset(dst) + 1
                if $(he.contcol) && sdst1 == 1
                    $(mapker!)(n, parent(dst), idst, $(he.pkerargs...))
                else
                    $(mapker!)(n, parent(dst), idst, sdst1, $(he.pkerargs1...))
                end

            elseif dim == 2
                m = siz[1]::Int
                n = siz[2]::Int
                $(he.getstrides2)
                sdst1, sdst2 = strides(dst)::@compat(Tuple{Int, Int})
                $(he.getoffsets)
                idst = offset(dst) + 1
                if $(he.contcol) && sdst1 == 1
                    for j = 1:n
                        $(mapker!)(m, parent(dst), idst, $(he.pkerargs...))
                        $(he.nextcol)
                        idst += sdst2
                    end
                else
                    for j = 1:n
                        $(mapker!)(m, parent(dst), idst, sdst1, $(he.pkerargs1...))
                        $(he.nextcol)
                        idst += sdst2
                    end
                end

            else
                n = siz[dim]::Int
                for i = 1:n
                    $(_fname!)(siz, dim-1, ellipview(dst, i), $(he.eviewargs...))
                end
            end
        end
    end
end

@compose_mapfuns map 1
@compose_mapfuns map 2
@compose_mapfuns map 3
@compose_mapfuns mapdiff (-2)

map1!(f::Functor{1}, a::DenseArray) = map!(f, a, a)
map1!(f::Functor{2}, a::DenseArray, b::DenseArray) = map!(f, a, a, b)
map1!(f::Functor{2}, a::DenseArray, b::Number) = map!(f, a, a, b)
map1!(f::Functor{3}, a::DenseArray, b::DenseArray, c::DenseArray) = map!(f, a, a, b, c)
map1!(f::Functor{3}, a::DenseArray, b::DenseArray, c::Number) = map!(f, a, a, b, c)
map1!(f::Functor{3}, a::DenseArray, b::Number, c::DenseArray) = map!(f, a, a, b, c)
map1!(f::Functor{3}, a::DenseArray, b::Number, c::Number) = map!(f, a, a, b, c)

#################################################
#
#   Some inplace mapping functions
#
#################################################

macro derive_mapfuns1(fname, Op, BT)
    fname! = symbol(string(fname, '!'))
    quote
        global $(fname!)
        $(fname!){T<:$(BT)}(a::DenseArray{T}) = map1!($(Op)(), a)
        $(fname!){T<:$(BT)}(dst::DenseArray, a::DenseArray{T}) = map!($(Op)(), dst, a)
    end
end

macro derive_mapfuns2(fname, Op, BT)
    fname! = symbol(string(fname, '!'))
    quote
        global $(fname!)
        $(fname!){T1<:$(BT)}(a::DenseArray{T1}, b::$(BT)) = map1!($(Op)(), a, b)
        $(fname!){T1<:$(BT),T2<:$(BT)}(a::DenseArray{T1}, b::DenseArray{T2}) = map1!($(Op)(), a, b)
        $(fname!){T1<:$(BT)}(dst::DenseArray, a::DenseArray{T1}, b::$(BT)) = map!($(Op)(), dst, a, b)
        $(fname!){T1<:$(BT),T2<:$(BT)}(dst::DenseArray, a::DenseArray{T1}, b::DenseArray{T2}) = map!($(Op)(), dst, a, b)
    end
end

@derive_mapfuns1 negate Negate Number

@derive_mapfuns2 add Add Number
@derive_mapfuns2 subtract Subtract Number
@derive_mapfuns2 multiply Multiply Number
@derive_mapfuns2 divide Divide Number

@derive_mapfuns1 abs AbsFun Number
@derive_mapfuns1 abs2 Abs2Fun Number
@derive_mapfuns1 rcp RcpFun Number
@derive_mapfuns1 sqrt SqrtFun Real
@derive_mapfuns2 pow Pow Real

@derive_mapfuns1 floor FloorFun Real
@derive_mapfuns1 ceil CeilFun Real
@derive_mapfuns1 round RoundFun Real
@derive_mapfuns1 trunc TruncFun Real

@derive_mapfuns1 exp ExpFun Real
@derive_mapfuns1 log LogFun Real

# extensions

absdiff(x::DenseArrOrNum, y::DenseArrOrNum) = mapdiff(AbsFun(), x, y)
sqrdiff(x::DenseArrOrNum, y::DenseArrOrNum) = mapdiff(Abs2Fun(), x, y)

absdiff!(r::NumericArray, x::DenseArrOrNum, y::DenseArrOrNum) = mapdiff!(AbsFun(), r, x, y)
sqrdiff!(r::NumericArray, x::DenseArrOrNum, y::DenseArrOrNum) = mapdiff!(Abs2Fun(), r, x, y)

fma!(a::NumericArray, b::DenseArrOrNum, c::DenseArrOrNum) = map1!(FMA(), a, b, c)
fma(a::DenseArrOrNum, b::DenseArrOrNum, c::DenseArrOrNum) = map(FMA(), a, b, c)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module ResourceManagementClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_DeploymentOperationsApi.jl")
include("api_DeploymentsApi.jl")
include("api_ProvidersApi.jl")
include("api_ResourceGroupsApi.jl")
include("api_ResourcesApi.jl")
include("api_TagsApi.jl")

# export models
export convert, AliasPathType
export convert, AliasType
export convert, BasicDependency
export convert, DebugSetting
export convert, Dependency
export convert, Deployment
export convert, DeploymentExportResult
export convert, DeploymentExtended
export convert, DeploymentExtendedFilter
export convert, DeploymentListResult
export convert, DeploymentOperation
export convert, DeploymentOperationProperties
export convert, DeploymentOperationsListResult
export convert, DeploymentProperties
export convert, DeploymentPropertiesExtended
export convert, DeploymentValidateResult
export convert, ExportTemplateRequest
export convert, GenericResourceFilter
export convert, HttpMessage
export convert, Identity
export convert, ParametersLink
export convert, Plan
export convert, Provider
export convert, ProviderListResult
export convert, ProviderResourceType
export convert, Resource
export convert, ResourceGroup
export convert, ResourceGroupExportResult
export convert, ResourceGroupFilter
export convert, ResourceGroupListResult
export convert, ResourceGroupPatchable
export convert, ResourceGroupProperties
export convert, ResourceListResult
export convert, ResourceManagementErrorWithDetails
export convert, ResourceProviderOperationDisplayProperties
export convert, ResourcesMoveInfo
export convert, Sku
export convert, SubResource
export convert, TagCount
export convert, TagDetails
export convert, TagValue
export convert, TagsListResult
export convert, TargetResource
export convert, TemplateLink
export convert, GenericResource

# export operations
export convert, DeploymentOperationsApi, DeploymentsApi, ProvidersApi, ResourceGroupsApi, ResourcesApi, TagsApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
module ReinforcementLearningEnvironmentGym
using ReinforcementLearningBase
import ReinforcementLearningBase:interact!, reset!, getstate, plotenv, actionspace
using PyCall
const gym = PyNULL()

function __init__()
    copy!(gym, pyimport_conda("gym", ""))
    pyimport("pybullet_envs")
end

function gymspace2jlspace(s::PyObject)
    spacetype = s[:__class__][:__name__]
    if     spacetype == "Box"           BoxSpace(s[:low], s[:high])
    elseif spacetype == "Discrete"      DiscreteSpace(s[:n], 1)
    elseif spacetype == "MultiBinary"   MultiBinarySpace(s[:n])
    elseif spacetype == "MultiDiscrete" MultiDiscreteSpace(s[:nvec], 1)
    elseif spacetype == "Tuple"         map(gymspace2jlspace, s[:spaces])
    elseif spacetype == "Dict"          Dict(map((k, v) -> (k, gymspace2jlspace(v)), s[:spaces]))
    else error("Don't know how to convert [$(spacetype)]")
    end
end

struct GymEnv{Ta<:AbstractSpace, To<:AbstractSpace} <: AbstractEnv
    pyobj::PyObject
    observationspace::To
    actionspace::Ta
    state::PyObject
end

function GymEnv(name::String)
    pyenv = gym[:make](name)
    obsspace = gymspace2jlspace(pyenv[:observation_space])
    actspace = gymspace2jlspace(pyenv[:action_space])
    state = PyNULL()
    env = GymEnv(pyenv, obsspace, actspace, state)
    reset!(env) # state needs to be set to call defaultbuffer in RL
    env
end

function interact!(env::GymEnv, action)
    pycall!(env.state, env.pyobj[:step], PyVector, action)
    (observation=env.state[1], reward=env.state[2], isdone=env.state[3])
end

function interact!(env::GymEnv{DiscreteSpace}, action::Int)
    pycall!(env.state, env.pyobj[:step], PyVector, action - 1)
    (observation=env.state[1], reward=env.state[2], isdone=env.state[3])
end

function interact!(env::GymEnv{MultiDiscreteSpace}, action::AbstractArray{Int})
    pycall!(env.state, env.pyobj[:step], PyVector, action .- 1)
    (observation=env.state[1], reward=env.state[2], isdone=env.state[3])
end

"Not very useful, kept for compat"
function getstate(env::GymEnv)
    if pyisinstance(env.state, PyCall.@pyglobalobj :PyTuple_Type)
        (observation=env.state[1], isdone=env.state[3])
    else
        # env has just been reseted
        (observation=Float64.(env.state), isdone=false)
    end
end

reset!(env::GymEnv) = (observation=Float64.(pycall!(env.state, env.pyobj[:reset], PyArray)),)
plotenv(env::GymEnv) = env.pyobj[:render]()
actionspace(env::GymEnv) = env.actionspace

"""
    listallenvs(pattern = r"")

List all registered gym environment names. The optional argument `pattern`
allows to list all environment  that contain the `pattern` in their name.
"""
function listallenvs(pattern = r"")
    envs = sort(py"[spec.id for spec in $gym.envs.registry.all()]")
    if pattern != ""
        envs[findall(x -> occursin(pattern, x), envs)]
    else
        envs
    end
end

export GymEnv, listallenvs
end # module


text:
@testset "Testing Reco submodule" begin

fnRecoV1 = "./data/mdf/reco_V1.mdf"
fnRecoV2 = "./data/mdf/reco_V2.mdf"
fnRecoV2b = "./data/mdf/reco_V2b.mdf"

mdfv1 = MPIFile(fnRecoV1)
@test typeof(mdfv1) <: MDFFileV1

saveasMDF(fnRecoV2, fnRecoV1)

mdfv2 = MPIFile(fnRecoV2)
@test typeof(mdfv2) <: MDFFileV2

c1 = loadRecoData(fnRecoV1)
c2 = loadRecoData(fnRecoV2)
saveRecoData(fnRecoV2b,c1)

mdfv2b = MPIFile(fnRecoV2b)
@test typeof(mdfv2b) <: MDFFileV2

for mdf in (mdfv1,mdfv2,mdfv2b)
  @info "Test $mdf"
  @test recoSize(mdf) == [40,40,1]
  @test norm(recoFovCenter(mdf) - [0.0,0.0,0.0]) < eps() # Image conversion fails
  @test recoFov(mdf) == [0.044,0.044,0.001]
  @test recoOrder(mdf) == "xyz"
  @test recoPositions(mdf) == nothing

  @test size( recoData(mdf) ) == (1, 1600, 1)
end

@test c1 == c2

end


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

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

# ╔═╡ b6ff5052-04ec-11eb-043b-458c42b19edb
begin
    import Pkg
    # activate the shared project environment
    Pkg.activate(Base.current_project())

	using Distributions
	using Plots
	using PlutoUI
	using Random
end

# ╔═╡ 87c1dc3b-d361-4f46-a68a-c28d7c7de0c3
md"""
# Equilibrium binding simulations

This notebook can either:

1. simulate the outcome of an equilibrium binding experiment given a hypothetical $K_D$ and desired experimental parameters (section [Simulate a binding curve](#d50721c7-34e4-411d-ae6b-dc6196b0f5a9)) or
2. calculate the fractional saturation of receptor given a $K_D$ and total concentrations of receptor and ligand (section [Predict saturation of receptor](#3954adfa-73cc-495a-93c1-bcde87a60ce8)).

The following reference explains very well the theory of equilibrium binding experiments, as well as many important practical considerations:

> Jarmoskaite I, AlSadhan I, Vaidyanathan PP & Herschlag D (2020) How to measure and evaluate binding affinities. *eLife* **9**: e57264 <https://doi.org/10.7554/eLife.57264>
"""

# ╔═╡ d50721c7-34e4-411d-ae6b-dc6196b0f5a9
md"""
## Simulate a binding curve

Enter values of $K_D$ and experimental parameters below to simulate the resulting data points (sampling the binding curve) as they would come out if you performed the experiment. The vertical line on the graph shows the value of $K_D$. The horizontal lines show zero and full saturation of the receptor.
"""

# ╔═╡ c4d66e7a-38e4-497d-b680-8c50c51fe146
md"""
Value of $K_D$: $(@bind kd PlutoUI.Slider(0.1:0.1:100.0, default = 10.0, show_value = true))

Concentration of ligand at highest titration point $L_{max}$: $(@bind lmax PlutoUI.Slider(kd:1.0:100000.0, default = 1000.0, show_value = true))

Dilution factor between titration points: $(@bind dilution PlutoUI.NumberField(1.1:0.1:10.0, default = 2.0))

Number of titration points: $(@bind points PlutoUI.NumberField(1:1:24, default = 15))

Concentration of receptor $R_{tot}$: $(@bind Rtot PlutoUI.Slider(0.001:1.0:50.0, default = 1.0, show_value = true))

Simulate random, normally distributed noise? $(@bind simulateNoise PlutoUI.CheckBox(default = false))
"""

# ╔═╡ b05df20b-353e-462d-aa13-534862f8f6cd
if simulateNoise
	md"""
	Mean (μ) and standard deviation (σ)

	μ = $(@bind μ PlutoUI.NumberField(-3.0:0.1:3.0, default = 0.0))
	σ = $(@bind σ PlutoUI.NumberField(0.0:0.001:0.1, default = 0.001))
	"""
end

# ╔═╡ d49dee36-c12d-4206-9699-a181570caa80
md"""
Ratio $\frac{L_{max}}{K_D} =$ $(lmax / kd)

Ratio $\frac{K_D}{R_{tot}} =$ $(kd / Rtot)
"""

# ╔═╡ 915b2a49-cb32-40a0-bfae-21380144ed3a
if simulateNoise
	md"""
	Signal-to-noise ratio $\frac{S}{σ} =$ $(1.0 / σ)
	"""
end

# ╔═╡ 3954adfa-73cc-495a-93c1-bcde87a60ce8
md"""
## Predict saturation of receptor

Value of $K_D$: $(@bind kd2 PlutoUI.NumberField(1.0:10.0:1000.0, default = 10.0))

Concentration of ligand $L_{tot}$: $(@bind ltot2 PlutoUI.NumberField(1.0:100.0:100000.0, default = 1000.0))

Concentration of receptor $R_{tot}$: $(@bind rtot2 PlutoUI.NumberField(0.001:0.01:100000.0, default = 1000.0))

Percent saturation of receptor $B =$
"""

# ╔═╡ ae1d6dac-92d7-4962-a641-63965f73a8f2
md"""
## Code

The quadratic model gives $B$, the fraction of bound receptor (or fractional saturation of receptor) for any values of receptor concentration $R_{tot}$ and ligand concentration $L_{tot}$, given a binding affinity determined by the value of $K_D$:

$B = \frac{(K_{D} + R_{tot} + L_{tot}) - \sqrt{(- K_{D} - R_{tot} - L_{tot})^2 - 4 \times R_{tot} \times L_{tot}}}{2 \times R_{tot}}$

The function below defines this binding model in code:
"""

# ╔═╡ 1eed6126-331b-4218-acc7-8b7c132ef355
function quadratic(Ltot)
	( (kd + Rtot + Ltot) - sqrt((- kd - Rtot - Ltot) ^ 2 - 4 * Rtot * Ltot) ) / (2 * Rtot)
end

# ╔═╡ a5ac735b-f748-4d37-9cf1-7e64e1ea4c53
function quadratic(kd, Rtot, Ltot)
	( (kd + Rtot + Ltot) - sqrt((kd + Rtot + Ltot) ^ 2 - 4 * Rtot * Ltot) ) / (2 * Rtot)
end

# ╔═╡ 4d816572-03a3-481c-a29c-522f0535a679
round(100 * quadratic(kd2, ltot2, rtot2), digits = 2)

# ╔═╡ 3794717b-1f61-4622-9a55-ad86ab947831
md"""
The code below calculates concentrations of the titration series ($X$ values for the plot) given the maximal ligand concentration $L_{max}$, the dilution factor and the number of titration points entered above:
"""

# ╔═╡ 7c43edeb-a5e6-4fb4-8494-e2072d0d56fd
x = [ lmax / dilution^i for i in 0:points-1 ]

# ╔═╡ 06f341e8-6d40-4eb7-b7e9-58c5d45f1f85
md"""
The code below uses the quadratic model to calculate fractional saturation of the receptor ($Y$ values for the plot) at each titration point for the values of $K_D$ and $R_{tot}$ entered above:
"""

# ╔═╡ 82af41ba-d9b9-4ce2-8736-6fe32c8e024e
if simulateNoise
	noise = rand(Normal(μ, σ), points)
	y = quadratic.(x) .+ noise
else
	y = quadratic.(x)
end

# ╔═╡ efc70ae1-abcc-472d-8df8-c55cbbd6529b
begin
	scatter(
		x,
		y,
		xscale = :log10,
		legend = :none,
		xlabel = "Ligand concentration",
		ylabel = "Predicted fractional saturation of receptor"
	)
	vline!([kd], color = :red)
	hline!([1.0], color = :black)
	hline!([0.0], color = :black)
end

# ╔═╡ 9c5021f3-67ab-464d-8f22-25e269c36083
md"""
## Required packages
"""

# ╔═╡ b0340ca3-02d5-49a9-841b-7827e4aee2bb
PlutoUI.TableOfContents()

# ╔═╡ Cell order:
# ╟─87c1dc3b-d361-4f46-a68a-c28d7c7de0c3
# ╟─d50721c7-34e4-411d-ae6b-dc6196b0f5a9
# ╟─c4d66e7a-38e4-497d-b680-8c50c51fe146
# ╟─b05df20b-353e-462d-aa13-534862f8f6cd
# ╟─d49dee36-c12d-4206-9699-a181570caa80
# ╟─915b2a49-cb32-40a0-bfae-21380144ed3a
# ╟─efc70ae1-abcc-472d-8df8-c55cbbd6529b
# ╟─3954adfa-73cc-495a-93c1-bcde87a60ce8
# ╠═4d816572-03a3-481c-a29c-522f0535a679
# ╟─ae1d6dac-92d7-4962-a641-63965f73a8f2
# ╠═1eed6126-331b-4218-acc7-8b7c132ef355
# ╠═a5ac735b-f748-4d37-9cf1-7e64e1ea4c53
# ╟─3794717b-1f61-4622-9a55-ad86ab947831
# ╠═7c43edeb-a5e6-4fb4-8494-e2072d0d56fd
# ╟─06f341e8-6d40-4eb7-b7e9-58c5d45f1f85
# ╠═82af41ba-d9b9-4ce2-8736-6fe32c8e024e
# ╟─9c5021f3-67ab-464d-8f22-25e269c36083
# ╠═b6ff5052-04ec-11eb-043b-458c42b19edb
# ╠═b0340ca3-02d5-49a9-841b-7827e4aee2bb

max_stars_repo_name: Guillawme/julia-curve-fitting

max_issues_repo_name: Guillawme/julia-curve-fitting

max_forks_repo_name: Guillawme/julia-curve-fitting


text:
const REQUIRED_ARGS = ["resistance_file",
                       "radius",
                       "project_name"]

const TRUELIST = ["true", "True", "TRUE"]

const SINGLE = ["single", "Single", "SINGLE"]


text:
export constraint_learning_by_obseration

"""
Compute matrix/image properties and return the results in a dictionary.
m is a tensor where the first index is the image index.
i.e., m[i,:,:] is one image
"""
function constraint_learning_by_obseration(comp_grid,m_train::Array{TF}) where {TF<:Real}

if TF==Float64
  TI=Int64
else
  TI=Int32
end

if length(size(m_train))==3
  (n_train_ex,t2,t3)=size(m_train)
else
  (t2,t3)=size(m_train)
  n_train_ex=1
end

#a non exhaustive list of things to learn/observe. Add some more...
observations=Dict()

observations["nuclear_norm"] = zeros(TF,n_train_ex)
observations["nuclear_Dx"]   = zeros(TF,n_train_ex)
observations["nuclear_Dz"]   = zeros(TF,n_train_ex)
observations["rank_095"]     = zeros(TI,n_train_ex)
observations["TV"]           = zeros(TF,n_train_ex)
observations["wavelet_l1"]   = zeros(TF,n_train_ex)
observations["Dx_l1"]        = zeros(TF,n_train_ex)
observations["Dz_l1"]        = zeros(TF,n_train_ex)
#observations["curvelet_l1"]  = zeros(TF,n_train_ex)
observations["DFT_l1"]       = zeros(TF,n_train_ex)
observations["DFT_card_095"] = zeros(TI,n_train_ex)
#observations["curvelet_card_095"]= zeros(TI,n_train_ex)
observations["TV_card_095"]  = zeros(TI,n_train_ex)
observations["annulus"]      = zeros(TF,n_train_ex)
observations["TV_annulus"]   = zeros(TF,n_train_ex)
observations["D_l2"]         = zeros(TF,n_train_ex)
observations["D_x_min"]      = zeros(TF,n_train_ex)
observations["D_x_max"]      = zeros(TF,n_train_ex)
observations["D_z_min"]      = zeros(TF,n_train_ex)
observations["D_z_max"]      = zeros(TF,n_train_ex)

observations["DCT_x_LB"]      = zeros(TF,t2) .+ 1e8
observations["DCT_x_UB"]      = zeros(TF,t2)

observations["DCT_y_LB"]      = zeros(TF,t3) .+ 1e8 #add because we take the min looping over every training image,
observations["DCT_y_UB"]      = zeros(TF,t3)

observations["hist_min"]      = zeros(TF,t2*t3) .+ 1e8
observations["hist_max"]      = zeros(TF,t2*t3)

observations["hist_TV_min"]   = zeros(TF,(t2-1)*t3+t2*(t3-1)) .+ 1e8
observations["hist_TV_max"]   = zeros(TF,(t2-1)*t3+t2*(t3-1))

#get transform-domain operators
#need to install CurveLab to use the Curvelet transform
(Dx_OP, dummy1, dummy2, Dx_TD_n) = get_TD_operator(comp_grid,"D_x",TF)
(Dz_OP, dummy1, dummy2, Dz_TD_n) = get_TD_operator(comp_grid,"D_z",TF)
(TV_OP, dummy1, dummy2, TD_n_TV) = get_TD_operator(comp_grid,"TV",TF)
#(C_OP, dummy1, dummy2, dummy3)   = get_TD_operator(comp_grid,"curvelet",TF)
(DFT, dummy1, dummy2, dummy3)    = get_TD_operator(comp_grid,"DFT",TF)
if comp_grid.n[1]==comp_grid.n[2]
  (DWT, dummy1, dummy2, dummy3)    = get_TD_operator(comp_grid,"wavelet",TF)
end
#DFT = x-> vec(fft(reshape(x,comp_grid.n)))
DCT = x-> vec(dct(reshape(x,comp_grid.n)))

DCT_x = x-> vec(dct(reshape(x,comp_grid.n),1))
DCT_y = x-> vec(dct(reshape(x,comp_grid.n),2))

for i=1:n_train_ex #can be changed to a parallel loop for larger datasets

    training_image = vec(m_train[i,:,:])

    #observe sorted values for mean histogram
    observations["hist_min"].=min.(observations["hist_min"],sort(training_image))
    observations["hist_max"].=max.(observations["hist_max"],sort(training_image))

    #observe sorted values for mean histogram
    observations["hist_TV_min"].=min.(observations["hist_TV_min"],sort(TV_OP*training_image))
    observations["hist_TV_max"].=max.(observations["hist_TV_max"],sort(TV_OP*training_image))

    #observe Nuclear norm
    sv=svdvals(reshape(training_image,comp_grid.n))
    observations["nuclear_norm"][i]=norm(sv,1)

    #observe nuclear norm in transform-domain (TV-operator, total-nuclear-variation)
    observations["nuclear_Dx"][i]=norm(svdvals(reshape(Dx_OP*training_image,Dx_TD_n)),1)
    observations["nuclear_Dz"][i]=norm(svdvals(reshape(Dz_OP*training_image,Dz_TD_n)),1)

    #observe rank that represents 95% of each image
    sv=cumsum(sv)
    sv.=sv./maximum(sv)
    observations["rank_095"][i]=findfirst(sv.>0.95)

    #observe slope in x direction
    observations["D_x_min"][i]=minimum(Dx_OP*training_image)
    observations["D_x_max"][i]=maximum(Dx_OP*training_image)

    #observe slope in z direction
    observations["D_z_min"][i]=minimum(Dz_OP*training_image)
    observations["D_z_max"][i]=maximum(Dz_OP*training_image)

    #observe anisotropic total-variation
    observations["TV"][i]=norm(TV_OP*training_image,1)

    #Wavelet l1
    if comp_grid.n[1]==comp_grid.n[2]
      observations["wavelet_l1"][i]=norm(DWT*training_image,1)
    end

    #different direction separately
    observations["Dx_l1"][i] = norm(Dx_OP*training_image,1)
    observations["Dz_l1"][i] = norm(Dz_OP*training_image,1)

    #observe l2 norm of discrete gradients (i.e., l2 version of anisotripic TV)
    observations["D_l2"][i]=norm(TV_OP*training_image,2)

    #observe annulus
    observations["annulus"][i]=norm(vec(training_image),2)
    observations["TV_annulus"][i]=norm(TV_OP*vec(training_image),2)

    #observe curvelet domain l1 norm
    #observations["curvelet_l1"][i]=norm(C_OP*training_image,1)

    #observe Fourier domain l1 norm
    observations["DFT_l1"][i]=norm(DFT*training_image,1)

    #observe nr or required DFT atoms to represent 95%
    temp=DFT*training_image; temp=abs.(temp); temp=sort(temp); temp=cumsum(temp); temp=temp./maximum(temp);
    observations["DFT_card_095"][i]=length(temp)-findfirst(temp.>0.05)

    #observe nr or required curvelet atoms to represent95%
    # temp=C_OP*training_image; temp=abs.(temp); temp=sort(temp); temp=cumsum(temp); temp=temp./maximum(temp);
    # observations["curvelet_card_095"][i]=length(temp)-findfirst(temp.>0.05)

    #observe nr or required TV atoms to represent95%
    temp=TV_OP*training_image; temp=abs.(temp); temp=sort(temp); temp=cumsum(temp); temp=temp./maximum(temp);
    observations["TV_card_095"][i]=length(temp)-findfirst(temp.>0.05)

    #observe DCT spectrum per column (x)
    temp_img=dct(reshape(training_image,t2,t3),1);
    (temp_vals,indx)=findmin(temp_img,dims=2)
    observations["DCT_x_LB"].=min.(temp_vals[:],observations["DCT_x_LB"])
    (temp_vals,indx)=findmax(temp_img,dims=2)
    observations["DCT_x_UB"].=max.(temp_vals[:],observations["DCT_x_UB"])

    #observe DCT spectrum per row (y)
    temp_img=dct(reshape(training_image,t2,t3),2);
    (temp_vals,indx)=findmin(temp_img,dims=1)
    observations["DCT_y_LB"].=min.(temp_vals[:],observations["DCT_y_LB"])
    (temp_vals,indx)=findmax(temp_img,dims=1)
    observations["DCT_y_UB"].=max.(temp_vals[:],observations["DCT_y_UB"])
end

return observations
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
@testset "NaturalPolicyGradient (PointMass)" begin
    tseed!(1)
    etype = LyceumMuJoCo.PointMass

    e = etype()
    dobs, dact = length(observationspace(e)), length(actionspace(e))

    DT = Float32
    Hmax, K = 256, 16
    N = Hmax * K

    policy = DiagGaussianPolicy(
        multilayer_perceptron(dobs, 16, 16, dact, σ=tanh),
        zeros(dact)
    )
    policy = Flux.paramtype(DT, policy)

    value = multilayer_perceptron(dobs, 32, 32, 1, σ=Flux.relu)
    valueloss(bl, X, Y) = mse(vec(bl(X)), vec(Y))

    valuetrainer = FluxTrainer(
        optimiser = RADAM(1e-3),
        szbatch = 64,
        lossfn = valueloss,
        stopcb = s->s.nepochs > 1
    )
    value = Flux.paramtype(DT, value)

    npg = NaturalPolicyGradient(
        n -> tconstruct(etype, n),
        policy,
        value,
        gamma = 0.97,
        gaelambda = 0.95,
        valuetrainer,
        Hmax=Hmax,
        norm_step_size=0.01,
        N=N,
    )

    envsampler = EnvironmentSampler(n -> tconstruct(etype, n))

    # Test that the terminal reward for the mean policy is > 0.95 for at least 5
    # iterations in a row, in at most 250 iterations.
    npasses = 0
    for (i, state) in enumerate(npg)
        (npasses > 5 || i > 250) && break
        batch = sample(envsampler, N, reset! = randreset!, Hmax=Hmax) do a, o
            a .= policy(o)
        end
        npasses = mean(τ -> τ.R[end], batch) > 0.95 ? npasses + 1 : 0
    end
    @test npasses > 5
end


text:
using NPZ
using CSV
using Dates
using DataFrames
using Statistics
using CartesianGeneticProgramming
using Cambrian
import Cambrian.mutate
using Plots
using StatsPlots
include("utils.jl")

struct Month
    topo0
    topo1
    d_topo

    tp
    tp_mean
    tp_std
    tp_max
    tp_med

    hs
    hs_mean
    hs_std
    hs_max
    hs_med

    direction
    dir_mean
    dir_std
    dir_max
    dir_med

    input_row
end

dataset_path = get_dataset_path()
n_months = 236
months = []

for i in 1:n_months
    topo_path = string(dataset_path, '/', i, "_topo.npy")
    topo = npzread(topo_path)
    topo0 = norm_to_range(topo[1, :], -10.0, 10.0, -1, 1)
    topo1 = norm_to_range(topo[2, :], -10.0, 10.0, -1, 1)

    waves_path = string(dataset_path, '/', i, "_waves.csv")
    waves_csv = CSV.read(waves_path, DataFrame)

    # t = min_max_norm(waves_csv.t)
    tp_month = norm_to_range(waves_csv.tp, 0, 25, -1, 1)
    hs_month = norm_to_range(waves_csv.hs, 0, 10, -1, 1)
    direction_month = norm_to_range(waves_csv.direction, 0, 360, -1, 1)

    # plot_tp = density(tp_month, label="Tp")
    # plot_hs = density(hs_month, label="Hs")
    # plot_dir = density(direction_month, label="Direction")
    # display(plot(plot_tp, plot_hs, plot_dir, layout = (3, 1)))
    # readline()

    input_row = vcat(topo0,
        mean(tp_month), std(tp_month), maximum(tp_month), median(tp_month),
        mean(hs_month), std(hs_month), maximum(hs_month), median(hs_month),
        mean(direction_month), std(direction_month), maximum(direction_month), median(direction_month))

    month = Month(topo0, topo1, topo1 - topo0,
        tp_month, mean(tp_month), std(tp_month), maximum(tp_month), median(tp_month),
        hs_month, mean(hs_month), std(hs_month), maximum(hs_month), median(hs_month),
        direction_month, mean(direction_month), std(direction_month), maximum(direction_month), median(direction_month),
        input_row)
    push!(months, month)
end

function evaluate(ind::CGPInd, months::Vector{Any})
    errors = Vector{Float64}(undef, size(months)[1])

    for m_i in 1:size(months)[1]
        month = months[m_i]
        d_topo = process(ind, month.input_row)
        errors[m_i] = mean(abs.(d_topo - month.d_topo))
    end
    [mean(errors) * -1]
end

cfg = get_config("topo.yaml")
fit(i::CGPInd) = evaluate(i, months)
mutate(i::CGPInd) = goldman_mutate(cfg, i)
e = CGPEvolution(cfg, fit)
run!(e)

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
function K_DKGQ(elm::Quad)::Matrix{Float64}
    return K_GQ12(elm)+K_DKQ(elm)
end


text:
using Test, YaoPlots
using ZXCalculus, Multigraphs

@testset "zx plot" begin
    g = Multigraph(6)
    for e in [[1,3],[2,3],[3,4],[4,5],[4,6]]
        add_edge!(g, e)
    end
    ps = [0, 0, 0//1, 2//1, 0, 0]
    v_t = [SpiderType.In, SpiderType.Out, SpiderType.X, SpiderType.Z, SpiderType.Out, SpiderType.In]
    zxd = ZXDiagram(g, v_t, ps)
    plot(zxd)
    replace!(Rule{:b}(), zxd)
    plt = plot(zxd)
    @test plt !== nothing
end


text:
using JSONSchemaGenerator

using JSON3
using JSONSchema
using OrderedCollections
using Test
using StructTypes

const INT_EXAMPLE = """
{
  "\$schema": "http://json-schema.org/draft-07/schema#",

  "type": "integer"
}
"""

# TODO: Test other types

struct Address
    street_address::String
    city::String
    state::String
end
StructTypes.StructType(::Type{Address}) = StructTypes.Struct()

struct Model
    billing_address::Address
    shipping_address::Address
end
StructTypes.StructType(::Type{Model}) = StructTypes.Struct()

const MODEL_EXAMPLE = """
{
  "\$schema": "http://json-schema.org/draft-07/schema#",

  "definitions": {
    "address": {
      "type": "object",
      "properties": {
        "street_address": { "type": "string" },
        "city":           { "type": "string" },
        "state":          { "type": "string" }
      },
      "required": ["street_address", "city", "state"],
      "additionalProperties": false
    }
  },

  "type": "object",

  "properties": {
    "billing_address": { "\$ref": "#/definitions/address" },
    "shipping_address": { "\$ref": "#/definitions/address" }
  },

  "required": ["billing_address", "shipping_address"],
  "additionalProperties": false
}
"""

abstract type AbstractFoo end
StructTypes.StructType(::Type{AbstractFoo}) = StructTypes.AbstractType()
StructTypes.subtypes(::Type{AbstractFoo}) = (foo=Foo, bar=Bar)

struct Foo <: AbstractFoo
    val::Int
end
StructTypes.StructType(::Type{Foo}) = StructTypes.Struct()

struct Bar <: AbstractFoo
    val::Float64
end
StructTypes.StructType(::Type{Bar}) = StructTypes.Struct()

struct Option
    val::Union{Nothing,Int}
end
StructTypes.StructType(::Type{Option}) = StructTypes.Struct()

struct Inner
    val::Float64
end
StructTypes.StructType(::Type{Inner}) = StructTypes.Struct()

struct Middle
    val::Inner
end
StructTypes.StructType(::Type{Middle}) = StructTypes.Struct()

struct Outer
    val::Middle
end
StructTypes.StructType(::Type{Outer}) = StructTypes.Struct()

const NESTED_EXAMPLE = """
{
  "\$schema": "http://json-schema.org/draft-07/schema#",

  "definitions": {
    "inner": {
      "type": "object",
      "properties": {
        "val": { "type": "number" }
      },
      "required": ["val"],
      "additionalProperties": false
    },
    "middle": {
      "type": "object",
      "properties": {
        "val": { "\$ref": "#/definitions/inner" }
      },
      "required": ["val"],
      "additionalProperties": false
    }
  },

  "type": "object",

  "properties": {
    "val": { "\$ref": "#/definitions/middle" }
  },

  "required": ["val"],
  "additionalProperties": false
}
"""

stripws(str) = replace(str, r"\s" => "")

@testset "JSONSchemaGenerator.jl" begin
    @testset "Schema" begin
        @test JSONSchemaGenerator.isrequired(Address)
        @test !JSONSchemaGenerator.isrequired(Union{Nothing,Address})
        @test JSONSchemaGenerator.Schema(String) == JSONSchemaGenerator.JSONString()
        @test JSONSchemaGenerator.Schema(Int) == JSONSchemaGenerator.JSONInt()
        @test JSONSchemaGenerator.Schema(Float64) == JSONSchemaGenerator.JSONNumber()
        @test JSONSchemaGenerator.Schema(Bool) == JSONSchemaGenerator.JSONBool()
        @test JSONSchemaGenerator.Schema(Vector{Float64}) ==
              JSONSchemaGenerator.JSONList(JSONSchemaGenerator.JSONNumber())
        @test JSONSchemaGenerator.Schema(typeof((1, 2, 3, 4))) ==
              JSONSchemaGenerator.JSONList(JSONSchemaGenerator.JSONInt(), 4)
        @test JSONSchemaGenerator.Schema(typeof(("tuple", 1))) == JSONSchemaGenerator.JSONTuple(
            [JSONSchemaGenerator.JSONString(), JSONSchemaGenerator.JSONInt()], 2
        )
        @test JSONSchemaGenerator.Schema(Address) == JSONSchemaGenerator.JSONObject(
            :Address,
            OrderedDict(
                :street_address => JSONSchemaGenerator.JSONString(),
                :city => JSONSchemaGenerator.JSONString(),
                :state => JSONSchemaGenerator.JSONString(),
            ),
            [:street_address, :city, :state],
        )
        @test JSONSchemaGenerator.Schema(AbstractFoo) == JSONSchemaGenerator.JSONOneOf([
            JSONSchemaGenerator.Schema(Foo), JSONSchemaGenerator.Schema(Bar)
        ])
        @test JSONSchemaGenerator.Schema(Option) == JSONSchemaGenerator.JSONObject(
            :Option, OrderedDict(:val => JSONSchemaGenerator.JSONInt()), []
        )
        @test string(JSONSchemaGenerator.JSONRef("blob")) ==
              "{\"\$ref\":\"#/definitions/blob\"}"
        @test generate(SchemaGenerator(Int)) == stripws(INT_EXAMPLE)
        @test generate(SchemaGenerator(Model)) == stripws(MODEL_EXAMPLE)
        @test generate(SchemaGenerator(Outer)) == stripws(NESTED_EXAMPLE)
    end
    @testset "Round Trip" begin
        schema_file = tempname()
        json_file = tempname()
        open(schema_file, "w") do f
            write(f, generate(SchemaGenerator(Model)))
        end
        model = Model(
            Address("1600 Pennsylvania Avenue NW", "Washington", "DC"),
            Address("1st Street SE", "Washington", "DC"),
        )
        open(json_file, "w") do f
            JSON3.write(f, model)
        end
        schema = Schema(read(schema_file, String))
        @test isvalid(JSON3.read(read(json_file, String), Dict), schema)
    end
end


text:
using Test

@testset "Some tests of our rate band and indirect tax calcs" begin

    # TODO

    calc_indirect(
           expenditure=100.0,
           selling_price=10,
           vat=0.20,
           advalorem = 0.4,
           specific = 2.0 )

    inc = 100.0
    rates :: RateBands = [0.2,0.4]
    bands :: RateBands = [50.0,99999]
    onethresh :: RateBands = [999.0]
    tr = calctaxdue( taxable=inc, rates=rates, thresholds=bands )
    @test tr.due ≈ 30.0
    tr = calctaxdue( taxable=inc, rates=RateBands([1.0]), thresholds=onethresh )
    @test tr.due ≈ 100.0

    # tests of length of bands
    tr = calctaxdue(taxable=10.0,rates=RateBands([0.12,0.21]), thresholds=RateBands([12.0,99999999999999]))
    @test tr.due ≈ 1.2

    tr = calctaxdue(taxable=10.0,rates=RateBands([0.12,0.21]), thresholds=RateBands([12.0]))
    @test tr.due ≈ 1.2

    tr = calctaxdue(taxable=10.0,rates=RateBands([0.12]), thresholds=RateBands(zeros(0)))
    @test tr.due ≈ 1.2

end # testset


text:
include("magpoly.jl")

using .magpoly
const mp = magpoly
using DelimitedFiles, Random, Distributions, LinearAlgebra
#= This file contains the structures and functions to do carry
   out the inference using synthetic likelihood =#

struct summary_stats{T1<:Int, T2<:Float64}
    n_samples::T1
    magnetizations::T2
    spin_feature_overlaps::T2
    lag_p_autocorrelation::T2
end

function summary_stats(n_samples::Int)
    summary_stats(n_samples,zeros(n_samples),zeros(n_samples),zeros(n_samples))
end

#####################################################################################
# This is the structure for the parameter samples generated through synthetic likelihood
struct SL_data{T1<:Int, T2<:Float64}
    n_samples::T1
    n_params::T1
    param_values::T2
end

function SL_data(n_samples::Int, n_params::Int)
    SL_data(n_samples, n_params, zeros(n_params, n_samples))
end

#####################################################################################
#####################################################################################
function lag_p_autocovariance(vect::Array{Int}, p::Int)
    avg = 0.0
    avg_lag = 0.0
    var = 0.0
    var_lag = 0.0
    corr = 0.0
    n = length(vect)
    for i in 1:(n-p)
        avg += vect[i]
        avg_lag += vect[i+p]
        corr += vect[i]*vect[i+p]
    end
    avg = avg/(n-p)
    avg_lag = avg_lag/(n-p)
    corr = corr/(n-p)
    for i in 1:(n-p)
        var += (vect[i]-avg)^2
        var_lag += (vect[i+p]-avg_lag)^2
    end
    var = var/(n-p-1)
    var_lag = var_lag/(n-p-1)
    return (corr - avg*avg_lag)/(sqrt(var*var_lag))
end

function compute_summary_stats!(ss::Matrix{Float64}, sp_conf::Matrix{Int}, feats::Matrix{Float64})
    n_s = size(ss,2) #number of samples
    n_m = size(sp_conf,1)
    n_feats = size(feats, 2)
    for i_sample in 1:n_s
        for i_feat in 1:n_feats
            ss[i_feat,i_sample]=0
        end
        for i_m in 1:n_m
            for i_feat in 1:n_feats
                ss[i_feat,i_sample] += sp_conf[i_m,i_sample]*feats[i_m, i_feat]
            end
        end
        for i_feat in 1:n_feats
            ss[i_feat,i_sample] = ss[i_feat,i_sample]/n_m
        end
    end
end

function log_synth_likelihood(mat::Matrix{Float64}, avg::Array{Float64}, data::Array{Float64})
    s = -0.5*(avg .- data)'*inv(mat)*(avg .- data) -0.5*log(abs(det(mat)))
    return s
end

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

## for now this is just the backbone of the function to carry out the likelihood-free inference task
function synthetic_likelihood_polymer(n_samples::Int, sample_lag::Int, stride::Int, n_params::Int, delta_w::Float64,initial_weights::Array{Float64}, features_file::String, data_file::String)
    io = open(features_file, "r")
    features = readdlm(io,Float64)
    close(io)
    n_mono = size(features, 1)
    n_feats = size(features, 2)

    io = open(data_file, "r")
    data = readdlm(io,Float64; header=true)[1][1:end,:]
    close(io)
    println(data,"\n\n")
    n_data = size(data, 1)
    n_ss = size(data, 2) # number of summary statistics
    println(n_ss, "\n\n")

    accepted_moves = 0
    n_strides = cld(n_samples*sample_lag, stride) # integer ceiling of the division
    spins_coupling = 1.0
    #spins_coupling = 0.0
    alpha = 0.5
    inv_temps = [1.0, 0.9, 0.8, 0.74, 0.65, 0.62, 0.6, 0.55, 0.5, 0.4]
    n_temps = length(inv_temps)

    spins_configs = zeros(Int, n_mono, n_samples)

    polymers = Array{mp.Magnetic_polymer}(undef, n_temps)
    trajs = Array{mp.MC_data}(undef, n_temps)
    for i_temp in 1:n_temps
        polymers[i_temp] = mp.Magnetic_polymer(n_mono, inv_temps[i_temp], spins_coupling, alpha)
        trajs[i_temp] = mp.MC_data(n_strides*stride)
        if isfile("simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
            mp.initialize_poly!(polymers[i_temp],"simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
        else
            mp.initialize_poly!(polymers[i_temp])
        end
    end

    cov_mat = zeros(n_ss,n_ss)
    #inv_cov_mat = zeros(3,3)
    #delta_w = 0.05
    param_series = zeros(n_feats,n_params)
    SL_series = zeros(n_params)
    ss_mean_series = zeros(n_ss, n_params)
    #ss_cov_determinant_series = zeros(n_params)
    summary_stats = zeros(n_ss, n_samples)
    weights = zeros(n_feats)
    trial_weights = zeros(n_feats)
    weights .= initial_weights

    fields = zeros(n_mono)
    for i in 1:n_feats
        fields .+= features[:,i] .* weights[i]
    end

    mp.set_fields!(polymers, fields)
    for i_eq in 1:20
        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps)
    end

    mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps,spins_configs,sample_lag,n_samples)
    compute_summary_stats!(summary_stats, spins_configs, features)
    ss_mean_series[:,1] .= vec(mean(summary_stats, dims=2))
    cov_mat .= cov(summary_stats, dims=2)

    syn_like = 0.0
    for i_data in 1:n_data
        syn_like += log_synth_likelihood(cov_mat, ss_mean_series[:,1], data[i_data,:])
    end

    param_series[1,1] = weights[1]
    param_series[2,1] = weights[2]
    param_series[3,1] = weights[3]
    SL_series[1] = syn_like
    w_acceptance = 0.0

    for i_param in 2:n_params
        #=trial_weights .= weights
        if rand((1,2))==1
            trial_weights[1] += (2*rand()-1)*delta_w
        else
            trial_weights[2] += (2*rand()-1)*delta_w
        end=#

        trial_weights .= weights .+ (2 .* rand(n_feats) .-1) .*delta_w

        fields = zeros(n_mono)
        for i in 1:n_feats
            fields .+= features[:,i] .* trial_weights[i]
        end
        mp.set_fields!(polymers, fields)

        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps) # This is an short equilibration run so that the expectations are a bit better when changing weights
        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps,spins_configs,sample_lag,n_samples)
        compute_summary_stats!(summary_stats, spins_configs, features)

        ss_mean_series[:,i_param] .= vec(mean(summary_stats, dims=2))
        cov_mat .= cov(summary_stats, dims=2)

        trial_syn_like = 0.0
        for i_data in 1:n_data
            trial_syn_like += log_synth_likelihood(cov_mat, ss_mean_series[:,i_param], data[i_data,:])
        end

        delta_syn_like = trial_syn_like -syn_like
        delta_syn_like>=0 ? w_acceptance=1 : w_acceptance=exp(delta_syn_like)
        if i_param%10 == 0
            println(i_param)
            println("w1: ",weights[1]," ---> ",trial_weights[1])
            println("w2: ",weights[2]," ---> ",trial_weights[2])
            println("w3: ",weights[3]," ---> ",trial_weights[3])
            println("delta_syn_like: ", delta_syn_like)
            println("acceptance: ", w_acceptance)
            println("avg acceptance: ", accepted_moves/i_param)
        end
        if w_acceptance > rand()
            weights .= trial_weights
            syn_like = trial_syn_like
            accepted_moves += 1
        end
        param_series[1,i_param] = weights[1]
        param_series[2,i_param] = weights[2]
        param_series[3,i_param] = weights[3]
        SL_series[i_param] = syn_like
    end
    println("Acceptance ratio: ", accepted_moves/n_params)

    !isdir("SL_data") && mkdir("SL_data")
    open("SL_data/weights.txt", "w") do io
        writedlm(io, param_series)
    end
    open("SL_data/syn_likes.txt", "w") do io
        writedlm(io, SL_series)
    end
    mp.write_results(polymers[1],trajs[1])

end

#####################################################################################
#####################################################################################
#=
# AMHI: Approximate metropolis hastings inference
function amhi_polymer(n_samples::Int, sample_lag::Int, stride::Int, n_params::Int, delta_w::Float64,initial_weights::Array{Float64}, features_file::String, data_file::String)
    # Read features from file
    io = open(features_file, "r")
    features = readdlm(io,Float64)
    close(io)
    n_mono = size(features, 1)
    n_feats = size(features, 2)

    # Read generated data. In this algo they're not summary stats but the full data!
    io = open(data_file,"r")
    data_spins = readdlm(io, Int64; header=true)[1][1:end,:]
    close(io)
    n_data = size(data_spins, 1)
    println(n_mono)
    println(n_data)
    println(typeof(data_spins))

    accepted_moves = 0
    n_strides = cld(n_samples*sample_lag, stride) # integer ceiling of the division
    spins_coupling = 1.0
    #spins_coupling = 0.0
    alpha = 0.5
    inv_temps = [1.0, 0.9, 0.8, 0.74, 0.65, 0.62, 0.6, 0.55, 0.5, 0.4]
    n_temps = length(inv_temps)

    spins_configs = zeros(Int, n_mono, n_samples) # where you store the generated data

    ########################################################### Initialize some quantitites in this section
    polymers = Array{mp.Magnetic_polymer}(undef, n_temps)
    trajs = Array{mp.MC_data}(undef, n_temps)
    for i_temp in 1:n_temps
        polymers[i_temp] = mp.Magnetic_polymer(n_mono, inv_temps[i_temp], spins_coupling, alpha)
        trajs[i_temp] = mp.MC_data(n_strides*stride)
        if isfile("simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
            mp.initialize_poly!(polymers[i_temp],"simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
        else
            mp.initialize_poly!(polymers[i_temp])
        end
    end

    #delta_w = 0.05
    param_series = zeros(n_feats,n_params)
    avg_spins = zeros(n_mono, n_samples)

    weights = zeros(n_feats)
    trial_weights = zeros(n_feats)
    weights .= initial_weights

    fields = zeros(n_mono)
    for i in 1:n_feats
        fields .+= features[:,i] .* weights[i]
    end
    #############################################################

    mp.set_fields!(polymers, fields)
    for i_eq in 1:20
        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps)
    end

    param_series[:,1] .= weights
    w_acceptance = 0.0

    energy = 0
    for i_data in 1:n_data
        for i_mono in 1:n_mono
            energy -= fields[i_mono]*data_spins[i_data,i_mono]
        end
    end
    energy = energy*(1-alpha)

    resample_needed = true
    energy_correction = 0.0

    ##
    w_variation = zeros(n_feats)
    ##

    for i_param in 2:n_params
        #trial_weights .= weights
        #w_component = rand(1:n_feats)
        #w_variation = (2*rand()-1)*delta_w
        #trial_weights[w_component] += w_variation

        w_variation .= (2 .* rand(n_feats) .-1) .*delta_w
        trial_weights .= weights .+ w_variation

        trial_fields = zeros(n_mono)
        for i in 1:n_feats
            trial_fields .+= features[:,i] .* trial_weights[i]
        end

        trial_energy = 0
        for i_data in 1:n_data
            for i_mono in 1:n_mono
                trial_energy -= trial_fields[i_mono]*data_spins[i_data,i_mono]
            end
        end
        trial_energy = trial_energy*(1-alpha)

        # If the previous step was rejected there is no need to estimate energy correction from scratch
        # but maybe not recomputing the estimates will induce some additional bias in the chain
        #if resample_needed
        mp.set_fields!(polymers, fields)

        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps) # This is an short equilibration run so that the expectations are a bit better when changing weights
        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps,spins_configs,sample_lag,n_samples)
        avg_spins .= vec(mean(spins_configs, dims=2))
        energy_correction = 0
        for j in 1:n_feats
            for i_mono in 1:n_mono
                #energy_correction += avg_spins[i_mono]*features[i_mono,w_component]
                energy_correction += avg_spins[i_mono] * features[i_mono, j] * w_variation[j]
            end
        end
        energy_correction = energy_correction * (1-alpha)
        #resample_needed = false
        #end

        #delta_acc = -(trial_energy - energy) - w_variation * energy_correction * n_data
        delta_acc = -(trial_energy - energy) - energy_correction * n_data
        delta_acc>=0 ? w_acceptance=1 : w_acceptance=exp(delta_acc)
        if i_param%10 == 0
            println(i_param)
            println("w1: ",weights[1]," ---> ",trial_weights[1])
            println("w2: ",weights[2]," ---> ",trial_weights[2])
            println("w3: ",weights[3]," ---> ",trial_weights[3])
            println("delta_energy: ", -(trial_energy-energy))
            println("delta_acc: ", delta_acc)
            println("order of error: ", - sum(w_variation.^2)*n_data)
            println("acceptance: ", w_acceptance)
        end
        if w_acceptance > rand()
            weights .= trial_weights
            energy = trial_energy
            resample_needed = true
            fields .= trial_fields
            accepted_moves += 1
        end
        param_series[1,i_param] = weights[1]
        param_series[2,i_param] = weights[2]
        param_series[3,i_param] = weights[3]
    end
    println("Acceptance ratio: ", accepted_moves/n_params)

    !isdir("AMHI_data") && mkdir("AMHI_data")
    open("AMHI_data/weights_$(n_data)amh$(n_samples)_$(delta_w).txt", "w") do io
        writedlm(io, param_series)
    end

    mp.write_results(polymers[1],trajs[1])

end
=#

function amhi_polymer(n_samples::Int, sample_lag::Int, stride::Int, n_params::Int, delta_w::Float64,initial_weights::Array{Float64}, features_file::String, data_file::String)
    # Read features from file
    io = open(features_file, "r")
    features = readdlm(io,Float64)
    close(io)
    n_mono = size(features, 1)
    n_feats = size(features, 2)

    # Read generated data. In this algo they're not summary stats but the full data!
    io = open(data_file,"r")
    data_spins = readdlm(io, Int64; header=true)[1][1:end,:]
    close(io)
    n_data = size(data_spins, 1)
    println(n_mono)
    println(n_data)
    println(typeof(data_spins))

    accepted_moves = 0
    n_strides = cld(n_samples*sample_lag, stride) # integer ceiling of the division
    spins_coupling = 1.0
    #spins_coupling = 0.0 ##############################################
    ####################################################################
    ###################################################################
    ## RESET COUPLING TO 1 AFTER. NOW IT'S JUST TO BENCHMARK WRT A
    ## CASE WHERE WE HAVE ACCESS TO THE REAL POSTERIOR WITH J = 0
    ###################################################################
    ###################################################################
    ###################################################################

    alpha = 0.5
    inv_temps = [1.0, 0.9, 0.8, 0.74, 0.65, 0.62, 0.6, 0.55, 0.5, 0.4]
    n_temps = length(inv_temps)

    spins_configs = zeros(Int, n_mono, n_samples) # where you store the generated data

    ########################################################### Initialize some quantitites in this section
    polymers = Array{mp.Magnetic_polymer}(undef, n_temps)
    trajs = Array{mp.MC_data}(undef, n_temps)
    for i_temp in 1:n_temps
        polymers[i_temp] = mp.Magnetic_polymer(n_mono, inv_temps[i_temp], spins_coupling, alpha)
        trajs[i_temp] = mp.MC_data(n_strides*stride)
        if isfile("simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
            mp.initialize_poly!(polymers[i_temp],"simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
        else
            mp.initialize_poly!(polymers[i_temp])
        end
    end

    param_series = zeros(n_feats,n_params)
    overlaps = zeros(n_feats, n_samples)
    avg_overlaps = zeros(n_feats)
    #cov_overlaps = zeros(n_feats, n_feats)

    weights = zeros(n_feats)
    trial_weights = zeros(n_feats)
    weights .= initial_weights

    fields = zeros(n_mono)
    for i in 1:n_feats
        fields .+= features[:,i] .* weights[i]
    end
    #############################################################

    mp.set_fields!(polymers, fields)
    for i_eq in 1:20
        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps)
    end

    param_series[:,1] .= weights
    w_acceptance = 0.0

    energy = 0
    for i_data in 1:n_data
        for i_mono in 1:n_mono
            energy -= fields[i_mono]*data_spins[i_data,i_mono]
        end
    end
    energy = energy*(1-alpha)

    #resample_needed = true
    energy_correction = 0.0

    ##
    w_variation = zeros(n_feats)
    ##

    for i_param in 2:n_params

        w_variation .= (2 .* rand(n_feats) .-1) .*delta_w
        trial_weights .= weights .+ w_variation

        trial_fields = zeros(n_mono)
        for i in 1:n_feats
            trial_fields .+= features[:,i] .* trial_weights[i]
        end

        trial_energy = 0
        for i_data in 1:n_data
            for i_mono in 1:n_mono
                trial_energy -= trial_fields[i_mono]*data_spins[i_data,i_mono]
            end
        end
        trial_energy = trial_energy*(1-alpha)

        # If the previous step was rejected there is no need to estimate energy correction from scratch
        # but maybe not recomputing the estimates will induce some additional bias in the chain
        #if resample_needed
        mp.set_fields!(polymers, fields)

        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps) # This is an short equilibration run so that the expectations are a bit better when changing weights
        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps,spins_configs,sample_lag,n_samples)

        #=fill!(overlaps, 0.0)
        for i_sample in 1:n_samples
            for i_feat in 1:n_feats
                for i_mono in 1:n_mono
                    overlaps[i_feat, i_sample] += spins_configs[i_mono, i_sample]*features[i_mono, i_feat]
                end
            end
        end=#

        overlaps = features'*spins_configs ## shorter but maybe less efficient way to write stuff in the previous nested for loops
        avg_overlaps = vec(mean(overlaps, dims=2))
        #cov_overlaps = cov(overlaps, dims=2)

        energy_correction = 0
        #linear correction
        for j in 1:n_feats
            energy_correction += w_variation[j]*avg_overlaps[j] * (1-alpha)
        end
        #quadratic correction
        #=quadratic_correction = 0
        for j1 in 1:n_feats
            for j2 in 1:n_feats
                quadratic_correction += w_variation[j1]*cov_overlaps[j1,j2]*w_variation[j2]*(1-alpha)^2 * 0.5
            end
        end
        energy_correction += quadratic_correction=#
        #resample_needed = false
        #end

        delta_acc = -(trial_energy - energy) - energy_correction * n_data
        delta_acc>=0 ? w_acceptance=1 : w_acceptance=exp(delta_acc)
        if i_param%10 == 0
            println(i_param)
            println("w1: ",weights[1]," ---> ",trial_weights[1])
            println("w2: ",weights[2]," ---> ",trial_weights[2])
            println("w3: ",weights[3]," ---> ",trial_weights[3])
            println("delta_energy: ", -(trial_energy-energy))
            println("delta_acc: ", delta_acc)
            println("order of error: ", - sum(w_variation.^2) *n_data)
            println("acceptance: ", w_acceptance)
            println("avg acceptance: ", accepted_moves/i_param)
        end
        if w_acceptance > rand()
            weights .= trial_weights
            energy = trial_energy
            resample_needed = true
            fields .= trial_fields
            accepted_moves += 1
        end
        param_series[1,i_param] = weights[1]
        param_series[2,i_param] = weights[2]
        param_series[3,i_param] = weights[3]
    end
    println("Acceptance ratio: ", accepted_moves/n_params)

    !isdir("AMHI_data") && mkdir("AMHI_data")
    open("AMHI_data/weights_$(n_data)amh$(n_samples)_$(delta_w).txt", "w") do io
        writedlm(io, param_series)
    end

    mp.write_results(polymers[1],trajs[1])

end

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

# QAMHI: Approximate metropolis hastings inference with Quadratic corrections
function Qamhi_polymer(n_samples::Int, sample_lag::Int, stride::Int, n_params::Int, delta_w::Float64,initial_weights::Array{Float64}, features_file::String, data_file::String)
    # Read features from file
    io = open(features_file, "r")
    features = readdlm(io,Float64)
    close(io)
    n_mono = size(features, 1)
    n_feats = size(features, 2)

    # Read generated data. In this algo they're not summary stats but the full data!
    io = open(data_file,"r")
    data_spins = readdlm(io, Int64; header=true)[1][1:end,:]
    close(io)
    n_data = size(data_spins, 1)
    println(n_mono)
    println(n_data)
    println(typeof(data_spins))

    accepted_moves = 0
    n_strides = cld(n_samples*sample_lag, stride) # integer ceiling of the division
    spins_coupling = 1.0
    #spins_coupling = 0.0 ##############################################
    ####################################################################
    ###################################################################
    ## RESET COUPLING TO 1 AFTER. NOW IT'S JUST TO BENCHMARK WRT A
    ## CASE WHERE WE HAVE ACCESS TO THE REAL POSTERIOR WITH J = 0
    ###################################################################
    ###################################################################
    ###################################################################

    alpha = 0.5
    inv_temps = [1.0, 0.9, 0.8, 0.74, 0.65, 0.62, 0.6, 0.55, 0.5, 0.4]
    n_temps = length(inv_temps)

    spins_configs = zeros(Int, n_mono, n_samples) # where you store the generated data

    ########################################################### Initialize some quantitites in this section
    polymers = Array{mp.Magnetic_polymer}(undef, n_temps)
    trajs = Array{mp.MC_data}(undef, n_temps)
    for i_temp in 1:n_temps
        polymers[i_temp] = mp.Magnetic_polymer(n_mono, inv_temps[i_temp], spins_coupling, alpha)
        trajs[i_temp] = mp.MC_data(n_strides*stride)
        if isfile("simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
            mp.initialize_poly!(polymers[i_temp],"simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
        else
            mp.initialize_poly!(polymers[i_temp])
        end
    end

    #delta_w = 0.05
    param_series = zeros(n_feats,n_params)
    #avg_spins = zeros(n_mono, n_samples)
    overlaps = zeros(n_feats, n_samples)
    avg_overlaps = zeros(n_feats)
    cov_overlaps = zeros(n_feats, n_feats)

    weights = zeros(n_feats)
    trial_weights = zeros(n_feats)
    weights .= initial_weights

    fields = zeros(n_mono)
    for i in 1:n_feats
        fields .+= features[:,i] .* weights[i]
    end
    #############################################################

    mp.set_fields!(polymers, fields)
    for i_eq in 1:20
        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps)
    end

    param_series[:,1] .= weights
    w_acceptance = 0.0

    energy = 0
    for i_data in 1:n_data
        for i_mono in 1:n_mono
            energy -= fields[i_mono]*data_spins[i_data,i_mono]
        end
    end
    energy = energy*(1-alpha)

    #resample_needed = true
    energy_correction = 0.0

    ##
    w_variation = zeros(n_feats)
    ##

    for i_param in 2:n_params

        w_variation .= (2 .* rand(n_feats) .-1) .*delta_w
        trial_weights .= weights .+ w_variation

        trial_fields = zeros(n_mono)
        for i in 1:n_feats
            trial_fields .+= features[:,i] .* trial_weights[i]
        end

        trial_energy = 0
        for i_data in 1:n_data
            for i_mono in 1:n_mono
                trial_energy -= trial_fields[i_mono]*data_spins[i_data,i_mono]
            end
        end
        trial_energy = trial_energy*(1-alpha)

        # If the previous step was rejected there is no need to estimate energy correction from scratch
        # but maybe not recomputing the estimates will induce some additional bias in the chain
        #if resample_needed
        mp.set_fields!(polymers, fields)

        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps) # This is an short equilibration run so that the expectations are a bit better when changing weights
        mp.MMC_run!(polymers,trajs,n_strides,stride,inv_temps,spins_configs,sample_lag,n_samples)

        #=fill!(overlaps, 0.0)
        for i_sample in 1:n_samples
            for i_feat in 1:n_feats
                for i_mono in 1:n_mono
                    overlaps[i_feat, i_sample] += spins_configs[i_mono, i_sample]*features[i_mono, i_feat]
                end
            end
        end=#

        overlaps = features'*spins_configs ## shorter but maybe less efficient way to write stuff in the previous nested for loops
        avg_overlaps = vec(mean(overlaps, dims=2))
        cov_overlaps = cov(overlaps, dims=2)

        energy_correction = 0
        #linear correction
        for j in 1:n_feats
            energy_correction += w_variation[j]*avg_overlaps[j] * (1-alpha)
        end
        #quadratic correction
        quadratic_correction = 0
        for j1 in 1:n_feats
            for j2 in 1:n_feats
                quadratic_correction += w_variation[j1]*cov_overlaps[j1,j2]*w_variation[j2]*(1-alpha)^2 * 0.5
            end
        end
        energy_correction += quadratic_correction
        #resample_needed = false
        #end

        delta_acc = -(trial_energy - energy) - energy_correction * n_data
        delta_acc>=0 ? w_acceptance=1 : w_acceptance=exp(delta_acc)
        if i_param%10 == 0
            println(i_param)
            println("w1: ",weights[1]," ---> ",trial_weights[1])
            println("w2: ",weights[2]," ---> ",trial_weights[2])
            println("w3: ",weights[3]," ---> ",trial_weights[3])
            println("delta_energy: ", -(trial_energy-energy))
            println("delta_acc: ", delta_acc)
            println("order of error: ", - sum(w_variation.^2)^1.5 *n_data)
            println("ratio quadratic_correction/tot_corr: ", quadratic_correction/energy_correction)
            println("acceptance: ", w_acceptance)
            println("avg acceptance: ", accepted_moves/i_param)
        end
        if w_acceptance > rand()
            weights .= trial_weights
            energy = trial_energy
            resample_needed = true
            fields .= trial_fields
            accepted_moves += 1
        end
        param_series[1,i_param] = weights[1]
        param_series[2,i_param] = weights[2]
        param_series[3,i_param] = weights[3]
    end
    println("Acceptance ratio: ", accepted_moves/n_params)

    !isdir("AMHI_data") && mkdir("AMHI_data")
    open("AMHI_data/weights_$(n_data)Qamh$(n_samples)_$(delta_w).txt", "w") do io
        writedlm(io, param_series)
    end

    mp.write_results(polymers[1],trajs[1])

end

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

function generate_data(features_file::String, n_strides::Int, weights::Array{Float64})
    io = open(features_file, "r")
    features = readdlm(io,Float64)
    close(io)
    n_mono = size(features, 1)
    n_feats = size(features, 2)

    spins_coupling = 1.0
    alpha = 0.5
    inv_temps = [1.0, 0.9, 0.8, 0.74, 0.65, 0.62, 0.6, 0.55, 0.5, 0.4]
    n_temps = length(inv_temps)
    stride = 500

    polymers = Array{mp.Magnetic_polymer}(undef, n_temps)
    trajs = Array{mp.MC_data}(undef, n_temps)
    for i_temp in 1:n_temps
        polymers[i_temp] = mp.Magnetic_polymer(n_mono, inv_temps[i_temp], spins_coupling, alpha)
        trajs[i_temp] = mp.MC_data(n_strides*stride)
        if isfile("simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
            mp.initialize_poly!(polymers[i_temp],"simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
        else
            mp.initialize_poly!(polymers[i_temp])
        end
    end

    fields = zeros(n_mono)
    for i in 1:n_feats
        fields .+= features[:,i] .* weights[i]
    end
    mp.set_fields!(polymers, fields)

    #here I run the simulation for a while in order to equilibrate the chain
    for i_burnin in 1:50
        println("burnin number: ", i_burnin)
        mp.MMC_run!(polymers, trajs, n_strides, stride, inv_temps)
    end

    n_data = 10
    summary_stats = zeros(n_feats,n_data)
    spins_conf = zeros(Int,n_mono,n_data)

    #
    #poly_confs = zeros(Int,n_data*3,n_mono)
    #

    for i_data in 1:n_data
        println("Number of data point:  ", i_data)
        mp.MMC_run!(polymers, trajs, n_strides, stride, inv_temps)
        #mp.write_results(polymers[1],trajs[1])
        spins_conf[:,i_data] .= polymers[1].spins
        #for j in 1:3
        #    poly_confs[(i_data-1)*3+j,:] .= polymers[1].coord[j,:]
        #end
    end
    #compute_summary_stats!(summary_stats,spins_conf,features)
    #println(summary_stats)
    #println(spins_conf)
    open("saw_conf_data.txt","w") do io
        writedlm(io,transpose(weights))
        writedlm(io,transpose(spins_conf))
    end
    #=open("data_file.txt","w") do io
        writedlm(io,transpose(weights))
        writedlm(io,transpose(summary_stats))
    end=#
    #=open("poly_confs.txt","w") do io
        writedlm(io,transpose(weights))
        writedlm(io, poly_confs)
    end=#

    mp.write_results(polymers[1])
end


text:
# Problem 1 - Subproblem 2

using ControlSystems;

A=[0 0;-1 sqrt(2)];
B=[0 -1;0 sqrt(2)];
Q=[0 0;0 1];
R=[1 0;0 1];

X=care(A,B,Q,R)

max_stars_repo_name: akshayhiregoudar/Control-System-Design

max_issues_repo_name: akshayhiregoudar/Control-System-Design

max_forks_repo_name: akshayhiregoudar/Control-System-Design


text:
using Plots
using Interact
# plotlyjs()
pyplot()

unique_vehicle_ids(high_level) = sort(unique(high_level[:vehicle_ID]))

function plot_number_of_trips_per_vehicle(high_level)
    vehicles = unique_vehicle_ids(high_level)
    counts = map(id->sum(high_level[:vehicle_ID] .== id), vehicles)
    bar(collect(1:length(vehicles)), counts, size=(800,400), xlabel="vehicle", ylabel="number of trips")
end

function plot_histogram_per_vehcile(high_level, target)
    vehicles = unique_vehicle_ids(high_level)
    xlims = extrema(high_level[target])
    p1 = histogram(high_level[find(high_level[:vehicle_ID] .== vehicles[1]), target], title=1, xlims=xlims)
    p2 = histogram(high_level[find(high_level[:vehicle_ID] .== vehicles[2]), target], title=2, xlims=xlims)
    p3 = histogram(high_level[find(high_level[:vehicle_ID] .== vehicles[3]), target], title=3, xlims=xlims)
    p4 = histogram(high_level[find(high_level[:vehicle_ID] .== vehicles[4]), target], title=4, xlims=xlims)
    p5 = histogram(high_level[find(high_level[:vehicle_ID] .== vehicles[5]), target], title=5, xlims=xlims)
    plot(p1, p2, p3, p4, p5, size=(950, 600))
end

function scatter_targets(high_level, A, B)
    vehicles = unique_vehicle_ids(high_level)
    p = plot(xlabel=string(A), ylabel=string(B))
    for (i,id) in enumerate(vehicles)
        x = mean(high_level[find(high_level[:vehicle_ID] .== id), A])
        y = mean(high_level[find(high_level[:vehicle_ID] .== id), B])
        scatter!(p, Float64[x], Float64[y], lab=string(i))
    end
    p
end

function plot_drive(index::Int; scatter::Bool=false)
    utmx = Float64[]
    utmy = Float64[]
    path = anon_paths[index]["utm"]
    for pt in path
        push!(utmx, pt[1])
        push!(utmy, pt[2])
    end
    lox, hix = extrema(utmx)
    loy, hiy = extrema(utmy)
    max_sep = max(hix - lox, hiy - loy)
    midx = (hix + lox)/2
    midy = (hiy + loy)/2
    half_sep = max_sep/2*1.1
    xlims = (midx-half_sep, midx+half_sep)
    ylims = (midy-half_sep, midy+half_sep)
    p = plot(utmx, utmy, xlabel="x position [m]", ylabel="y position [m]", xlims=xlims, ylims=ylims, size=(800,800))

    if scatter == true
        scatter!(p, utmx, utmy)
    end

    p
end


text:
include("StatisticsEstimation/MarkovDiscreteEstimation.jl")
include("StatisticsEstimation/MonteCarloEstimation.jl")
include("StatisticsEstimation/TailEstimation.jl")


text:
using MakieGallery, AbstractPlotting, GLMakie, Test
using MakieGallery: @block, @cell
empty!(MakieGallery.plotting_backends)
push!(MakieGallery.plotting_backends, "GLMakie", "AbstractPlotting")
database = MakieGallery.load_database()

exclude = (
    "Cobweb plot", # has some weird scaling issue on CI
    "Colormap collection", # has one size different...
    # doesn't match 0.035520551315007046 <= 0.032. Looked at the artifacts and it looks fairly similar
    # so blaming video compression
    "Interaction with Mouse"
)
# Download is broken on CI
filter!(entry-> !("download" in entry.tags) && !(entry.title in exclude), database)

tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")
for path in (tested_diff_path, test_record_path)
    rm(path, force = true, recursive = true)
    mkpath(path)
end
recordings = MakieGallery.record_examples(test_record_path)
@test length(recordings) == length(database)
MakieGallery.run_comparison(test_record_path, tested_diff_path)
empty!(database) # remove other examples
include("glmakie_tests.jl") # include GLMakie specific tests
# THese examples download additional data - don't want to deal with that!
for path in (tested_diff_path, test_record_path)
    rm(path, force = true, recursive = true)
    mkpath(path)
end
examples = MakieGallery.record_examples(test_record_path)
MakieGallery.run_comparison(test_record_path, tested_diff_path, maxdiff = 0.00001)
# MakieGallery.generate_preview(test_record_path)
# repo = joinpath(homedir(), "ReferenceImages", "gallery")
# cp(test_record_path, repo, force = true)


text:
#!/usr/bin/env julia
# A very thin wrapper to start a GeMM simulation. If you need something fancier,
# have a look at `rungemmparallel.jl`, or import the GeMM module into your own
# wrapper script.

thisDir = joinpath(pwd(), "src")
any(path -> path == thisDir, LOAD_PATH) || push!(LOAD_PATH, thisDir)
using Pkg
Pkg.activate(".")
using GeMM

rungemm()

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
"""
    RepeatedRangeMatrix(r::AbstractRange{T}, at::AbstractVector{T}) where T

A RepeatedRange is like a RangeMatrix, except that it only stores one range and
a vector of offsets, at which the range repeats. For now, both the range and
vector of offsets must have the same element type.
"""
struct RepeatedRangeMatrix{T,R,A} <: AbstractMatrix{T}
    r::R # <: Range{T}
    at::A #<: AbstractVector{T}
end
RepeatedRangeMatrix(r::AbstractRange{T}, at::AbstractVector{T}) where T =
    RepeatedRangeMatrix{T, typeof(r), typeof(at)}(r, at)

Base.size(R::RepeatedRangeMatrix) = (length(R.r), length(R.at))
Base.IndexStyle(::Type{<:RepeatedRangeMatrix}) = IndexCartesian()

if VERSION < v"0.7.0-DEV.5126"
# This coupled iteration over the two fields is 10-20x faster than Cartesian iteration
# TODO: re-implement in the new iteration protocol
@inline function Base.start(R::RepeatedRangeMatrix)
    is = start(R.r)
    idone = done(R.r, is)
    js = start(R.at)
    jdone = done(R.at, js)
    return (idone | jdone) ? ((one(eltype(R.r)), is), (one(eltype(R.at)), js), true) :
                             (next(R.r, is), next(R.at, js), false)
end
@inline function Base.next(R::RepeatedRangeMatrix, state)
    (i, is), (j, js), _ = state
    val = i + j
    if done(R.r, is)
        if done(R.at, js)
            return (val, ((i, is), (j, js), true))
        end
        is = start(R.r)
        j, js = next(R.at, js)
    end
    return (val, (next(R.r, is), (j, js), false))
end
@inline Base.done(R::RepeatedRangeMatrix, state) = state[end]
end

# Scalar indexing
@inline function Base.getindex(R::RepeatedRangeMatrix, i::Int, j::Int)
    @boundscheck checkbounds(R, i, j)
    @inbounds return R.r[i] + R.at[j]
end

# For non-scalar indexing, only specialize with inner Ranges and Colons to
# return Ranges or RangeMatrixes. For everything else, we can use the fallbacks.
if VERSION < v"0.7.0-alpha"
    @inline function Base.getindex(R::RepeatedRangeMatrix, I::Union{AbstractRange, Colon}, j::Real)
        @boundscheck checkbounds(R, I, j)
        @inbounds return R.r[I] + R.at[j]
    end
else
    @inline function Base.getindex(R::RepeatedRangeMatrix, I::Union{AbstractRange, Colon}, j::Real)
        @boundscheck checkbounds(R, I, j)
        @inbounds return R.r[I] .+ R.at[j]
    end
end
@inline function Base.getindex(R::RepeatedRangeMatrix, I::Union{AbstractRange, Colon}, J)
    @boundscheck checkbounds(R, I, J)
    @inbounds return RepeatedRangeMatrix(R.r[I], R.at[J])
end

# We can also optimize bounds checks to only look at the range's endpoints
function Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, R::RepeatedRangeMatrix)
    b = true
    @inbounds for a in R.at
        b &= checkindex(Bool, inds, R.r[1] + a)
        b &= checkindex(Bool, inds, R.r[end] + a)
    end
    b
end


text:
export PPOTrajectory, PPOActionMaskTrajectory

using MacroTools

const PPOTrajectory = CombinedTrajectory{
    <:SharedTrajectory{
        <:CircularArrayBuffer,
        <:NamedTuple{(:action_log_prob, :next_action_log_prob, :full_action_log_prob)},
    },
    <:CircularCompactSARTSATrajectory,
}

function PPOTrajectory(;
    capacity,
    action_log_prob_size = (),
    action_log_prob_type = Float32,
    kw...,
)
    CombinedTrajectory(
        SharedTrajectory(
            CircularArrayBuffer{action_log_prob_type}(
                action_log_prob_size...,
                capacity + 1,
            ),
            :action_log_prob,
        ),
        CircularCompactSARTSATrajectory(; capacity = capacity, kw...),
    )
end

const PPOActionMaskTrajectory = CombinedTrajectory{
    <:SharedTrajectory{
        <:CircularArrayBuffer,
        <:NamedTuple{(:action_log_prob, :next_action_log_prob, :full_action_log_prob)},
    },
    <:CircularCompactSALRTSALTrajectory,
}

function PPOActionMaskTrajectory(;
    capacity,
    action_log_prob_size = (),
    action_log_prob_type = Float32,
    kw...,
)
    CombinedTrajectory(
        SharedTrajectory(
            CircularArrayBuffer{action_log_prob_type}(
                action_log_prob_size...,
                capacity + 1,
            ),
            :action_log_prob,
        ),
        CircularCompactSALRTSALTrajectory(; capacity = capacity, kw...),
    )
end