text:
using DSGE, Test

writing_output = false

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:forecast_horizons, 12)

# Define dummy scenario
scen = Scenario(:testscen, "Test Scenario", [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF")

# Check compute_scenario_system throws error before zeroing out measurement error
@testset "Check compute_scenario_system error throwing" begin
    @test_throws ErrorException compute_scenario_system(m, scen)
end

# Zero out measurement error (or else filtering shocks won't work exactly)
for para in [:e_y, :e_π, :e_R]
    m[para].value = 0
end

# Test compute_scenario_system
@testset "Test compute_scenario_system, filter_shocks!, and shock_scaling" begin
    global sys = compute_scenario_system(m, scen)
    @test all(x -> x == 0, sys[:CCC])
    @test all(x -> x == 0, sys[:DD])
    @test all(x -> x == 0, sys[:DD_pseudo])
    for shock in keys(m.exogenous_shocks)
        ind = m.exogenous_shocks[shock]
        if shock in scen.instrument_names
            @test sys[:QQ][ind, ind] != 0
        else
            @test sys[:QQ][ind, ind] == 0
        end
    end

    # Test filter_shocks!
    for var in scen.target_names
        scen.targets[!, var] = rand(forecast_horizons(m))
    end
    global forecastshocks = filter_shocks!(m, scen, sys)
    for var in scen.instrument_names
        ind = m.exogenous_shocks[var]
        @test scen.instruments[!, var] == forecastshocks[ind, :]
    end

    s_T = zeros(n_states_augmented(m))
    global _, forecastobs,_ ,_  = forecast(m, sys, s_T, shocks = forecastshocks)
    for var in scen.target_names
        ind = m.observables[var]
        @test scen.targets[!, var] ≈ forecastobs[ind, :]
    end

    # Test scenario with shock scaling
    global scale_scenario = Scenario(:scaledscen, "Test Shock Scaling Scenario",
                            [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF",
                            shock_scaling = 2.0)
    scale_scenario.targets = scen.targets

    global forecastshocks = filter_shocks!(m, scale_scenario, sys)
    for var in scale_scenario.target_names
        ind = m.observables[var]
        @test scale_scenario.targets[!, var] == scen.targets[!, var]
    end
    for var in scale_scenario.instrument_names
        ind = m.exogenous_shocks[var]
        @test scale_scenario.instruments[!, var] == forecastshocks[ind, :]
    end
end

###########################
# forecast_scenario
###########################
#=
file = "$path/../reference/forecast_scenario_draw.jld2"
sys = compute_scenario_system(m, scen)

forecast_output = DSGE.forecast_scenario_draw(m, scen, sys, 1)

if writing_output
    JLD2.jldopen(file, true, true, true, IOStream) do file
        file["forecast_output"] = forecast_output
    end
end
forecast_output_test = load(file, "forecast_output")

@testset "Test forecast_scenario_draw, forecast_scenario" begin
    @test forecast_output == forecast_output_test
end
=#


text:
using Dates
using PyPlot
using Plots
#using ProfileView

include("../btg.jl")
include("../datasets/load_abalone.jl")

#ind = 350:370
#ind = 3300:3359
ind = 3300:3500
#posx = 1:3 #
#posx = [1;4]
posx = 1:7
#posc = 1:1
posc = 1:1
x = data[ind, posx]
#choose a subset of variables to be regressors for the mean
Fx = data[ind, posc]
y = float(target[ind])/29
pind = 10:10 #prediction index
#pind = i:i
trainingData1 = trainingData(x, Fx, y) #training data used for testing various functions
#prediction data
d = getDimension(trainingData1); n = getNumPts(trainingData1); p = getCovDimension(trainingData1)
Fx0 = reshape(data[pind, posc], 1, length(posc))
x0 = reshape(data[pind, posx], 1, length(posx))
#rangeθ = [100.0 200]
#rangeθ = [100.0 200]
#rangeθ = select_single_theta_range(x)
rangeθ = [500.0 1000]
rangeλ = [-2.0 2.0] #we will always used 1 range scale for lambda
btg1 = btg(trainingData1, rangeθ, rangeλ; quadtype = ["Gaussian", "Gaussian"])
#θ1 = btg1.nodesWeightsθ.nodes[1, 6] #pick some theta value which doubles as quadrature node
#λ1 = btg1.nodesWeightsλ.nodes[3]
##################################################################################################
PyPlot.close("all") #close existing windows

#(pdf, cdf, dpdf) = solve(btg1); #initialize training_buffer_dicts, solve once so can use fast techiques to extrapolate submatrix determinants, etc.
#pdf(x0, Fx0, .4);
#(pdf, cdf, dpdf) = solve(btg1; validate = i) #do LOOCV this time around
#a1 = y-> pdf(x0, Fx0, y)

println("recommended theta range:", select_single_theta_range(x))

"""
leave-one-out-training-data
"""
function lootd(td::AbstractTrainingData, i::Int64)
    x = getPosition(td)
    Fx = getCovariates(td)
    z = getLabel(td)
    x_minus_i = x[[1:i-1;i+1:end], :]
    Fx_minus_i = Fx[[1:i-1;i+1:end], :]
    z_minus_i = z[[1:i-1;i+1:end]]
    x_i = x[i:i, :]
    Fx_i = Fx[i:i, :]
    z_i = z[i:i, :]
    return trainingData(x_minus_i, Fx_minus_i, z_minus_i), x_i, Fx_i, z_i
end
"""
Run
"""
function test_LOOCV(btg1::btg, m::Int64, n::Int64; fast)
    fast = fast
    m = m; n=n
    #plt, axs = PyPlot.subplots(m, n,figsize=[20,11])
    plt, axs = PyPlot.subplots(m, n)
    #figure(1)
    before = Dates.now()
    println("LOOCV type is fast?: $fast")
    for j = 1:m*n
        ind1 = Int64(ceil(j/n))
        ind2 = Int64(j - n*(floor((j-.1)/n)))
        println("iteration $j")
        if j <m*n
            if fast
                (pdf, cdf, dpdf) = solve(btg1, validate = j);
                b1 = y0 -> pdf(x0, Fx0, y0);
                c1 = y0 -> cdf(x0, Fx0, y0);
                (x, y) = plt_data(b1, .01, 1.2, 100)
                (xc, yc) = plt_data(c1, .01, 1.2, 100)
            else
                (trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, j);
                #btg2 = btg(trainingdata_minus_i, rangeθ, rangeλ; quadtype = ["MonteCarlo", "MonteCarlo"])
                btg2 = btg(trainingdata_minus_i, rangeθ, rangeλ; quadtype = ["Gaussian", "Gaussian"]);
                (pdf_minus_i, cdf_minus_i, dpdf_minus_i) = solve(btg2);
                b1 = y -> pdf_minus_i(x_i, Fx_i, y);
                c1 = y -> cdf_minus_i(x_i, Fx_i, y);
                (x, y) = plt_data(b1, .01, 1.2, 100);
                (xc, yc) = plt_data(c1, .01, 1.2, 100);
            end
            axs[ind1, ind2].plot(x, y, color = "red", linewidth = 2.0, linestyle = ":")
            axs[ind1, ind2].plot(xc, yc, color = "orange", linewidth = 2.0, linestyle = ":")
            #axs[ind1, ind2].plot(x1, y1, color = "blue", linewidth = 1.0, linestyle = "--")
            axs[ind1, ind2].axvline(x =  getLabel(btg1.trainingData)[j])
        else #do annotation in last box
            #println("annotating...")
            after = Dates.now()
            elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)
            qtype =  btg1.quadType
            axs[ind1, ind2].annotate(
            L"$\theta: $" * "$rangeθ" * "\n" *
            L"$\lambda: $" * "$rangeλ" * "\n" *
            "$qtype \n" *
            "ind: " * "$ind \n" *
            "x: " * "$posx " *", cov: " * "$posc\n" *
            "time: " * "$elapsedmin",
            xy=[.03;.03],
            fontsize = 5.0)
            # *
            #L"$\lambda: $" * "$rangeλ" * ".\n" *
            #"quad: " * "$btg1.quadType",
        end
    end
    #for ax in axs
    #    ax.set(xlabel="x-label", ylabel="y-label")
    #end
    # Hide x labels and tick labels for top plots and y ticks for right plots.
    for ax in axs
        ax.label_outer()
    end
end

if true
    test_LOOCV(btg1, 4, 5; fast = false);
end

###
### reserved REPL testing block
###

if true
i=3
#function isolate(i)
    (trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, i);
    btg2 = btg(trainingdata_minus_i, rangeθ, rangeλ; quadtype = ["Gaussian", "Gaussian"]);
    (pdf, cdf, dpdf, augmented_cdf_deriv, augmented_cdf_hess, quantInfo, tgridpdf, tgridcdf, tgridm, tgridsigma_m, weightsTensorGrid)  = solve(btg2);

    function func_grid_eval(funcgrid, x_i, Fx_i, y)
        g = similar(weightsTensorGrid)
        R = CartesianIndices(weightsTensorGrid)
        for I in R
            g[I] = funcgrid[I](x_i, Fx_i, y)
        end
        return g
    end

    b1 = y -> pdf_minus_i(x_i, Fx_i, y);
    c1 = y -> cdf_minus_i(x_i, Fx_i, y);
    d1 = y->  func_grid_eval(tgridpdf, x_i, Fx_i, y);
    e1 = y-> func_grid_eval(tgridcdf, x_i, Fx_i, y);
    f1 = () -> func_grid_eval(tgridm, x_i, Fx_i, y);
    g1 = () -> func_grid_eval(tgridsigma, x_i, Fx_i, y);
    (x, y) = plt_data(b1, .01, 1.2, 100);
    (xc, yc) = plt_data(c1, .01, 1.2, 100);
    #return (b1, c1, d1, e1, f1, g1, weightsTensorGrid)
#end

if false
    (trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData1, 3);
    (b1, c1, d1, e1, f1, g1, w) = isolate(3);
end

#manual computation for checking

#btg2.θλbuffer.θ
#btg2.θλbuffer.λ
#btg2.θλbuffer.Σθ_inv_y

# testing data
x0_cur = btg2.testingData.x0;
Fx0_cur = btg2.testingData.Fx0;
#fix theta and lambda - run thru entire computation to get m, and see if it's negative...
theta_cur = btg2.nodesWeightsθ.nodes[1];
lambda_cur = btg2.nodesWeightsλ.nodes[1];
#computation by hand
Σθ_cur = correlation(Gaussian(), theta_cur, trainingdata_minus_i.x);
Σθ_inv_X_cur = Σθ_cur\trainingdata_minus_i.Fx;
Bθ_cur = cross_correlation(Gaussian(), theta_cur, x0_cur, trainingdata_minus_i.x);
Fx_cur = trainingdata_minus_i.Fx;
glz_cur = copy((trainingdata_minus_i.y));
for i = 1:length(glz_cur)
    glz_cur[i] = btg2.g(glz_cur[i], lambda_cur);
end
beta_hat_cur = (Fx_cur' * (Σθ_cur\Fx_cur))\(Fx_cur'*(Σθ_cur\glz_cur));
H_theta_cur = Fx_i - Bθ_cur*(Σθ_cur\Fx_cur);
m_cur = Bθ_cur*(Σθ_cur\glz_cur) + H_theta_cur*beta_hat_cur

#how sigma theta is computed in code
m=200;
Σθ = Array{Float64}(undef, 1000, 1000);
Σθ[1:m, 1:m] = correlation(Gaussian(), theta_cur, trainingdata_minus_i.x[1:m, :]; jitter = 1e-8);
choleskyΣθ = incremental_cholesky!(Σθ, m);

#stuff computed using btg code
Σθ_inv_X = btg2.train_buffer_dict[theta_cur].Σθ_inv_X;
cholΣθ = btg2.train_buffer_dict[theta_cur].choleskyΣθ;
Bθ = btg2.test_buffer_dict[theta_cur].Bθ;
gλz = btg2.λbuffer_dict[lambda_cur].gλz;
βhat = btg2.θλbuffer_dict[(theta_cur, lambda_cur)].βhat;
Hθ =  btg2.test_buffer_dict[theta_cur].Hθ;

id = Diagonal(ones(trainingData1.n-1));
rhs = rand(trainingData1.n-1, 1);

@assert norm(trainingdata_minus_i.x - btg2.trainingData.x) <0.0001 #training data matches up
@info "Σθ_inv_X_cur - Σθ_inv_X", norm(Σθ_inv_X_cur - Σθ_inv_X);
@info "rand Σθ_cur norm diff", norm(Σθ_cur\rhs - cholΣθ\rhs);
@info "Bθ_cur - Bθ", norm(Bθ - Bθ_cur);
#@info "rand Σθ_cur norm diff", norm(cholΣθ\rhs - choleskyΣθ\rhs);
@info "sigma theta inv B theta", norm(Σθ_cur\Bθ' - choleskyΣθ\Bθ');
@info "sigma theta inv B same matrix", ( ccc = cholesky(Σθ_cur);  norm(Σθ_cur\Bθ' - ccc\Bθ'));
@info "glz nrom diff", norm(glz_cur - gλz);
@info "betahat norm diff", norm(βhat - beta_hat_cur);
#@info btg2.nodesWeightsθ.nodes
#@info btg2.nodesWeightsλ.nodes
#@info btg2.train_buffer_dict[btg2.nodesWeightsθ.nodes[1]].choleskyΣθ
end

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# mpiexec.exe -n 4 julia .\mpisum.jl
using ADCME

mpi_init()
r = mpi_rank()
a = constant(ones(10) * r)
b = mpi_bcast(a, 3)
L = sum(b^2)
L = mpi_sum(L)
g = gradients(L, a)

sess = Session(); init(sess)
v, G = run(sess, [b, G])

@info r, v

# gradient test

# expected [ Info: (3, [24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0])
G = run(sess, g)

@info r, G

mpi_finalize()


text:
function find_pos_val(val, input_string)
    t= split(split(input_string)[8],";")
    ind = [occursin(val,i) for i in t]
    y= replace(t[ind][1],val => "")
    return y
end

"""
    get_sv_attrs(in_str::String)::Dict

Parse a line from the vcf file and return a dictionary with the following attributes as keys
    - start: SV start position
    - end: SV end position
    - chr: SV Chromosome
    - svlen: Length tof the SV
    - svtype: The SV type

"""
function get_sv_attrs(in_str)
    dict = Dict()
    dict["chr"] = split(in_str)[1]
    dict["start"] = split(in_str)[2]
    dict["end"] = find_pos_val("END=",in_str)
    dict["svlen"] = find_pos_val("SVLEN=",in_str)
    dict["svtype"] = find_pos_val("SVTYPE=",in_str)
    dict["seq"] = split(in_str)[5]
    return dict
end

function annotate_sv(sv,x, y)
    if sv =="Unique" || sv == "unclassified"
        return x*"\t"*sv*"\tNA\tNA"

    else
        return x*"\t"*sv*"\t"*replace(y,"\t" => "|")
    end
end

"""
    compare_locations(x::String, y:String)::String

Compare two lines each from a single vcf file and return what the status of the location is. The function checks for the following
  - Unique: start and end not equal
  - Merged: Start and end equal and of same SV type
  - singleton: Either one of start or end matches and of the same SV type
  - unclassified: All other cases
"""
function compare_locations(x,y)
    x_attrs = get_sv_attrs(x)
    y_attrs = get_sv_attrs(y)
    if x_attrs["chr"] != y_attrs["chr"]
        return "Unique"
    end
    #println(x)
    consensus = x_attrs["seq"] == y_attrs["seq"]
    if x_attrs["start"] != y_attrs["start"] &&  x_attrs["end"] != y_attrs["end"]
        return "Unique"

    elseif x_attrs["start"]== y_attrs["start"] && x_attrs["end"] == y_attrs["end"] && x_attrs["svtype"] == y_attrs["svtype"]
        return "Both"*"\t"*string(consensus)

    elseif  x_attrs["start"]== y_attrs["start"] && x_attrs["end"] != y_attrs["end"] && x_attrs["svtype"] == y_attrs["svtype"]
        return "Left"*"\t"*string(consensus)

    elseif  x_attrs["start"]!= y_attrs["start"] && x_attrs["end"] == y_attrs["end"] && x_attrs["svtype"] == y_attrs["svtype"]
        return "Right"*"\t"*string(consensus)

    else
        return "unclassified"
    end
end

"""
    check_pos(x:String, y::String)::Bool

Compare two input SV lines to make sure that the position of the first line is lower than or equal to the other
"""
function check_pos(x,y)
    chr_start_x = split(x)[1]
    pos_x = split(x)[2]
    chr_start_y = split(y[1])[1]
    pos_y = split(y[1])[2]
    if parse(Int,pos_x) <= parse(Int,pos_y) && chr_start_x == chr_start_y
        return true
    elseif parse(Int,pos_x ) >  parse(Int,pos_y) && chr_start_x != chr_start_y
        return true
    else
        return false
    end
end

function annotate_pindel(v1,v2)
    check_line = iterate(v2)
    v3=[]
    i = 1
    j = 1
    for line1 in v1
        if check_pos(line1,check_line)
            sv = compare_locations(line1,check_line[1])
        elseif !check_pos(line1,check_line) && iterate(v2, check_line[2]) != nothing
            while !check_pos(line1,check_line)
                check_line = iterate(v2,check_line[2])
            end
            sv= compare_locations(line1,check_line[1])
        else
            sv= compare_locations(line1,check_line[1])
        end
        val = annotate_sv(sv,line1,check_line[1])
        #println(val)
        push!(v3,val)
        #println(string(i)*" " *string(j))
        i +=1
    end
    return v3
end

function print_res(filePath,vec)
    fheader = "CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tGT\tAD_REF\tAD_ALT"
    fheader = fheader*"\tStatus_of_merge\tSeq_overlap_status\tResults_of_merge"
    fheader = fheader*"\tSV_Type\tSV_length\tHompolymer_length\tfrac_ALT"
    f = open(filePath,"w")
    println(f,fheader)
    for line in vec
        tmpline = split(line)
        tmpline[10] =replace(tmpline[10],":" => "\t")
        tmpline[10] =replace(tmpline[10],"," => "\t")
        AD_ref = parse(Float64,split(tmpline[10])[2])
        AD_alt = parse(Float64,split(tmpline[10])[3])

        tmpline_to_print = tmpline
        tmpline_to_print[4] = first(tmpline_to_print[4],10)
        tmpline_to_print[5] = first(tmpline_to_print[5],10)
        tmpline_to_print[13] = first(tmpline_to_print[13],10)

        push!(tmpline_to_print, find_pos_val("SVTYPE=",line))
        push!(tmpline_to_print, find_pos_val("SVLEN=",line))
        push!(tmpline_to_print, find_pos_val("HOMLEN=",line))
        push!(tmpline_to_print, string(AD_alt/(AD_alt+AD_ref)))
        println(f,join(tmpline_to_print,"\t"))
    end
    close(f)
end

function run_pindel_compare(samp1, samp2; datdir="/Users/aragaven/scratch/test_julia/test_pindel/pindel_vcf",
                            outdir= "/Users/aragaven/Documents/Research/CBC/ene/")
    #datdir = "/Users/aragaven/scratch/test_julia/test_pindel/pindel_vcf";
    #samp1 = "DSY303_OP";
    #samp2 = "DSY303_WH";

    vec1 = readlines(joinpath(datdir,samp1*".vcf"));
    vec2 = readlines(joinpath(datdir,samp2*".vcf"));

    vec1_1 = filter(x -> !startswith(x, "#"), vec1);
    vec2_2 = filter(x -> !startswith(x, "#"), vec2);

    vec4 = annotate_pindel(vec1_1,vec2_2);
    vec5 = annotate_pindel(vec2_2,vec1_1);

    filep1 = outdir*"/"*samp1*".results.tsv"
    print_res(filep1,vec4)
    filep2 = outdir*"/"*samp2*".results.tsv"
    print_res(filep2,vec5)
end

#run_pindel_compare("DSY303_OP","DSY303_WH")
run_pindel_compare("RBY1118_OP","RBY1118_WH")

max_stars_repo_name: compbiocore/BioJulia-cookbook

max_issues_repo_name: compbiocore/BioJulia-cookbook

max_forks_repo_name: compbiocore/BioJulia-cookbook


text:
function make_S_N_NG86_table(code::GeneticCode)
    f = i -> S_N_NG86(i, code)
    return S_N_NG86_LOOKUP(f)
end

function make_DS_DN_NG86_table(code::GeneticCode)
    f = (i, j) -> DS_DN_NG86(i, j, code)
    return DS_DN_NG86_LOOKUP(f)
end

isinteractive() && info("Compiling lookup tables for NG86...")
const S_N_NG86_LOOKUPS = (function ()
    lookups = Dict{Int, S_N_NG86_LOOKUP}()
    for (i, code) in ncbi_trans_table.tables
        lookups[i] = make_S_N_NG86_table(code)
    end
    return lookups
end)()

const DS_DN_NG86_LOOKUPS = (function ()
    lookups = Dict{Int, DS_DN_NG86_LOOKUP}()
    for (i, code) in ncbi_trans_table.tables
        lookups[i] = make_DS_DN_NG86_table(code)
    end
    return lookups
end)()


text:
struct SubSetSimulation
    n::Integer
    target::Float64
    levels::Integer
    proposal::Sampleable{Univariate}

    function SubSetSimulation(
        n::Integer, target::Float64, levels::Integer, proposal::Sampleable{Univariate}
    )
        skewness(proposal) != 0.0 && error("proposal must be a symmetric distribution")
        mean(proposal) != median(proposal) &&
            error("proposal must be a symmetric distribution")
        return new(n, target, levels, proposal)
    end
end

function sample(inputs::Array{<:UQInput}, sim::SubSetSimulation)
    random_inputs = filter(i -> isa(i, RandomUQInput), inputs)
    deterministic_inputs = filter(i -> isa(i, DeterministicUQInput), inputs)

    n_rv = count_rvs(random_inputs)
    rv_names = names(random_inputs)

    samples = DataFrame(rv_names .=> eachcol(rand(Normal(), sim.n, n_rv)))

    to_physical_space!(random_inputs, samples)

    if !isempty(deterministic_inputs)
        samples = hcat(samples, sample(deterministic_inputs, sim.n))
    end

    return samples
end

function probability_of_failure(
    models::Union{Array{<:UQModel},UQModel},
    performancefunction::Function,
    inputs::Union{Array{T},T} where {T<:UQInput},
    sim::SubSetSimulation,
)
    random_inputs = filter(i -> isa(i, RandomUQInput), inputs)
    rvs = names(random_inputs)

    samples = [sample(inputs, sim)]

    evaluate!(models, samples[end])

    performance = [performancefunction(samples[end])]

    number_of_chains = Int64(max(1, ceil(sim.n * sim.target)))
    samples_per_chain = floor(sim.n / number_of_chains)

    threshold = zeros(sim.levels, 1)
    pf = ones(sim.levels, 1)

    for i in 1:(sim.levels)
        sorted_performance = sort(performance[end])
        sorted_indices = sortperm(performance[end])

        threshold[i] = sorted_performance[number_of_chains]
        pf[i] = if threshold[i] <= 0
            mean(performance[end] .<= 0)
        else
            mean(performance[end] .<= threshold[i])
        end

        if threshold[i] <= 0 || i == sim.levels
            break
        end

        nextlevelsamples = [samples[end][sorted_indices[1:number_of_chains], :]]
        nextlevelperformance = [sorted_performance[1:number_of_chains]]

        # Modified metropolis hastings to generate samples for the next intermediate failure region
        for c in 1:samples_per_chain
            chainsamples = copy(nextlevelsamples[end])

            to_standard_normal_space!(inputs, chainsamples)

            chainsamples[:, rvs] = candidatesamples(
                Matrix{Float64}(chainsamples[:, rvs]), sim.proposal
            )

            to_physical_space!(inputs, chainsamples)

            evaluate!(models, chainsamples)

            chainperformance = performancefunction(chainsamples)

            reject = chainperformance .> threshold[i]

            chainperformance[reject] = nextlevelperformance[end][reject]
            chainsamples[reject, rvs] = nextlevelsamples[end][reject, rvs]

            if c == 1
                nextlevelsamples = [chainsamples]
                nextlevelperformance = [chainperformance]
            else
                push!(nextlevelsamples, chainsamples)
                push!(nextlevelperformance, chainperformance)
            end
        end

        nextlevelsamples = reduce(vcat, nextlevelsamples)
        nextlevelperformance = reduce(vcat, nextlevelperformance)

        push!(samples, nextlevelsamples)
        push!(performance, nextlevelperformance)
    end

    # add level to each dataframe
    for i in eachindex(samples)
        samples[i][!, :level] .= i
    end
    # merge (vcat) all samples
    samples = reduce(vcat, samples)

    pf = prod(pf)

    return pf, samples
end

function candidatesamples(θ::AbstractMatrix, proposal::Sampleable{Univariate})
    Φ = MvNormal(size(θ, 2), 1.0)

    ξ = θ + rand(proposal, size(θ)...)
    α = pdf(Φ, transpose(ξ)) ./ pdf(Φ, transpose(θ))

    accept = α .>= rand(size(α)...)
    θ[accept, :] = ξ[accept, :]

    return θ
end


text:
function format_number_test()
    @testset "Format number" begin
        y = format_number(123456.789);
        @test startswith(y,  "123,456")
        y = format_number(123456.789; format = :none);
        @test startswith(y, "123456.789")

        y = format_number(123.456);
        @test startswith(y, "123.4")

        y = format_number(0.0012345);
        @test startswith(y, "0.00123")
    end
end

function format_vector_test()
    s = format_vector("abc", [1.2, 2.3, 34.56789], 2);
    @test startswith(s, "abc:")
    @test endswith(s, "34.57]")
end

@testset "Display" begin
    format_number_test();
    format_vector_test();
end

# -------------


text:
function falseposition(f::Function, xl::Number, xu::Number, ea::Number, maxIter::Int64)
  # Estimate roots of functions using the false position method
  xr = NaN; # Initial estimate of the root
  fxu = f(xu);
  fxl = f(xl);

  for i = 1:maxIter
    xrold = xr;
    xr = xu - fxu*(xl-xu) / (fxl-fxu);
    fxr = f(xr);
    relativeError = abs((xr - xrold) / xr) * 100;

    if fxr * fxl == 0
      break;
    elseif fxr * fxl > 0
      xl = xr;
      fxl = fxr;
    else
      xu = xr;
      fxu = fxr;
    end

    if relativeError < ea
      break;
    end
  end
  return xr;
end

max_stars_repo_path: Julia/Root Estimation/falseposition.jl

max_stars_repo_name: alexjohnj/numerical-methods

max_issues_repo_path: Julia/Root Estimation/falseposition.jl

max_issues_repo_name: alexjohnj/numerical-methods

max_forks_repo_path: Julia/Root Estimation/falseposition.jl

max_forks_repo_name: alexjohnj/numerical-methods


text:
##############################
#
# obj-Files
#
##############################

FileIO.load(fn::File{format"OBJ"}, MeshType=GLNormalMesh) = open(fn) do s
    load(s, MeshType)
end

function FileIO.load{MT <: Mesh}(io::Stream{format"OBJ"}, MeshType::Type{MT}=GLNormalMesh)
    io           = stream(io)
    lineNumber   = 1
    mesh         = MeshType()
    last_command = ""
    attrib_type  = nothing
    for line in eachline(io)
        # read a line, remove newline and leading/trailing whitespaces
        line = strip(chomp(line))
        !isvalid(line) && error("non valid ascii in obj")

        if !startswith(line, "#") && !isempty(line) && !iscntrl(line) #ignore comments
            lines   = split(line)
            command = shift!(lines) #first is the command, rest the data

            if "v" == command # mesh always has vertices
                push!(vertices(mesh), Point{3, Float32}(lines)) # should automatically convert to the right type in vertices(mesh)
            elseif "vn" == command && has_normals(mesh)
                push!(normals(mesh), Normal{3, Float32}(lines))
            elseif "vt" == command && has_texturecoordinates(mesh)
                length(lines) == 2 && push!(lines, "0.0") # length can be two, but obj normaly does three coordinates with the third equals 0.
                push!(texturecoordinates(mesh), UVW{Float32}(lines))
            elseif "f" == command #mesh always has faces
                if any(x->contains(x, "/"), lines)
                    fs = process_face_uv_or_normal(lines)
                elseif any(x->contains(x, "//"), lines)
                    fs = process_face_normal(lines)
                else
                    return push!(faces(mesh), Triangle{Uint32}(lines))
                end
                push!(faces(mesh), Triangle{Uint32}(map(first, fs)))
            else
                #TODO
            end
        end
        # read next line
        lineNumber += 1
    end
    if isempty(mesh.normals) || length(mesh.normals) != length(mesh.vertices) # silly way of dealing with normals that are not per vertex.
        empty!(mesh.normals)
        append!(mesh.normals, normals(mesh.vertices, mesh.faces))
    end
    return mesh
end

# face indices are allowed to be negative, this methods handles this correctly
function handle_index{T <: Integer}(bufferlength::Integer, s::AbstractString, index_type::Type{T})
    i = parse(T, s)
    i < 0 && return convert(T, bufferlength) + i + one(T) # account for negative indexes
    return i
end
function push_index!{T}(buffer::Vector{T}, s::AbstractString)
    push!(buffer, handle_index(length(buffer), s, eltype(T)))
end

immutable SplitFunctor <: Base.Func{1}
    seperator::ASCIIString
end
call(s::SplitFunctor, array) = split(array, s.seperator)

# of form "f v1 v2 v3 ....""
process_face{S <: AbstractString}(lines::Vector{S}) = (lines,) # just put it in the same format as the others
# of form "f v1//vn1 v2//vn2 v3//vn3 ..."
process_face_normal{S <: AbstractString}(lines::Vector{S}) = map(SplitFunctor("//"), lines)
# of form "f v1/vt1 v2/vt2 v3/vt3 ..." or of form "f v1/vt1/vn1 v2/vt2/vn2 v3/vt3/vn3 ...."
process_face_uv_or_normal{S <: AbstractString}(lines::Vector{S}) = map(SplitFunctor("/"), lines)

function process{F <: Face, S <: AbstractString}(::Type{F}, lines::Vector{S}, mesh::Mesh, line::Int)
    if any(x->contains(x, "/"), lines)
        faces = process_face_uv_or_normal(lines)
    elseif any(x->contains(x, "//"), lines)
        faces = process_face_normal(lines)
    else
        return push!(mesh[F], F(Triangle{Uint32}(lines)))
    end
    push!(mesh[F], F(Triangle{Uint32}(map(first, faces))))
end


text:
using Random, Test, CommonLH

function compare_test()
    @testset "Comparisons" begin
        rng = MersenneTwister(12);
        x = rand(rng, 4,3);
        xMin, xMax = extrema(x);
        @test all_at_least(x, xMin)
        @test !all_at_least(x, xMin + 1e-6)
        @test all_at_most(x, xMax)
        @test !all_at_most(x, xMax - 1e-6)

        @test all_at_least(x, x .- 1e-5)
        @test all_at_most(x, x .+ 1e-5)
    end
end

function scale_array_test()
    @testset "Scale array" begin
        x = rand(4,3,2);
        for d = 1 : 3
            totalV = collect(range(1.0, 2.0, length = size(x, d)));
            scale_array!(x, d, totalV);
            for j = 1 : size(x, d)
                @test isapprox(sum(selectdim(x, d, j)),  totalV[j])
            end
        end
	end
end

function bracket_array_test()
    @testset "Bracket array" begin
        x = [1.0 1.1 1.2; 2.0 2.1 2.2];
        ub = x .+ 0.01;
        lb = x .- 0.01;
        z = copy(x);
        bracket_array!(z, lb, ub);
        @test isequal(z, x)

        ub[3] = x[3] - 0.001;
        lb[4] = x[4] + 0.001;
        z = copy(x);
        bracket_array!(z, lb, ub);
        @test !isequal(z, x)
        @test all(z .<= ub)
        @test all(z .>= lb)
	end
end

function bracket_array_scalar_test()
    @testset "Bracket array in scalars" begin
        x = [1.0 1.1 1.2; 2.0 2.1 2.2];
        ub = 2.099;
        lb = 1.099;
        bracket_array!(x, lb, ub);
        @test all(x .<= ub);
        @test all(x .>= lb);
    end
end

@testset "Arrays" begin
    compare_test();
    scale_array_test()
    bracket_array_test();
    bracket_array_scalar_test();
end

# -------------


text:
"""
    struct CellQuadrature <: GridapType
      array
    end
"""
struct CellQuadrature <: GridapType
  array
  coords
  weights
  @doc """
      CellQuadrature(array::AbstractArray{<:Quadrature})
  """
  function CellQuadrature(array::AbstractArray{<:Quadrature})
    coords = _get_coordinates(array)
    weights = _get_weights(array)
    new(array,coords,weights)
  end
end

"""
    CellQuadrature(degree,polytopes::Vector{<:Polytope}, cell_types::AbstractVector)
"""
function CellQuadrature(degree,polytopes::Vector{<:Polytope}, cell_types::AbstractVector)
  f = (p) -> Quadrature(p,degree)
  quads = map(f,polytopes)
  array = CompressedArray(quads,cell_types)
  CellQuadrature(array)
end

function CellQuadrature(degree,polytopes::Vector{<:Polytope}, cell_types::Fill)
  ctype = cell_types.value
  p = polytopes[ctype]
  quad = Quadrature(p,degree)
  array = Fill(quad,length(cell_types))
  CellQuadrature(array)
end

"""
    get_array(quad::CellQuadrature)
"""
get_array(quad::CellQuadrature) = quad.array

"""
    get_coordinates(q::CellQuadrature)
"""
function get_coordinates(quad::CellQuadrature)
  quad.coords
end

"""
    get_weights(q::CellQuadrature)
"""
function get_weights(quad::CellQuadrature)
  quad.weights
end

function _get_coordinates(q::AbstractArray{<:Quadrature})
  @notimplemented "Not implemented, since we dont need it"
end

function _get_coordinates(q::CompressedArray{<:Quadrature})
  coords = map(get_coordinates,q.values)
  CompressedArray(coords,q.ptrs)
end

function _get_coordinates(q::Fill{<:Quadrature})
  coords = get_coordinates(q.value)
  Fill(coords,length(q))
end

function _get_coordinates(q::AppendedArray)
  a = _get_coordinates(q.a)
  b = _get_coordinates(q.b)
  lazy_append(a,b)
end

function _get_weights(q::AbstractArray{<:Quadrature})
  @notimplemented "Not implemented, since we dont need it"
end

function _get_weights(q::CompressedArray{<:Quadrature})
  w = map(get_weights,q.values)
  CompressedArray(w,q.ptrs)
end

function _get_weights(q::Fill{<:Quadrature})
  w = get_weights(q.value)
  Fill(w,length(q))
end

function _get_weights(q::AppendedArray)
  a = _get_weights(q.a)
  b = _get_weights(q.b)
  lazy_append(a,b)
end

"""
    integrate(cell_field,cell_map::AbstractArray{<:Field},quad::CellQuadrature)
"""
function integrate(cell_field,cell_map::AbstractArray{<:Field},quad::CellQuadrature)
  q = get_coordinates(quad)
  w = get_weights(quad)
  j = gradient(cell_map)
  _f = convert_to_cell_field(cell_field,cell_map)
  f = to_ref_space(_f)
  @assert length(f) == length(cell_map) "Are you using the right triangulation to integrate?"
  @assert length(f) == length(w) "Are you using the right quadrature to integrate?"
  integrate(get_array(f),q,w,j)
end

function lazy_append(quad1::CellQuadrature,quad2::CellQuadrature)
  array = lazy_append(quad1.array,quad2.array)
  CellQuadrature(array)
end


text:
module MLJ

## METHOD IMPORT

# from the Standard Library:
import Distributed: @distributed, nworkers, pmap
import Pkg
import Pkg.TOML

# from the MLJ universe:
using MLJBase
import MLJBase.save
using MLJEnsembles
using MLJTuning
using MLJModels
using OpenML
using MLJSerialization
import MLJSerialization.save # but not re-exported
using MLJIteration
import MLJIteration.IterationControl

using Tables, CategoricalArrays
import Distributions
import Distributions: pdf, mode
import Statistics, StatsBase, LinearAlgebra, Random
import Random: AbstractRNG, MersenneTwister
using ProgressMeter
using ComputationalResources
using ComputationalResources: CPUProcesses

# to be extended:
import MLJBase: fit, update, clean!, fit!, predict, fitted_params,
                show_as_constructed, ==
import MLJModels: models
import ScientificTypes

## METHOD EXPORT

export MLJ_VERSION

## METHOD RE-EXPORT

# re-export from Random, Statistics, Distributions, CategoricalArrays:
export pdf, logpdf, mode, median, mean, shuffle!, categorical, shuffle,
    levels, levels!, std, support, sampler

# re-exports from (MLJ)ScientificTypesBase via MLJBase
export Scientific, Found, Unknown, Known, Finite, Infinite,
    OrderedFactor, Multiclass, Count, Continuous, Textual,
    Binary, ColorImage, GrayImage, Image, Table
export scitype, scitype_union, elscitype, nonmissing, trait
export coerce, coerce!, autotype, schema, info

# re-export from MLJBase:
export nrows, color_off, color_on,
    selectrows, selectcols, restrict, corestrict, complement,
    target_scitype, input_scitype, output_scitype, load_path, training_losses,
    predict, predict_mean, predict_median, predict_mode, predict_joint,
    transform, inverse_transform, evaluate, fitted_params, params,
    @constant, @more, HANDLE_GIVEN_ID, UnivariateFinite,
    classes, table, report, rebind!,
    partition, unpack,
    default_measure, measures,
    @load_boston, @load_ames, @load_iris, @load_reduced_ames, @load_crabs,
    load_boston, load_ames, load_iris, load_reduced_ames, load_crabs,
    Machine, machine, AbstractNode, @node,
    source, node, fit!, freeze!, thaw!, Node, sources, origins,
    machines, sources, anonymize!, @from_network, fitresults,
    @pipeline, Stack,
    ResamplingStrategy, Holdout, CV, TimeSeriesCV,
    StratifiedCV, evaluate!, Resampler, iterator,
    default_resource, pretty,
    make_blobs, make_moons, make_circles, make_regression,
    fit_only!, return!, int, decoder

# MLJBase/composition/abstract_types.jl:
for T in vcat(MLJBase.MLJModelInterface.ABSTRACT_MODEL_SUBTYPES,
    MLJBase.COMPOSITE_TYPES,
    MLJBase.SURROGATE_TYPES)
    @eval(export $T)
end
export Surrogate, Composite
export DeterministicNetwork, ProbabilisticNetwork, UnsupervisedNetwork

# MLJBase/measure/measures.jl:
export orientation, reports_each_observation,
    is_feature_dependent, aggregation,
    aggregate, default_measure, value,
    supports_class_weights, prediction_type, human_name

# MLJBase/measures/continuous.jl:
export mav, mae, mape, rms, rmsl, rmslp1, rmsp, l1, l2, log_cosh,
    MAV, MAE, MeanAbsoluteError, mean_absolute_error, mean_absolute_value,
    LPLoss, RootMeanSquaredProportionalError, RMSP,
    RMS, rmse, RootMeanSquaredError, root_mean_squared_error,
    RootMeanSquaredLogError, RMSL, root_mean_squared_log_error, rmsl, rmsle,
    RootMeanSquaredLogProportionalError, rmsl1, RMSLP,
    MAPE, MeanAbsoluteProportionalError, log_cosh_loss, LogCosh, LogCoshLoss

# MLJBase/measures/confusion_matrix.jl:
export confusion_matrix, confmat, ConfusionMatrix

# MLJBase/measures/finite.jl:
export cross_entropy, BrierScore, brier_score,
    BrierLoss, brier_loss,
    LogLoss, log_loss,
    misclassification_rate, mcr, accuracy,
    balanced_accuracy, bacc, bac, BalancedAccuracy,
    matthews_correlation, mcc, MCC, AUC, AreaUnderCurve,
    MisclassificationRate, Accuracy, MCR, BACC, BAC,
    MatthewsCorrelation

# MLJBase/measures/finite.jl -- Multiclass{2} (order independent):
export auc, area_under_curve, roc_curve, roc

# MLJBase/measures/finite.jl -- OrderedFactor{2} (order dependent):
export TruePositive, TrueNegative, FalsePositive, FalseNegative,
    TruePositiveRate, TrueNegativeRate, FalsePositiveRate,
    FalseNegativeRate, FalseDiscoveryRate, Precision, NPV, FScore,
    NegativePredictiveValue,
    # standard synonyms
    TPR, TNR, FPR, FNR, FDR, PPV,
    Recall, Specificity, BACC,
    # instances and their synonyms
    truepositive, truenegative, falsepositive, falsenegative,
    true_positive, true_negative, false_positive, false_negative,
    truepositive_rate, truenegative_rate, falsepositive_rate,
    true_positive_rate, true_negative_rate, false_positive_rate,
    falsenegative_rate, negativepredictive_value,
    false_negative_rate, negative_predictive_value,
    positivepredictive_value, positive_predictive_value,
    tpr, tnr, fpr, fnr,
    falsediscovery_rate, false_discovery_rate, fdr, npv, ppv,
    recall, sensitivity, hit_rate, miss_rate,
    specificity, selectivity, f1score, fallout

# MLJBase/measures/finite.jl -- Finite{N} - multiclass generalizations of
# above OrderedFactor{2} measures (but order independent):
export MulticlassTruePositive, MulticlassTrueNegative, MulticlassFalsePositive,
      MulticlassFalseNegative, MulticlassTruePositiveRate,
      MulticlassTrueNegativeRate, MulticlassFalsePositiveRate,
      MulticlassFalseNegativeRate, MulticlassFalseDiscoveryRate,
      MulticlassPrecision, MulticlassNegativePredictiveValue, MulticlassFScore,
      # standard synonyms
      MTPR, MTNR, MFPR, MFNR, MFDR, MPPV,
      MulticlassRecall, MulticlassSpecificity,
      # instances and their synonyms
      multiclass_truepositive, multiclass_truenegative,
      multiclass_falsepositive,
      multiclass_falsenegative, multiclass_true_positive,
      multiclass_true_negative, multiclass_false_positive,
      multiclass_false_negative, multiclass_truepositive_rate,
      multiclass_truenegative_rate, multiclass_falsepositive_rate,
      multiclass_true_positive_rate, multiclass_true_negative_rate,
      multiclass_false_positive_rate, multiclass_falsenegative_rate,
      multiclass_negativepredictive_value, multiclass_false_negative_rate,
      multiclass_negative_predictive_value, multiclass_positivepredictive_value,
      multiclass_positive_predictive_value, multiclass_tpr, multiclass_tnr,
      multiclass_fpr, multiclass_fnr, multiclass_falsediscovery_rate,
      multiclass_false_discovery_rate, multiclass_fdr, multiclass_npv,
      multiclass_ppv, multiclass_recall, multiclass_sensitivity,
      multiclass_hit_rate, multiclass_miss_rate, multiclass_specificity,
      multiclass_selectivity, macro_f1score, micro_f1score,
      multiclass_f1score, multiclass_fallout, multiclass_precision,
      # averaging modes
      no_avg, macro_avg, micro_avg

# MLJBase/measures/loss_functions_interface.jl
export dwd_margin_loss, exp_loss, l1_hinge_loss, l2_hinge_loss, l2_margin_loss,
    logit_margin_loss, modified_huber_loss, perceptron_loss, sigmoid_loss,
    smoothed_l1_hinge_loss, zero_one_loss, huber_loss, l1_epsilon_ins_loss,
    l2_epsilon_ins_loss, lp_dist_loss, logit_dist_loss, periodic_loss,
    quantile_loss

# MLJBase/measures/loss_functions_interface.jl
export DWDMarginLoss, ExpLoss, L1HingeLoss, L2HingeLoss, L2MarginLoss,
    LogitMarginLoss, ModifiedHuberLoss, PerceptronLoss, SigmoidLoss,
    SmoothedL1HingeLoss, ZeroOneLoss, HuberLoss, L1EpsilonInsLoss,
    L2EpsilonInsLoss, LPDistLoss, LogitDistLoss, PeriodicLoss,
    QuantileLoss

# re-export from MLJEnsembles:
export EnsembleModel

# re-export from MLJTuning:
export Grid, RandomSearch, Explicit, TunedModel, LatinHypercube,
    learning_curve!, learning_curve

# re-export from MLJModels:
export models, localmodels, @load, @iload, load, info,
    ConstantRegressor, ConstantClassifier,     # builtins/Constant.jl
    FeatureSelector, UnivariateStandardizer,   # builtins/Transformers.jl
    Standardizer, UnivariateBoxCoxTransformer,
    OneHotEncoder, ContinuousEncoder, UnivariateDiscretizer,
    FillImputer, matching, BinaryThresholdPredictor

# re-export from MLJIteration:
export MLJIteration
export IteratedModel
for control in MLJIteration.CONTROLS
    eval(:(export $control))
end
export IterationControl

# re-export from MLJOpenML
const MLJOpenML = OpenML

export OpenML, MLJOpenML

# re-export from MLJSerialization:
export Save # control, not method

# re-export from ComputaionalResources:
export CPU1, CPUProcesses, CPUThreads

## CONSTANTS

const srcdir = dirname(@__FILE__)

## INCLUDE FILES

include("version.jl")       # defines MLJ_VERSION constant
include("scitypes.jl")      # extensions to ScientificTypesBase.scitype

end # module


text:
# Checks whether x[1]^z[1] * ... x[n]^z[n]*y^[n+1] is a torsion unit
# This can be improved
function _check_relation_mod_torsion(x::Vector{FacElem{nf_elem, AnticNumberField}}, y::FacElem{nf_elem, AnticNumberField}, z::Vector{fmpz}, p::Int = 16)
  (length(x) + 1 != length(z)) && error("Lengths of arrays does not fit")
  r = x[1]^z[1]

  for i in 2:length(x)
    r = r*x[i]^z[i]
  end

  w = r*y^z[length(z)]

  b, _ = is_torsion_unit(w, false, p)
#  b, _ = is_torsion_unit(w)
  return b
end

function _find_rational_relation!(rel::Vector{fmpz}, v::arb_mat, bound::fmpz)
  #push!(_debug, (deepcopy(rel), deepcopy(v), deepcopy(bound)))
  @vprint :UnitGroup 2 "Finding rational approximation in $v\n"
  r = length(rel) - 1

  z = Array{fmpq}(undef, r)

  # Compute an upper bound in the denominator of an entry in the relation
  # using Cramer's rule and lower regulator bounds

  # Now comes the rational approximation phase

  # First a trivial check:
  # If the relation contains only integers, it does not yield any information

  i = 0

  is_integer = true

  while is_integer && i < r
    i = i + 1
    b, o = unique_integer(v[1, i])
    if b
      rel[i] = o
    end
    is_integer = is_integer && b
  end

  if is_integer
    rel[r + 1] = -1
    @vprint :UnitGroup 2 "Found rational relation: $rel.\n"
    return true
  end

  for i in 1:r
    if radius(v[1, i]) > 1
      # This is a strange case I cannot handle at the moment
      return false
    end

    app = simplest_inside(v[1, i], bound)

    if app[1]
      z[i] = app[2]
    else
      @vprint :UnitGroup 2 "Something went wrong with the approximation.\n"
      return false
    end
  end

  dlcm = denominator(z[1])

  for i in 2:length(z)
    dlcm = lcm(dlcm, denominator(z[i]))
  end

  if dlcm > bound
    return false
  end

  for i in 1:r
    rel[i] = numerator(z[i]*dlcm)
  end

  rel[r + 1] = -dlcm

  # Check that relation is primitive
  g = rel[1]

  for i in 2:length(rel)
    g = gcd(g, rel[i])
    if g == 1
      break
    end
  end

  @assert g == 1

  @vprint :UnitGroup 2 "Found rational relation: $rel.\n"
  return true
end

# Given r elements x_1,...,x_r, where r is the unit rank, and y an additional unit,
# compute a basis z_1,...z_r such that <x_1,...x_r,y,T> = <z_1,...,z_r,T>,
# where T are the torsion units

function _find_relation(x::Vector{S}, y::T, p::Int = 64) where {S, T}

  K = _base_ring(x[1])

  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  R = ArbField(p, cached = false)

  zz = Array{fmpz}(undef, r + 1)

  @vprint :UnitGroup 1 "Computing conjugates log matrix ... \n"
  A = _conj_log_mat_cutoff(x, p)

  Ar = base_ring(A)

  b = ArbMatSpace(Ar, 1, r)()

  conlog = conjugates_arb_log(y, p)

  for i in 1:r
    b[1, i] = conlog[i]
  end

  B = parent(A)()

  # I should do this using lu decomposition and caching
  # The inversion could go wrong,
  # Then we again have to increase the precision

  inv_succesful = false

  try
    @vprint :UnitGroup 1 "Inverting matrix ... \n"
    B = inv(A)
    inv_succesful = true
  catch e
    @vprint :UnitGroup 1 "Cannot invert matrix ... \n"
    rethrow(e)
  end

  v = b*B

  z = Array{fmpq}(undef, r)

  rreg = det(A)

  bound = _denominator_bound_in_relation(rreg, K)

  # Compute an upper bound in the denominator of an entry in the relation
  # using Cramer's rule and lower regulator bounds

  rel = Array{fmpz}(undef, r + 1)
  for i in 1:r+1
    rel[i] = zero(FlintZZ)
  end

  while !inv_succesful || !_find_rational_relation!(rel, v, bound)
    p =  2*p

    inv_succesful = false

    A = _conj_log_mat_cutoff(x, p)

    Ar = base_ring(A)

    b = ArbMatSpace(Ar, 1, r)()

    conlog = conjugates_arb_log(y, p)

    for i in 1:r
      b[1, i] = conlog[i]
    end

    if !inv_succesful
      try
        @vprint :UnitGroup 1 "Inverting matrix ... \n"
        B = inv(A)
        inv_succesful = true
      catch
        @vprint :UnitGroup 1 "Cannot invert matrix. Increasing precision to $(2*p)\n"
      end
    end

    v = b*B
  end

  # Check if it is a relation modulo torsion units!
  @vprint :UnitGroup 1 "Checking relation $rel \n"

  if !_check_relation_mod_torsion(x, y, rel)
    #error("Dirty approximation did not work")
    return _find_relation(x, y, 2*p)
    #rel[r + 1 ] = 0
    #return rel
  end

  @vprint :UnitGroup 1 "Found a valid relation!\n"
  return rel
end

function _denominator_bound_in_relation(rreg::arb, K::AnticNumberField)
  # Compute an upper bound in the denominator of an entry in the relation
  # using Cramer's rule and lower regulator bounds

  arb_bound = rreg * inv(lower_regulator_bound(K))

  # I want to get an upper bound as an fmpz
  return abs_upper_bound(fmpz, arb_bound)
end

function simplest_inside(x::arb, B::fmpz)
  q = simplest_inside(x)
  if denominator(q) < B
    return true, q
  else
    return false, q
  end
end


text:
### Data Whitening Functions

"""
    whiten_data!(X; dims, df)

Compute a whitening transform matrix for centered data matrix `X`. Use `dims=1` for
row-based observations and `dims=2` for column-based observations. The `df` parameter
specifies the effective degrees of freedom.
"""
function whiten_data_chol!(X::Matrix{T}; dims::Integer, df::Integer=size(X,dims)-1) where T
    df > 0 || error("degrees of freedom must be greater than 0")

    n, p = check_dims(X, dims=dims)

    n > p || error("insufficient number of within-class observations to produce a full " *
                   "rank covariance matrix ($(n) observations, $(p) predictors)")

    if dims == 1
        # X = QR ⟹ S = XᵀX = RᵀR
        W⁻¹ = UpperTriangular(qr!(X, Val(false)).R)
    else
        # Xᵀ = LQ ⟹ S = XXᵀ = LLᵀ = RᵀR
        W⁻¹ = LowerTriangular(lq!(X).L)
    end

    broadcast!(/, W⁻¹, W⁻¹, √(df))

    detΣ = det(W⁻¹)^2

    W = try
        inv(W⁻¹)
    catch err
        if isa(err, LAPACKException) || isa(err, SingularException)
            if err.info ≥ 1
                error("rank deficiency detected (collinearity in predictors)")
            end
        end
        throw(err)
    end

    return (W, detΣ)
end

function whiten_data_svd!(X::Matrix{T}, γ::Union{Nothing,T}; dims::Integer,
                          df::Integer=size(X,dims)-1) where T
    n, p = check_dims(X, dims=dims)

    n > p || error("insufficient number of within-class observations to produce a full " *
                   "rank covariance matrix ($(n) observations, $(p) predictors)")

    0 ≤ γ ≤ 1 || throw(DomainError(γ, "γ must be in the interval [0,1]"))

    tol = eps(one(T))*p*maximum(X)

    UDVᵀ = svd!(X, full=false)

    D = UDVᵀ.S

    if γ !== nothing && γ ≠ zero(T)
        # Regularize: Σ = VD²Vᵀ ⟹ Σ(γ) = V((1-γ)D² + (γ/p)trace(D²)I)Vᵀ
        broadcast!(σᵢ -> abs2(σᵢ)/df, D, D)  # Convert data singular values to Σ eigenvalues
        broadcast!(regularize, D, D, mean(D), γ)
        detΣ = prod(D)
        broadcast!(√, D, D)
    else
        detΣ = prod(σᵢ -> abs2(σᵢ)/df, D)
        broadcast!(/, D, D, √(df))
    end

    all(D .> tol) || error("rank deficiency (collinearity) detected with tolerance $(tol)")

    # Whitening matrix
    if dims == 1
        Vᵀ = UDVᵀ.Vt
        Wᵀ = broadcast!(/, Vᵀ, Vᵀ, D)  # in-place diagonal matrix multiply DVᵀ
    else
        U = UDVᵀ.U
        Wᵀ = broadcast!(/, U, U, transpose(D))
    end

    return (copy(transpose(Wᵀ)), detΣ)
end

function whiten_cov_chol!(Σ::AbstractMatrix{T}, γ::Union{Nothing,T};
                          dims::Integer=1) where T
    p, p₂ = check_dims(Σ, dims=dims)

    p == p₂ || throw(DimensionMismatch("Σ must be square"))

    if γ !== nothing && γ ≠ zero(T)
        zero(T) ≤ γ ≤ one(T) || throw(DomainError(γ, "γ must be in the interval [0,1]"))
        regularize!(Σ, γ)
    end

    UᵀU = cholesky!(Σ, Val(false); check=true)

    if dims == 1
        Lᵀ = UᵀU.U
        detΣ = det(Lᵀ)^2

        return (inv(Lᵀ), detΣ)
    else
        L = UᵀU.L
        detΣ = det(L)^2

        return (inv(L), detΣ)
    end
end


text:
PressureField(model, data=model.pressures.pHY′.data; kw...) =
    Field(model.pressures.pHY′ + model.pressures.pNHS; data, kw...)


text:
# # Solid body rotation of a meridional sector on the sphere
#
# This script implements the "Global Steady State Nonlinear Zonal Geostrophic Flow"
# validation experiment from
#
# > Williamson et al., "A Standard Test Set for Numerical Approximations to the Shallow
#   Water Equations in Spherical Geometry", Journal of Computational Physics, 1992.
#
# The problem is posed in spherical strip between 60ᵒS and 60ᵒN latitude on a sphere with
# unit radius.
#
# # Dependencies
#
# The validation experiment depends on Oceananigans, JLD2, Printf, and GLMakie for visualization

using Oceananigans
using Oceananigans.Grids

using Oceananigans.Coriolis: HydrostaticSphericalCoriolis

using Oceananigans.Advection:
    EnergyConservingScheme,
    EnstrophyConservingScheme

using Oceananigans.Models.HydrostaticFreeSurfaceModels:
    HydrostaticFreeSurfaceModel,
    VectorInvariant,
    ExplicitFreeSurface,
    PrescribedVelocityFields

using Oceananigans.Utils: prettytime, hours
using Oceananigans.OutputWriters: JLD2OutputWriter, TimeInterval, IterationInterval

using JLD2
using Printf
using GLMakie

# # The geostrophic flow
#
# ```math
# u = U \cos ϕ
# v = 0
# η = - g^{-1} \left (R Ω U + \frac{U^2}{2} \right ) \sin^2 ϕ
# ```
#
# is a steady nonlinear flow on a sphere of radius ``R`` with gravitational
# acceleration ``g``, corresponding to solid body rotation
# in the same direction as the "background" rotation rate ``\Omega``.
# The velocity ``U`` determines the magnitude of the additional rotation.

const U = 0.1

solid_body_rotation(λ, ϕ) = U * cosd(ϕ)
solid_body_geostrophic_height(λ, ϕ, R, Ω, g) = (R * Ω * U + U^2 / 2) * sind(ϕ)^2 / g

# In addition to the solid body rotation solution, we paint a Gaussian tracer patch
# on the spherical strip to visualize the rotation.

northern_boundary = 80 # degrees
Ω = 1 # rad / s
g = 1 # m s⁻²

function run_solid_body_tracer_advection(; architecture = CPU(),
                                           Nx = 360,
                                           Ny = 8,
                                           super_rotations = 4)

    # A spherical domain
    @show grid = LatitudeLongitudeGrid(architecture, size = (Nx, Ny, 1),
                                       radius = 1,
                                       latitude = (-northern_boundary, northern_boundary),
                                       longitude = (-180, 180),
                                       z = (-1, 0))

    uᵢ(λ, ϕ, z, t=0) = solid_body_rotation(λ, ϕ)

    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        tracers = (:c, :d, :e),
                                        velocities = PrescribedVelocityFields(u=uᵢ),
                                        coriolis = nothing,
                                        buoyancy = nothing,
                                        closure = nothing)

    # Tracer patch for visualization
    Gaussian(λ, ϕ, L) = exp(-(λ^2 + ϕ^2) / 2L^2)

    # Tracer patch parameters
    L = 24 # degree
    ϕ₀ = 0 # degrees

    cᵢ(λ, ϕ, z) = Gaussian(λ, 0, L)
    dᵢ(λ, ϕ, z) = Gaussian(0, ϕ - ϕ₀, L)
    eᵢ(λ, ϕ, z) = Gaussian(λ, ϕ - ϕ₀, L)

    set!(model, c=cᵢ, d=dᵢ, e=eᵢ)

    ϕᵃᶜᵃ_max = maximum(abs, ynodes(Center, grid))
    Δx_min = grid.radius * cosd(ϕᵃᶜᵃ_max) * deg2rad(grid.Δλ)
    Δy_min = grid.radius * deg2rad(grid.Δϕ)
    Δ_min = min(Δx_min, Δy_min)

    # Time-scale for tracer advection across the smallest grid cell
    @show advection_time_scale = Δ_min / U
    super_rotation_period = 2π * grid.radius / U

    simulation = Simulation(model,
                            Δt = 0.1advection_time_scale,
                            stop_time = super_rotations * super_rotation_period,
                            iteration_interval = 100,
                            progress = s -> @info "Time = $(s.model.clock.time) / $(s.stop_time)")

    output_fields = model.tracers

    output_prefix = "solid_body_tracer_advection_Nx$(grid.Nx)"

    simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                          schedule = TimeInterval(super_rotation_period / 20),
                                                          prefix = output_prefix,
                                                          overwrite_existing = true)

    run!(simulation)

    return simulation.output_writers[:fields].filepath
end

function visualize_solid_body_tracer_advection(filepath)

    @show output_prefix = basename(filepath)[1:end-5]

    file = jldopen(filepath)

    iterations = parse.(Int, keys(file["timeseries/t"]))

    Nx = file["grid/Nx"]
    Ny = file["grid/Ny"]

    grid = LatitudeLongitudeGrid(size = (Nx, Ny, 1),
                                 radius = 1,
                                 latitude = (-northern_boundary, northern_boundary),
                                 longitude = (-180, 180),
                                 z = (-1, 0))

    super_rotation_period = 2π * grid.radius / U

    λ = xnodes(Face, grid)
    ϕ = ynodes(Center, grid)

    λ = repeat(reshape(λ, Nx, 1), 1, Ny)
    ϕ = repeat(reshape(ϕ, 1, Ny), Nx, 1)

    λ_azimuthal = λ .+ 180  # Convert to λ ∈ [0°, 360°]
    ϕ_azimuthal = 90 .- ϕ   # Convert to ϕ ∈ [0°, 180°] (0° at north pole)

    iter = Node(0)

    plot_title = @lift @sprintf("Tracer advection by geostrophic solid body rotation: rotations = %.3f",
                                file["timeseries/t/" * string($iter)] / super_rotation_period)

    c = @lift file["timeseries/c/" * string($iter)][:, :, 1]
    d = @lift file["timeseries/d/" * string($iter)][:, :, 1]
    e = @lift file["timeseries/e/" * string($iter)][:, :, 1]

    # Plot on the unit sphere to align with the spherical wireframe.
    x = @. cosd(λ_azimuthal) * sind(ϕ_azimuthal)
    y = @. sind(λ_azimuthal) * sind(ϕ_azimuthal)
    z = @. cosd(ϕ_azimuthal)

    fig = Figure(resolution = (1080, 1080))

    titles = ["c", "d", "e"]

    for (n, var) in enumerate((c, d, e))
        ax = fig[n, 1] = Axis(fig, xlabel = "λ", ylabel = "ϕ", title = titles[n])
        heatmap!(ax, var)
    end

    supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)

    record(fig, output_prefix * ".mp4", iterations, framerate=30) do i
        @info "Animating iteration $i/$(iterations[end])..."
        iter[] = i
    end

    return nothing
end

filepath = run_solid_body_tracer_advection(Nx=180, Ny=8, super_rotations=2)
visualize_solid_body_tracer_advection(filepath)


text:
using Thebes, Luxor

function makecube()

    cube = [
    Point3D(1,   1, -1),
    Point3D(1,  -1, -1),
    Point3D(-1, -1, -1),
    Point3D(-1,  1, -1),
    Point3D(1,   1,  1),
    Point3D(1,  -1,  1),
    Point3D(-1, -1,  1),
    Point3D(-1,  1,  1)]

    r = Point3D[]

    for e in [
        [1, 2], [2, 3], [3, 4], [4, 1],
        [5, 6], [6, 7], [7, 8], [8, 5],

        [1, 5], [5, 8], [8, 4], [4, 1],
        [1, 2], [2, 6], [6, 5], [5, 1],
        [3, 2], [2, 6], [6, 7], [7, 3],
        [4, 8], [8, 7], [7, 3], [3, 4],
        ]
        push!(r, cube[e][1])
        push!(r, cube[e][2])
    end
    return r
end

@draw begin
    setline(1.5)
    background("midnightblue")
    sethue("darkorange")
    eyepoint(Point3D(500, 500, 500))
    perspective(0)
    for θ in 0:π/48:4π
        p = Point3D(100cos(θ), 100sin(θ), 20θ)

        #=
        pin(p, gfunction = (p3, p2) -> begin
            d = distance(p3, eyepoint())
            circle(p2, rescale(d, 0, 300, 20, 5), :fill)
        end
        )
        =#

        #=
        pin(p, Point3D(50cos(θ), 50sin(θ), p.z), gfunction = (p3s, p2s) -> begin
            line(p2s..., :stroke)
        end)
        =#

        # pin(p, between(p, eyepoint(), 0.8))

        helix = [Point3D(100cos(θ), 100sin(θ), 20θ) for θ in 0:π/12:4π]
        eyepoint(Point3D(120, 120, 120))
        a_box =
        pin(helix,
        gfunction = (p3list, p2list) -> prettypoly(p2list, :stroke, () ->
        begin

            randomhue()
            circle(O, 2, :fill)
        end)
        )
    end

    axes3D()

    setline(0.5)

    #draw cubes, moving each one and scaling it
    for x in -200:40:200
        for y in -200:40:200
            randomhue()
            pin(moveby!(15makecube(), Point3D(x, y, 0)))
        end
    end

end


text:
using TrajectoryOptimization
using RigidBodyDynamics
using Test
using BenchmarkTools
using LinearAlgebra
using Random
using SparseArrays
using ForwardDiff
using Logging
using PartedArrays
using Ipopt
const TO = TrajectoryOptimization

disable_logging(Logging.Debug)
@testset "Logging" begin
    include("logger_tests.jl")
end
disable_logging(Logging.Info)

@testset "Constraints" begin
    include("constraint_tests.jl")
end
@testset "Model" begin
    include("model_tests.jl")
end
@testset "Modified Model" begin
    include("modified_model_tests.jl")
end
@testset "Problems" begin
    include("problem_tests.jl")
end
@testset "Square Root Backward Pass" begin
    include("sqrt_bp_tests.jl")
end
@testset "Infeasible State Trajectory Initialization" begin
    include("infeasible_tests.jl")
end
@testset "Minimum Time" begin
    include("minimum_time_tests.jl")
end
@testset "Costs" begin
    include("cost_tests.jl")
end
@testset "Utils" begin
    include("test_utils.jl")
end

# Systems
@testset "Pendulum" begin
    include("pendulum_tests.jl")
end
@testset "Car" begin
    include("car_tests.jl")
end
@testset "Quadrotor" begin
    include("quadrotor_tests.jl")
end

# Direct Methods
@testset "DIRCOL" begin
    include("dircol_test.jl")
end
disable_logging(Logging.Debug)

#= Tests needed
* Generic Cost
* Expansion multiply and copy
* Quad Cost constructor exceptions
* AugmentedLagrangianObjective constructors
* Sqrt bp exceptions
* Dynamics module?
* iLQR not inplace solves
* iLQR live plotting
* feedforward gradient
* convergence criteria (max iters, dJ zero, max reg)
* infeasible mintime
* logger: trim entry
* logger: print header
* evaluate continuous dynamics
* simultaneous jacobian! function in model
* jacobian! on continuous dynamics for entire trajectory
* model from urdf string (normal and underactuated)
* wrap_inplace
* Objective from cost and terminal cost
* Objective from cost trajectory and terminal cost
* Pass in costfuntion directly to problem
* _validate_time
* not inplace rollout
* copy iLQR solver
=#

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using ODEs
using Base.Test

# write your own tests here
f(t,y) = [4*exp(0.8*t)-0.5*y[1],4*exp(0.8*t)-0.5*y[2]]

println("Testing solver: rk4")
t, y = rk4(f, [0 4], [2 2])

@test((y[end,1]-75.33896)/75.33896 < 0.01)
print_with_color(:green,"Solver rk4 passed\n")

println("Testing solver: rk5")
t, y = rk4(f, [0 4], [2 2])

@test((y[end,1]-75.33896)/75.33896 < 0.01)
print_with_color(:green,"Solver rk5 passed\n")

println("All looks OK")


text:
using DiskCaches
using Test
using Base.Threads
using Pidfile

@testset "DiskCaches.jl" begin
    path = mktempdir()
    c1 = DiskCache(joinpath(path, "path1"))
    @test length(c1) == 0
    @test isempty(c1) == true
    @test_throws ErrorException empty!(c1)
    @test (c1[1] = 1) == 1
    @test_throws ErrorException c1[1] = 2
    @test_throws ErrorException pop!(c1)
    @test_throws ErrorException delete!(c1, 1)
    @test_throws KeyError c1[2]
    @test get(c1, 2, "hello") == "hello"
    @test get(c1, 2, "goodbye") == "goodbye"
    @test get!(c1, 2, "hello") == "hello"
    @test get!(c1, 3) do
        [1, 2, [3, "goodbye"]]
    end == [1, 2, [3, "goodbye"]]
    @test c1[2] == "hello"
    @test c1[3] == [1, 2, [3, "goodbye"]]
    @test length(c1) == 3
    @test !isempty(c1)

    c1_copy = DiskCache(joinpath(path, "path1"))
    @test c1_copy[1] == 1
    @test c1_copy[2] == "hello"
    @test c1_copy[3] == [1, 2, [3, "goodbye"]]
    @test_throws ErrorException c1_copy[3] = 4
    c1_copy[4] = 4
    @test c1_copy[4] == 4
    @test c1[4] == 4
    c1_copy[5] = 5
    @test length(c1) == 5

    c2 = DiskCache(joinpath(path, "path2"))
    c2_copy = DiskCache(joinpath(path, "path2"))
    c2[1] = 4
    c2[2] = 5
    c2[3] = 6
    @test Set([k for k in keys(c2_copy)]) == Set([1, 2, 3])
    @test Set([v for v in values(c2_copy)]) == Set([4, 5, 6])
    @test Set([(k, v) for (k, v) in c2_copy]) == Set([(1, 4), (2, 5), (3, 6)])
    @test Set([(k, v) for (k, v) in pairs(c2_copy)]) == Set([(1, 4), (2, 5), (3, 6)])

    C3 = [DiskCache(joinpath(path, "path3")) for _ in 1:16]
    ncalls = Atomic{Int}(0)
    @threads for i = 1:16
        get!(C3[i], "$(i)_get!", "get!")
        get!(C3[i], "$(rand(1:16))_get_call!") do
            "get_call!"
            atomic_add!(ncalls, 1)
        end
        get!(C3[i], "$(i)_get_call!") do
            "get_call!"
            atomic_add!(ncalls, 1)
        end

        C3[i]["$(i)_set!"] = "set!"
    end
    @test all(length.(C3) .== 3 * 16)
    @test ncalls[] == 16
end


text:
using Random, Distributions, DelimitedFiles
Random.seed!()

# Here I define the global quantities used by all pivot MCMC run which are the pivot moves
const pivot_moves = Array{Int,3}(zeros(47, 2, 3)) # Tensor that defines all the pivot moves of the octahedral symmetry group
const perms = [[1 2 3]; [1 3 2]; [2 1 3]; [2 3 1]; [3 1 2]; [3 2 1]]
const tr_signs = [[1 1 1]; [1 1 -1]; [1 -1 1]; [1 -1 -1]; [-1 1 1]; [-1 1 -1]; [-1 -1 1]; [-1 -1 -1]]
p_moves_count = 1
for i in 1:6
    for j in 1:8
        if i==1 && j==1
            continue
        else
            pivot_moves[p_moves_count, 1, :] .= perms[i,:]
            pivot_moves[p_moves_count, 2, :] .= tr_signs[j,:]

            global p_moves_count += 1
        end
    end
end

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

# This structure contains all I need to compute the Hamiltonian
struct Magnetic_polymer{T1<:Int, T2<:AbstractFloat}
    n_mono::T1
    inv_temp::T2
    J::T2
    alpha::T2
    energy::T2
    spins::Array{T1}
    fields::Array{T2}
    coord::Array{T1,2}         # Coordinates of the monomers
    trial_coord::Array{T1,2}
    hash_saw::Dict{Int64, Int64}
    neighbours::Array{T1, 2}
    trial_neighbours::Array{T1, 2}
end

function Magnetic_polymer(n_mono::T1, inv_temp::T2, J::T2, alpha::T2) where {T1<:Int, T2<:AbstractFloat}
    Magnetic_polymer(n_mono, inv_temp, J, alpha, 0.0, zeros(Int, n_mono), zeros(n_mono), zeros(Int, 3,n_mono),
    zeros(Int, 3,n_mono), Dict{Int, Int}(), zeros(Int, 7,n_mono), zeros(Int, 7,n_mono))
end

function initialize_poly!(poly::Magnetic_polymer, spin_values::Array{Int})
    poly.spins .= rand(spin_values, poly.n_mono)
    #poly.fields .= (rand(poly.n_mono).*2) .-1
    for i_mono in 1:poly.n_mono
        poly.coord[1, i_mono] = i_mono-1
    end
    for i_mono in 1:poly.n_mono
        for j in 1:3
            poly.trial_coord[j, i_mono] = poly.coord[j, i_mono]
        end
    end
end

function initialize_poly!(poly::Magnetic_polymer)
    initialize_poly!(poly,[0,1])
end

#The following method of the "initialize_poly!()" function allows initialization from file
function initialize_poly!(poly::Magnetic_polymer, file_name::String)
    io = open(file_name, "r")
    data = readdlm(io,Int64)
    close(io)
    for i_mono in 1:poly.n_mono
        for j in 1:3
            poly.coord[j,i_mono] = data[i_mono,j]
        end
        poly.spins[i_mono] = data[i_mono, 4]
    end
end

function set_fields!(polymers::Array{Magnetic_polymer}, ff::Array{Float64})
    for polymer in polymers
        for i_mono in 1:polymer.n_mono
            polymer.fields[i_mono] = ff[i_mono]
        end
    end
end

# This structure stores series of observables I'm computing along the MCMC
struct MC_data{T1<:Int, T2<:AbstractFloat}
    n_steps::T1
    energies::Array{T2}
    magnetization::Array{T2}
    rg2::Array{T2}
end

function MC_data(n_steps::Int)
    MC_data(n_steps, zeros(n_steps), zeros(n_steps), zeros(n_steps))
end

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

function try_pivot!(k::Int, move::Int, coo::Array{Int,2}, t_coo::Array{Int,2}, n_mono::Int)
    for i_mono in 1:k
        for j in 1:3
            #pol.trial_coord[j, i_mono] = pol.coord[j, i_mono]
            t_coo[j, i_mono] = coo[j, i_mono]
        end
    end
    for i_mono in k+1:n_mono
        for j in 1:3
            t_coo[j,i_mono] = (pivot_moves[move,2,j]*(coo[pivot_moves[move,1,j],i_mono]
                                         - coo[pivot_moves[move,1,j],k]) + coo[j,k])
        end
    end
end

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

function checksaw!(k::Int, pol::Magnetic_polymer, a::Int)
    L = max(k, pol.n_mono-k) + 1
    is_saw = true
    t = 0
    while(t<L && is_saw)
        if k+t <= pol.n_mono
            saw_key = pol.trial_coord[1,k+t]*a^2 + pol.trial_coord[2,k+t]*a + pol.trial_coord[3,k+t]
            if haskey(pol.hash_saw, saw_key)
                is_saw = false
            else
                pol.hash_saw[saw_key] = k+t
            end
        end
        if k-t >= 1 && t!=0
            saw_key = pol.trial_coord[1,k-t]*a^2 + pol.trial_coord[2,k-t]*a + pol.trial_coord[3,k-t]
            if haskey(pol.hash_saw, saw_key)
                is_saw = false
            else
                pol.hash_saw[saw_key] = k-t
            end
        end
        t += 1
    end
    return is_saw
end

#########################################################################
## poi quando si faranno le mosse locali non si userà questa funnzione perché
## non ci sarà bisogno di ricalcolare TUTTA la lista dei vicini
function compute_neighbours!(coo::Array{Int,2}, near::Array{Int,2}, dic::Dict{Int,Int}, a::Int, n_mono::Int)
    for i_mono in 1:n_mono
        near[1, i_mono] = 0
        for j in 2:7
            near[j, i_mono] = -1
        end
        saw_key_base = coo[1,i_mono]*a^2 + coo[2,i_mono]*a + coo[3,i_mono]
        for j in 1:3
            for k in -1:2:1
                #coo[j,i_mono] += k
                #saw_key = coo[1,i_mono]*a^2 + coo[2,i_mono]*a + coo[3,i_mono]
                saw_key = saw_key_base + k*a^(3-j)
                if haskey(dic, saw_key)
                    near[1, i_mono] += 1
                    n_neigh = near[1, i_mono]+1
                    near[n_neigh, i_mono] = dic[saw_key]
                end
                #coo[j,i_mono] -= k
            end
        end
    end
end

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

function gyration_radius_squared(pol::Magnetic_polymer)
    rcm = zeros(3)
    rg2 = 0.0
    for i in 1:pol.n_mono
        for j in 1:3
            #rg2 += sum((view(pol.coord,:,i) .- rcm).^2)
            rg2 += pol.coord[j,i]^2
            rcm[j] += pol.coord[j,i]
        end
    end
    r = 0.0
    for j in 1:3
        r += rcm[j]^2
    end
    return rg2/pol.n_mono - r/(pol.n_mono)^2
end

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

function compute_new_energy(pol::Magnetic_polymer, near::Array{Int,2}, ff::Array{Float64}, s::Array{Int})
    ene = 0.0
    ene_J = 0.0
    for i_mono in 1:pol.n_mono
        ene -= ff[i_mono] * s[i_mono]
        for j in 1:near[1, i_mono]
            ene_J -= s[i_mono] * s[near[j+1,i_mono]] * pol.J * 0.5
        end
    end
    return ene*(1-pol.alpha) + ene_J*pol.alpha
end

#########################################################################
#########################################################################
### In this section I there are the functions to perform local moves on the self-avoiding walk

function single_bead_flip!(t_coo::Array{Int,2}, dic::Dict{Int, Int}, n_mono::Int, a::Int)
    mono = rand(2:n_mono-1)
    dist = 0
    for j in 1:3
        dist += (t_coo[j,mono+1]-t_coo[j,mono-1])^2
    end
    if dist == 2
        saw_key = 0
        for j in 1:3
            saw_key += (t_coo[j,mono-1]+t_coo[j,mono+1]-t_coo[j,mono])*a^(3-j)
        end

        if !haskey(dic,saw_key)
            delete!(dic, t_coo[1,mono]*a^2 + t_coo[2,mono]*a + t_coo[3,mono])
            for j in 1:3
                t_coo[j,mono] = t_coo[j,mono-1]+t_coo[j,mono+1]-t_coo[j,mono]
            end
            dic[saw_key] = mono
        end
    end
end

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

function crankshaft_180!(t_coo::Array{Int,2}, dic::Dict{Int, Int}, n_mono::Int, a::Int)
    mono = rand(1:n_mono-3)
    dist = 0
    for j in 1:3
        dist += (t_coo[j,mono+3]-t_coo[j,mono])^2
    end
    if dist==1
        saw_key1 = 0
        saw_key2 = 0
        for j in 1:3
            saw_key1 += (2*t_coo[j,mono] - t_coo[j,mono+1])*a^(3-j)
            saw_key2 += (2*t_coo[j,mono+3] - t_coo[j,mono+2])*a^(3-j)
        end
        if !haskey(dic,saw_key1) && !haskey(dic,saw_key2)
            delete!(dic, t_coo[1,mono+1]*a^2 + t_coo[2,mono+1]*a + t_coo[3,mono+1])
            delete!(dic, t_coo[1,mono+2]*a^2 + t_coo[2,mono+2]*a + t_coo[3,mono+2])
            for j in 1:3
                t_coo[j,mono+1] = 2*t_coo[j,mono] - t_coo[j,mono+1]
                t_coo[j,mono+2] = 2*t_coo[j,mono+3] - t_coo[j,mono+2]
            end
            dic[saw_key1] = mono+1
            dic[saw_key2] = mono+2
        end
    end
end

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

function crankshaft_90_270!(t_coo::Array{Int,2}, dic::Dict{Int, Int}, n_mono::Int, a::Int, new_coord1::Array{Int}, new_coord2::Array{Int})
    mono = rand(1:n_mono-3)
    dist = 0
    axis = 0
    for j in 1:3
        dist += (t_coo[j,mono+3]-t_coo[j,mono])^2
    end
    if dist == 1
        orient = 0
        for j in 1:3
            if (t_coo[j,mono+3]-t_coo[j,mono]) != 0
                axis = j
                orient = t_coo[j,mono+3]-t_coo[j,mono]
                break
            end
        end
        p=0
        t=0
        if axis==1
            p=2
            t=rand((3,2))
        elseif axis==2
            p=6
            t=rand((2,5))
        elseif axis==3
            p=3
            t=rand((5,3))
        else
            println("Invalid Rotation Axis!")
        end

        for j in 1:3
            new_coord1[j] = t_coo[j,mono] + tr_signs[t,j]*(t_coo[perms[p,j],mono+1]-t_coo[perms[p,j],mono])*orient
            new_coord2[j] = t_coo[j,mono+3] + tr_signs[t,j]*(t_coo[perms[p,j],mono+2]-t_coo[perms[p,j],mono+3])*orient
        end

        saw_key1 = new_coord1[1]*a^2 + new_coord1[2]*a + new_coord1[3]
        saw_key2 = new_coord2[1]*a^2 + new_coord2[2]*a + new_coord2[3]
        if !haskey(dic,saw_key1) && !haskey(dic,saw_key2)
            delete!(dic, t_coo[1,mono+1]*a^2 + t_coo[2,mono+1]*a + t_coo[3,mono+1])
            delete!(dic, t_coo[1,mono+2]*a^2 + t_coo[2,mono+2]*a + t_coo[3,mono+2])
            for j in 1:3
                t_coo[j,mono+1] = new_coord1[j]
                t_coo[j,mono+2] = new_coord2[j]
            end
            dic[saw_key1] = mono+1
            dic[saw_key2] = mono+2
        end
    end
end

#########################################################################
#########################################################################
function spins_MC!(pol::Magnetic_polymer, n_flips::Int, ff::Array{Float64}, s::Array{Int}, near::Array{Int,2})
    delta_tot = 0.0
    for i_flip in 1:n_flips
        flip_candidate = rand(1:pol.n_mono)
        local_ene = 0.0
        local_ene -= ff[flip_candidate] * s[flip_candidate] * (1-pol.alpha)
        for j in 1:near[1,flip_candidate]
            local_ene -= pol.J * pol.alpha * s[flip_candidate] * s[near[j+1,flip_candidate]]
        end

        #s[flip_candidate]==1 ? trial_spin_value=-1 : trial_spin_value=1
        s[flip_candidate]==1 ? trial_spin_value=0 : trial_spin_value=1

        trial_local_ene = 0.0
        trial_local_ene -= ff[flip_candidate] * trial_spin_value * (1-pol.alpha)
        for j in 1:near[1,flip_candidate]
            trial_local_ene -= pol.J * pol.alpha * trial_spin_value * s[near[j+1,flip_candidate]]
        end

        delta_ene = trial_local_ene - local_ene
        acc = 0.0

        delta_ene<=0 ? acc=1.0 : acc=exp(-delta_ene*pol.inv_temp)

        if acc > rand()
            delta_tot += delta_ene
            s[flip_candidate] = trial_spin_value
        end
    end
    return delta_tot
end

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

#= This is the version of MC_run where I don't save the configurations of the spins (only the last one)
   This is useful if I just want to perform the simulation and I am interested only in storing the
    values of some observables along the markov chain. If instead I want to carry out some simulator
    based inference task, I will also need to save spin configurations. Expoiting julia's
    multiple dispatch feature I can always use the same "MC_run!()" function, the specific method
    is inferred by the compiler depending on the list of argument types. "MC_run_save!()" is the function
    that save the configurations and it will also be made into a MC_run!() module=#

function MC_run_base!(pol::Magnetic_polymer, traj::MC_data, start::Int, finish::Int)
    for i_mono in 1:pol.n_mono
        for j in 1:3
            pol.trial_coord[j,i_mono] = pol.coord[j,i_mono]
        end
    end
    new_coord1 = zeros(Int,3)
    new_coord2 = zeros(Int,3)
    n_acc = 0
    n_pivots = 0
    #=
      I represent the possible positions of the monomers in the cubic lattice
      as seeing the vector of coordinates as the representation in base "hash_base".
      e.g. (xi,yi,zi) --> xi*hash_base^2 + yi*hash_base + zi*1
      Since the polymer is fixed in the origin the coordinates are ≤ n_mono in modulus, so if I take hash_base
      >= 2*n_mono + 1, the mapping I think it should be 1-to-1. This speeds up the code significantly as the
      hashing of integer valued keys for the dictionary containing a certain polymer config. is much
      cheaper than hashing keys which are arrays sliced from a matrix
    =#
    hash_base = 2*pol.n_mono + 1

    checksaw!(1, pol, hash_base)
    compute_neighbours!(pol.coord, pol.neighbours, pol.hash_saw, hash_base, pol.n_mono)
    energy = compute_new_energy(pol,pol.neighbours,pol.fields, pol.spins)

    traj.energies[start] = energy
    traj.magnetization[start] = mean(pol.spins)
    traj.rg2[start] = gyration_radius_squared(pol)

    empty!(pol.hash_saw)

    for i_step in (start+1):finish
        pivot = rand(2:pol.n_mono-1)
        p_move = rand(1:47)
        if i_step%50000 == 0
            println("Pivot: ", pivot)
            println("Attmpted move: ", p_move)
            println("step number: ", i_step, "\n")
        end
        try_pivot!(pivot, p_move, pol.coord, pol.trial_coord, pol.n_mono)
        still_saw = checksaw!(pivot, pol, hash_base)

        # If the pivot move was unsuccessful try next move on the previous configuration
        if !still_saw
            empty!(pol.hash_saw)
            for i_mono in 1:pol.n_mono
                for j in 1:3
                    pol.trial_coord[j,i_mono] = pol.coord[j,i_mono]
                end
                pol.hash_saw[pol.trial_coord[1,i_mono]*hash_base^2+pol.trial_coord[2,i_mono]*hash_base+pol.trial_coord[3,i_mono]] = i_mono
            end
        else
            n_pivots += 1
        end

        for i_local in 1:pol.n_mono
            mv = rand(1:4)
            if mv==1
                single_bead_flip!(pol.trial_coord,pol.hash_saw,pol.n_mono,hash_base)
            elseif mv==2
                crankshaft_180!(pol.trial_coord,pol.hash_saw,pol.n_mono,hash_base)
            else
                crankshaft_90_270!(pol.trial_coord,pol.hash_saw,pol.n_mono,hash_base,new_coord1,new_coord2)
            end
        end

        compute_neighbours!(pol.trial_coord, pol.trial_neighbours, pol.hash_saw, hash_base, pol.n_mono)
        trial_energy = compute_new_energy(pol,pol.trial_neighbours,pol.fields, pol.spins)
        delta_energy = trial_energy - energy
        acc = 0.0
        delta_energy<=0 ? acc=1.0 : acc=exp(-delta_energy*pol.inv_temp)
        if acc > rand()
            for i_mono in 1:pol.n_mono
                for j in 1:3
                    pol.coord[j,i_mono] = pol.trial_coord[j,i_mono]
                end
            end
            for i_mono in 1:pol.n_mono
                for j in 1:7
                    pol.neighbours[j,i_mono] = pol.trial_neighbours[j,i_mono]
                end
            end
            energy = trial_energy
            n_acc += 1
        end

        energy += spins_MC!(pol, pol.n_mono,pol.fields,pol.spins,pol.neighbours) ## spins_MC! return the total energy variation of the accepted spin n_flips
        traj.rg2[i_step] = gyration_radius_squared(pol)
        traj.magnetization[i_step] = mean(pol.spins)
        traj.energies[i_step] = energy

        # This next if ... end can be removed. I kept as a sanity check in an early stage
        # when I had problems with the computation of the nearest neighbours
        #=if !isinteger(energy)
            println("Non integer energy:  ", energy)
            println(pol.coord)
            println(pol.neighbours .- pol.trial_neighbours)
            println(pol.trial_neighbours)
            break
        end=#

        empty!(pol.hash_saw)
    end
end

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

# This function is like MC_run_base! but also save configurations
function MC_run_save!(pol::Magnetic_polymer, traj::MC_data, start::Int, finish::Int, spins_configs::Matrix{Int}, sample_lag::Int, n_samples::Int)
    for i_mono in 1:pol.n_mono
        for j in 1:3
            pol.trial_coord[j,i_mono] = pol.coord[j,i_mono]
        end
    end
    new_coord1 = zeros(Int,3)
    new_coord2 = zeros(Int,3)
    n_acc = 0
    n_pivots = 0

    hash_base = 2*pol.n_mono + 1

    checksaw!(1, pol, hash_base)
    compute_neighbours!(pol.coord, pol.neighbours, pol.hash_saw, hash_base, pol.n_mono)
    energy = compute_new_energy(pol,pol.neighbours,pol.fields, pol.spins)

    traj.energies[start] = energy
    traj.magnetization[start] = mean(pol.spins)
    traj.rg2[start] = gyration_radius_squared(pol)

    empty!(pol.hash_saw)

    if start%sample_lag == 0
        i_sample = div(start,sample_lag)
        for i_mono in 1:pol.n_mono
            spins_configs[i_mono, i_sample]  = pol.spins[i_mono]
        end
    end

    for i_step in (start+1):min(finish,sample_lag*n_samples)
        pivot = rand(2:pol.n_mono-1)
        p_move = rand(1:47)
        if i_step%50000 == 0
            println("Pivot: ", pivot)
            println("Attmpted move: ", p_move)
            println("step number: ", i_step, "\n")
        end
        try_pivot!(pivot, p_move, pol.coord, pol.trial_coord, pol.n_mono)
        still_saw = checksaw!(pivot, pol, hash_base)

        # If the pivot move was unsuccessful try next move on the previous configuration
        if !still_saw
            empty!(pol.hash_saw)
            for i_mono in 1:pol.n_mono
                for j in 1:3
                    pol.trial_coord[j,i_mono] = pol.coord[j,i_mono]
                end
                pol.hash_saw[pol.trial_coord[1,i_mono]*hash_base^2+pol.trial_coord[2,i_mono]*hash_base+pol.trial_coord[3,i_mono]] = i_mono
            end
        else
            n_pivots += 1
        end

        for i_local in 1:pol.n_mono
            mv = rand(1:4)
            if mv==1
                single_bead_flip!(pol.trial_coord,pol.hash_saw,pol.n_mono,hash_base)
            elseif mv==2
                crankshaft_180!(pol.trial_coord,pol.hash_saw,pol.n_mono,hash_base)
            else
                crankshaft_90_270!(pol.trial_coord,pol.hash_saw,pol.n_mono,hash_base,new_coord1,new_coord2)
            end
        end

        compute_neighbours!(pol.trial_coord, pol.trial_neighbours, pol.hash_saw, hash_base, pol.n_mono)
        trial_energy = compute_new_energy(pol,pol.trial_neighbours,pol.fields, pol.spins)
        delta_energy = trial_energy - energy
        acc = 0.0
        delta_energy<=0 ? acc=1.0 : acc=exp(-delta_energy*pol.inv_temp)
        if acc > rand()
            for i_mono in 1:pol.n_mono
                for j in 1:3
                    pol.coord[j,i_mono] = pol.trial_coord[j,i_mono]
                end
            end
            for i_mono in 1:pol.n_mono
                for j in 1:7
                    pol.neighbours[j,i_mono] = pol.trial_neighbours[j,i_mono]
                end
            end
            energy = trial_energy
            n_acc += 1
        end

        energy += spins_MC!(pol, pol.n_mono,pol.fields,pol.spins,pol.neighbours) ## spins_MC! return the total energy variation of the accepted spin n_flips
        traj.rg2[i_step] = gyration_radius_squared(pol)
        traj.magnetization[i_step] = mean(pol.spins)
        traj.energies[i_step] = energy

        #=if !isinteger(energy)
            println("Non integer energy:  ", energy)
            println(pol.coord)
            println(pol.neighbours .- pol.trial_neighbours)
            println(pol.trial_neighbours)
            break
        end=#
        empty!(pol.hash_saw)
        if i_step%sample_lag == 0
            i_sample = div(i_step,sample_lag)
            for i_mono in 1:pol.n_mono
                spins_configs[i_mono, i_sample]  = pol.spins[i_mono]
            end
        end
    end
end

#########################################################################
#########################################################################
#########################################################################
# Multiple Markov Chains

function MMC_run_base!(polymers::Array{Magnetic_polymer}, trajs::Array{MC_data},
                    n_strides::Int, stride::Int, inv_temps::Array{Float64})
    n_temps = length(inv_temps)
    accepted_swaps = 0
    temp_coord = zeros(Int,3)

    for i_temp in 1:n_temps
        MC_run!(polymers[i_temp], trajs[i_temp],1,stride)
    end

    for i_strides in 2:n_strides
        swap = rand(1:n_temps-1)
        delta_ene = (inv_temps[swap] - inv_temps[swap+1]) * (trajs[swap+1].energies[(i_strides-1)*stride] - trajs[swap].energies[(i_strides-1)*stride])
        alpha_swap = 0.0
        delta_ene<=0 ? alpha_swap=1.0 : alpha_swap=exp(-delta_ene)
        if alpha_swap >= rand()
            for i_mono in 1:polymers[1].n_mono
                for j in 1:3
                   temp_coord[j] = polymers[swap].coord[j,i_mono]
                   polymers[swap].coord[j,i_mono] = polymers[swap+1].coord[j,i_mono]
                   polymers[swap+1].coord[j,i_mono] =  temp_coord[j]
                end
                temp_spin = polymers[swap].spins[i_mono]
                polymers[swap].spins[i_mono] = polymers[swap+1].spins[i_mono]
                polymers[swap+1].spins[i_mono] = temp_spin
            end
            accepted_swaps += 1
        end

        for i_temp in 1:n_temps
            MC_run!(polymers[i_temp], trajs[i_temp],(i_strides-1)*stride+1,i_strides*stride)
        end
    end
    #println("Accepted_swaps: ", accepted_swaps)
end

function MMC_run_save!(polymers::Array{Magnetic_polymer}, trajs::Array{MC_data},
                    n_strides::Int, stride::Int, inv_temps::Array{Float64}, spins_configs::Matrix{Int}, sample_lag::Int, n_samples::Int)
    n_temps = length(inv_temps)
    accepted_swaps = 0
    temp_coord = zeros(Int,3)

    MC_run!(polymers[1], trajs[1],1,stride, spins_configs,sample_lag, n_samples)
    # Only the lowest temperature is the system we're using for our likelihood approx
    # the chains at higher temps are only used to "fluidify" the chain of interest
    for i_temp in 2:n_temps
        MC_run!(polymers[i_temp], trajs[i_temp],1,stride, sample_lag, n_samples)
    end

    for i_strides in 2:n_strides
        swap = rand(1:n_temps-1)
        delta_ene = (inv_temps[swap] - inv_temps[swap+1]) * (trajs[swap+1].energies[i_strides*stride] - trajs[swap].energies[i_strides*stride])
        alpha_swap = 0.0
        delta_ene<=0 ? alpha_swap=1.0 : alpha_swap=exp(-delta_ene)
        if alpha_swap >= rand()
            for i_mono in 1:polymers[1].n_mono
                for j in 1:3
                   temp_coord[j] = polymers[swap].coord[j,i_mono]
                   polymers[swap].coord[j,i_mono] = polymers[swap+1].coord[j,i_mono]
                   polymers[swap+1].coord[j,i_mono] =  temp_coord[j]
                end
                temp_spin = polymers[swap].spins[i_mono]
                polymers[swap].spins[i_mono] = polymers[swap+1].spins[i_mono]
                polymers[swap+1].spins[i_mono] = temp_spin
            end
            accepted_swaps += 1
        end
        MC_run!(polymers[1], trajs[1],(i_strides-1)*stride+1,i_strides*stride, spins_configs,sample_lag, n_samples)
        # Only the lowest temperature is the system we're using for our likelihood approx
        # the chains at higher temps are only used to "fluidify" the chain of interest
        for i_temp in 2:n_temps
            MC_run!(polymers[i_temp], trajs[i_temp],(i_strides-1)*stride+1,i_strides*stride, sample_lag, n_samples)
        end
    end
    #println("Accepted_swaps: ", accepted_swaps)
end


text:
# @testset "Resampling" begin
#     spec = mock_spectrum()
#     new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) ÷ 2.4))
#     res_spec = resample(spec, new_wave)

#     @test res_spec.wave == new_wave
#     @test length(res_spec.flux) == length(new_wave)

#     resample!(spec, new_wave)
#     @test res_spec.wave == spec.wave
#     @test res_spec.flux == spec.flux

#     # Unitful
#     spec = mock_spectrum(use_units=true)
#     new_wave = range(minimum(spec.wave), maximum(spec.wave), length=Integer(length(spec.wave) ÷ 2.4))
#     @assert unit(eltype(new_wave)) == unit(eltype(spec.wave))
#     res_spec = resample(spec, new_wave)

#     @test res_spec.wave == new_wave
#     @test length(res_spec.flux) == length(new_wave)

#     resample!(spec, new_wave)
#     @test res_spec.wave == spec.wave
#     @test res_spec.flux == spec.flux

#     # Test resampling to another Spectrum
#     spec1 = mock_spectrum(Integer(1e4))
#     spec2 = mock_spectrum(Integer(1e3))
#     res_spec = resample(spec1, spec2)
#     @test res_spec.wave == spec2.wave

#     resample!(spec1, spec2)
#     @test spec1.wave == spec2.wave
#     @test spec1.flux == res_spec.flux

#     # Unitful
#     spec1 = mock_spectrum(Integer(1e4), use_units=true)
#     spec2 = mock_spectrum(Integer(1e3), use_units=true)
#     res_spec = resample(spec1, spec2)

#     @test res_spec.wave == spec2.wave
#     @test unit(eltype(spec1.flux)) == unit(eltype(spec2.flux)) == unit(eltype(res_spec.flux))

#     # Test when spectra2 has different units
#     spec1 = mock_spectrum(Integer(1e4), use_units=true)
#     spec2 = mock_spectrum(Integer(1e3), use_units=true)
#     spec1.wave = uconvert.(u"cm", spec1.wave)
#     resample!(spec1, spec2)
#     @test ustrip.(spec1.wave) ≈ ustrip.(unit(eltype(spec1.wave)), spec2.wave)

#     # address bug when doing the same thing but affecting spec2
#     spec1 = mock_spectrum(Integer(1e4), use_units=true)
#     spec2 = mock_spectrum(Integer(1e3), use_units=true)
#     spec2.wave = uconvert.(u"cm", spec2.wave)
#     resample!(spec1, spec2)
#     @test ustrip.(spec1.wave) ≈ ustrip.(unit(eltype(spec1.wave)), spec2.wave)

# end


text:
using StatsBase, Distributions
using MDP

# Based on the description in Sec 8 of the following Technical Report
# > Shalabh Bhatnagar, Richard S. Sutton, Mohammad Ghavamzadeh, and Mark Lee
# > "Natural Actor–Critic Algorithms",
# > University of Alberta Department of Computing Science Technical Report TR09-10
# > June 2009.

# In the description of Garnet, there is a parameter τ that determines
# the non-stationarity of the mode. Here, we only implement the case when
# τ = 0. We will also assume that b << n, so we generate a sparse P matrix.

# The description of the garnet model says that the rewards should have normal
# distribution. We let the user to explicitly specify the reward distribution.

function garnet(n, m, b; rewards=Uniform(0,1), objective=:Max)
   P = SparseMatrixCSC[ spzeros(n,n) for a in 1:m ]
   S = 1:n
   for a in 1:m
       for i in 1:n
           nextState = sample(S, b; replace=false)
           nextProb  = rand(b)
           nextProb /= sum(nextProb)
           P[a][i, nextState] = nextProb
       end
   end
   c = rand(rewards, n,m)
   model = ProbModel(c, P; objective=objective)
end

function randomPolicy(model :: ProbModel; stochastic = true)
  (n, m) = (model.stateSize, model.actionSize)

  Q = Matrix[zero(model.P[1])]
  c = zeros(n,1)

  if !stochastic
    a = sample(1:m, n; replace=true)
    for i in 1:n
      Q[1][i,:] = model.P[a[i]][i,:]
      c[i,1] = model.cost[i, a[i]]
    end
  else
    for i in 1:n
      p  = rand(m)
      p /= sum(p)

      for a in 1:m
        Q[1][i,:] += model.P[a][i,:]*p[a]
        c[i,1]    += model.cost[i,a]*p[a]
      end
    end
  end

  return ProbModel(c, Q)
end


text:
@inline CI(a...) = CartesianIndex(a...)
@inline δ(a,d::Int) = CI(ntuple(i -> i==a ? 1 : 0, d))
@inline δ(a,I::CartesianIndex{N}) where {N} = δ(a,N)

@inline CR(a...) = CartesianIndices(a...)
@inline inside(M::NTuple{N,Int}) where {N} = CR(ntuple(i-> 2:M[i]-1,N))
@inline inside(a; reverse::Bool=false) =
        reverse ? Iterators.reverse(inside(size(a))) : inside(size(a))
function inside_u(N::NTuple{n,Int},j::Int)::CartesianIndices{n} where n
    CartesianIndices(ntuple( i-> i==j ? (3:N[i]-1) : (2:N[i]), n))
end
splitn(n) = Base.front(n),n[end]
size_u(u) = splitn(size(u))

import Base.mapreduce
@fastmath function mapreduce(f,op,R::CartesianIndices;init=0.)
    val = init
    @inbounds @simd for I ∈ R
        val = op(val,f(I))
    end
    val
end
L₂(a) = mapreduce(I->@inbounds(abs2(a[I])),+,inside(a))

macro inside(ex)
    @assert ex.head==:(=)
    a,I = Meta.parse.(split(string(ex.args[1]),union("[","]")))
    return quote
        @inbounds @simd for $I ∈ inside($a)
            $ex
        end
    end |> esc
end

function median(a,b,c)
    if a>b
        b>=c && return b
        a>c && return c
    else
        b<=c && return b
        a<c && return c
    end
    return a
end

"""
    apply!(f, c)

Apply a vector function `f(i,x)` to the faces of a uniform staggered array `c`.
"""
function apply!(f,c)
    N = size(c)
    for b ∈ 1:N[end]
        @simd for I ∈ CR(N[1:end-1])
            x = collect(Float16, I.I) # location at cell center
            x[b] -= 0.5               # location at face
            @inbounds c[I,b] = f(b,x) # apply function to location
        end
    end
end

"""
    slice(N,s,dims) -> R

Return `CartesianIndices` slicing through an array of size `N`.
"""
function slice(N::NTuple{n,Int},s::Int,dims::Int,low::Int=1)::CartesianIndices{n} where n
    CartesianIndices(ntuple( i-> i==dims ? (s:s) : (low:N[i]), n))
end

function BC!(a::AbstractArray{T,m},A,f=1) where {T,m}
    n = m-1
    N = ntuple(i -> size(a,i), n)
    for j ∈ 1:n, i ∈ 1:n
        if i==j # Inline direction
            for s ∈ (1,2,N[j]); @simd for I ∈ slice(N,s,j)
                a[I,i] = f*A[i] # Dirichlet
            end; end
        else    # Perpendicular directions
            @simd for I ∈ slice(N,1,j)
                a[I,i] = a[I+δ(j,I),i] # Neumann
            end
            @simd for I ∈ slice(N,N[j],j)
                a[I,i] = a[I-δ(j,I),i] # Neumann
            end
        end
    end
end
function BC!(a::AbstractArray{T,n}) where {T,n}
    N = size(a)
    for j ∈ 1:n
        @simd for I ∈ slice(N,1,j)
            a[I] = a[I+δ(j,I)] # Neumann
        end
        @simd for I ∈ slice(N,N[j],j)
            a[I] = a[I-δ(j,I)] # Neumann
        end
    end
end


text:
###Hastings & Powell 1991
### Reproducing figure 3

tspan = (0.0, 10000) # timespan
p = [5.0, 0.1, 3.0, 2.0, 0.4, 0.01] # parameter values

# Initial conditions serie 1
u0_serie1 = [0.77, 0.16, 9.9]   #initial conditions for x,y,z
prob_serie1 = ODEProblem(parameterized_model!, u0_serie1, tspan, p)
sol_serie1 = solve(prob_serie1)

plot(sol_serie1, vars = 1, xlim = (0, 500), ylim = (0, 1), c=:black, linewidth=2)

# Initial conditions changing x by 0.01
u0_serie2 = [0.78, 0.16, 9.9]

prob_serie2 = ODEProblem(parameterized_model!, u0_serie2, tspan, p)
sol_serie2 = solve(prob_serie2, )

fig3 = plot!(sol_serie2, vars = 1, xlim = (0, 500), ylim = (0, 1),
             c = :darkgrey, linestyle = :dash, linewidth = 2,
             xlabel = "time", ylabel = "x", guidefontsize = 10,
             grid = :none, legend = false, framestyle = :box)

## Export figures
savefig(fig3, joinpath("..", "article", "figures", "fig3.pdf"))