text:
module BPLDenoising

export generate_scalar_tv_cost, generate_cost_plot, validate_tv_parameter
export generate_2d_tv_cost, generate_2d_cost_plot
export scalar_bilevel_tv_learn, patch_bilevel_tv_learn
export scalar_bilevel_sumregs_learn, patch_bilevel_sumregs_learn, validate_sumregs_parameter

using Printf
using SparseArrays
using JLD2
using FileIO
using PGFPlots
#pushPGFPlotsOptions("scale=0.8")
using ColorTypes: Gray
import ColorVectorSpace
using ImageContrastAdjustment
using ImageQualityIndexes

using AlgTools.Util
using AlgTools.LinOps
using AlgTools.LinkedLists
using ImageTools.Visualise

using VariationalImaging
using VariationalImaging.GradientOps
using VariationalImaging.OpDenoise

include("BilevelVisualise.jl")
include("TRBox.jl")
#include("TVLearningFunctionOp.jl")
include("TVLearningFunctionVec.jl")
include("SumRegsLearningFunction.jl")
include("Datasets.jl")

using BPLDenoising.BilevelVisualise
using BPLDenoising.Datasets

const default_save_prefix = "output"

function TVDenoise(data,parameter::Real;visualize=false)
    denoise_params = (
        ρ = 0,
        α = parameter,
        op = FwdGradientOp(),
        # PDPS
        τ₀ = 5,
        σ₀ = 0.99/5,
        accel = true,
        save_results = false,
        maxiter = 10000,
        verbose_iter = 10001,
        save_iterations = false
    )
    st, iterate = initialise_visualisation(visualize)
    out = op_denoise_pdps(data; iterate=iterate, params=denoise_params)
    finalise_visualisation(st)
    return out
end

function TVDenoise(data,parameter::AbstractArray;visualize=false)
    p = PatchOp(parameter,data[:,:,1]) # Adjust parameter size
	x̄ = zeros(p.size_out)
	inplace!(x̄,p,parameter)
    denoise_params = (
        ρ = 0,
        α = x̄,
        op = FwdGradientOp(),
        # PDPS
        τ₀ = 5,
        σ₀ = 0.99/5,
        accel = true,
        save_results = false,
        maxiter = 10000,
        verbose_iter = 10001,
        save_iterations = false
    )
    st, iterate = initialise_visualisation(visualize)
    out = op_denoise_pdps(data; iterate=iterate, params=denoise_params)
    finalise_visualisation(st)
    return out
end

function L2CostFunction(u,true_)
    return 0.5*norm₂²(u-true_)
end

################################
# Scalar Cost Function Plotting
################################

function generate_cost(dataset_name, parameter_range, cost_function, denoise_function;freq=10,num_samples=1)
    true_,data = testdataset(dataset_name)
    true_ = Float64.(Gray{Float64}.(true_))[:,:,1:num_samples]
    data = Float64.(Gray{Float64}.(data))[:,:,1:num_samples]
    costs = zeros(size(parameter_range))
    iter = 1
    for i = 1:length(parameter_range)
        u = denoise_function(data,parameter_range[i])
        costs[i] = cost_function(u,true_)
        if iter%freq == 0
            @info "Denoising parameter $(parameter_range[i]): cost = $(costs[i])"
        end
        iter += 1
    end
    output_dir = joinpath(default_save_prefix, dataset_name)
    if isdir(output_dir) == false
        mkpath(output_dir)
    end
    @save joinpath(output_dir, dataset_name*"_cost.jld2") parameter_range costs
end

function generate_cost_plot(dataset_name)
    cost_path = joinpath(default_save_prefix,dataset_name)

    if isdir(cost_path) == false
        @error "No cost calculation found at $cost_path"
        return false
    end

    @load joinpath(cost_path,dataset_name*"_cost.jld2") parameter_range costs
    p = Axis(Plots.Linear(parameter_range,costs,mark="none"),style="grid=both", xlabel=L"$\alpha$", ylabel=L"$\|u-\bar{u}\|^2$", title="Scalar Cost", xmode="log", ymode="log")
    PGFPlots.save(joinpath(cost_path,dataset_name*"_cost_plot.tex"),p,include_preamble=false)
    PGFPlots.save(joinpath(cost_path,dataset_name*"_cost_plot.pdf"),p)
end

function generate_scalar_tv_cost(dataset_name, parameter_range; num_samples=1)
    return generate_cost(dataset_name,parameter_range,L2CostFunction,TVDenoise; num_samples)
end

################################
# 2D Cost Function Plotting
################################

function generate_2d_cost(dataset_name, parameter_range_1, parameter_range_2, cost_function, denoise_function;freq=10,num_samples=1)
    true_,data = testdataset(dataset_name)
    true_ = Float64.(Gray{Float64}.(true_))[:,:,1:num_samples]
    data = Float64.(Gray{Float64}.(data))[:,:,1:num_samples]
    costs = zeros(length(parameter_range_1),length(parameter_range_2))
    iter = 1
    for i = 1:length(parameter_range_1)
        for j = 1:length(parameter_range_2)
            α = [parameter_range_1[i];parameter_range_2[j]] .* ones(2,1)
            u = denoise_function(data,α)
            costs[i,j] = cost_function(u,true_)
            if iter%freq == 0
                @info "Denoising parameter $α: cost = $(costs[i,j])"
            end
            iter += 1
        end
    end
    output_dir = joinpath(default_save_prefix, dataset_name)
    if isdir(output_dir) == false
        mkpath(output_dir)
    end
    @save joinpath(output_dir, dataset_name*"_cost_2d.jld2") parameter_range_1 parameter_range_2 costs
end

function generate_2d_cost_plot(dataset_name)
    cost_path = joinpath(default_save_prefix,dataset_name)

    if isdir(cost_path) == false
        @error "No cost calculation found at $cost_path"
        return false
    end

    @load joinpath(cost_path,dataset_name*"_cost_2d.jld2") parameter_range_1 parameter_range_2 costs
    c = reshape(costs,length(parameter_range_1),length(parameter_range_2))
    p = Axis(Plots.Contour(c,parameter_range_1,parameter_range_2,style="dashed",levels=150.0:20.0:300.0),style="grid=both", xlabel=L"$\alpha_1$", ylabel=L"$\alpha_2$", title="2D Cost")
    #p = Axis(Plots.Image(costs,(0.005,0.03),(0.005,0.03)),style="grid=both", xlabel=L"$\alpha_1$", ylabel=L"$\alpha_2$", title="2D Cost"),levels=30.0:0.75:100.0
    PGFPlots.save(joinpath(cost_path,dataset_name*"_cost_plot_2d.tex"),p,include_preamble=false)
    PGFPlots.save(joinpath(cost_path,dataset_name*"_cost_plot_2d.pdf"),p)
end

function generate_2d_tv_cost(dataset_name, parameter_range_1, parameter_range_2;num_samples=1)
    return generate_2d_cost(dataset_name,parameter_range_1,parameter_range_2,L2CostFunction,TVDenoise;num_samples)
end

################################
# Save Experiments Results
################################

function save_results(params, b, b_data, x::Union{Real,AbstractVector{Float64}}, opt_img, st)
    if params.save_results
        out_path = joinpath(default_save_prefix,params.dataset_name)
        if isdir(out_path) == false
            mkpath(out_path)
        end
        perffile = joinpath(out_path,params.save_prefix * ".txt")
        qualityfile = joinpath(out_path,params.save_prefix * "_quality.txt")
        println("Saving " * perffile)
        write_log(perffile, st.log, "# params = $(params), x = $x\n")
        open(qualityfile,"w") do io
            write(io,"img_num \t orig_ssim \t orig_psnr \t out_ssim \t out_psnr\n")
            M,N,O = size(b)
            mean_ssim=0
            mean_psnr=0
            for i = 1:O
                noisy_ssim = assess_ssim(b[:,:,i],b_data[:,:,i])
                noisy_psnr = assess_psnr(b[:,:,i],b_data[:,:,i])
                out_ssim = assess_ssim(b[:,:,i],opt_img[:,:,i])
                out_psnr = assess_psnr(b[:,:,i],opt_img[:,:,i])
                write(io,"$i\t $noisy_ssim \t $noisy_psnr \t $out_ssim \t $out_psnr\n")
                mean_ssim += out_ssim
                mean_psnr += out_psnr

                fn = (t, ext, i) -> "$(joinpath(out_path,params.save_prefix))_$(t)_$(i).$(ext)"
                FileIO.save(File(format"PNG", fn("true", "png", i)), grayimg(b[:,:,i]))
                FileIO.save(File(format"PNG", fn("data", "png", i)), grayimg(b_data[:,:,i]))
                FileIO.save(File(format"PNG", fn("reco", "png", i)), grayimg(opt_img[:,:,i]))
            end
            write(io,"\t\t\t\t\t $(mean_ssim/O)\t $(mean_psnr/O)\n")
        end
    end
end

function save_results(params, b, b_data, x::AbstractArray{T,2}, opt_img, st) where T
    if params.save_results
        out_path = joinpath(default_save_prefix,params.dataset_name)
        if isdir(out_path) == false
            mkpath(out_path)
        end
        perffile = joinpath(out_path,params.save_prefix * ".txt")
        qualityfile = joinpath(out_path,params.save_prefix * "_quality.txt")
        println("Saving " * perffile)
        write_log(perffile, st.log, "# params = $(params), x = $x\n")
        open(qualityfile,"w") do io
            write(io,"img_num \t orig_ssim \t orig_psnr \t out_ssim \t out_psnr\n")
            M,N,O = size(b)
            mean_ssim=0
            mean_psnr=0
            for i = 1:O
                noisy_ssim = assess_ssim(b[:,:,i],b_data[:,:,i])
                noisy_psnr = assess_psnr(b[:,:,i],b_data[:,:,i])
                out_ssim = assess_ssim(b[:,:,i],opt_img[:,:,i])
                out_psnr = assess_psnr(b[:,:,i],opt_img[:,:,i])
                write(io,"$i\t $noisy_ssim \t $noisy_psnr \t $out_ssim \t $out_psnr\n")
                mean_ssim += out_ssim
                mean_psnr += out_psnr

                fn = (t, ext, i) -> "$(joinpath(out_path,params.save_prefix))_$(t)_$(i).$(ext)"
                FileIO.save(File(format"PNG", fn("true", "png", i)), grayimg(b[:,:,i]))
                FileIO.save(File(format"PNG", fn("data", "png", i)), grayimg(b_data[:,:,i]))
                FileIO.save(File(format"PNG", fn("reco", "png", i)), grayimg(opt_img[:,:,i]))
            end
            write(io,"\t\t\t\t\t $(mean_ssim/O)\t $(mean_psnr/O)\n")
        end
        p = PatchOp(x,b[:,:,1]) # Adjust parameter size
        x̄ = zeros(p.size_out)
        inplace!(x̄,p,x)
        fn_par = (t, ext) -> "$(joinpath(out_path,params.save_prefix))_$(t).$(ext)"
        adjust_histogram!(x̄,LinearStretching())
        FileIO.save(File(format"PNG", fn_par("par", "png")), grayimg(x̄))
    end
end

function save_results(params, b, b_data, x::AbstractArray{T,3}, opt_img, st) where T
    if params.save_results
        out_path = joinpath(default_save_prefix,params.dataset_name)
        if isdir(out_path) == false
            mkpath(out_path)
        end
        perffile = joinpath(out_path,params.save_prefix * ".txt")
        qualityfile = joinpath(out_path,params.save_prefix * "_quality.txt")
        println("Saving " * perffile)
        write_log(perffile, st.log, "# params = $(params), x = $x\n")
        open(qualityfile,"w") do io
            write(io,"img_num \t orig_ssim \t orig_psnr \t out_ssim \t out_psnr\n")
            M,N,O = size(b)
            mean_ssim=0
            mean_psnr=0
            for i = 1:O
                noisy_ssim = assess_ssim(b[:,:,i],b_data[:,:,i])
                noisy_psnr = assess_psnr(b[:,:,i],b_data[:,:,i])
                out_ssim = assess_ssim(b[:,:,i],opt_img[:,:,i])
                out_psnr = assess_psnr(b[:,:,i],opt_img[:,:,i])
                write(io,"$i\t $noisy_ssim \t $noisy_psnr \t $out_ssim \t $out_psnr\n")
                mean_ssim += out_ssim
                mean_psnr += mean_psnr

                fn = (t, ext, i) -> "$(joinpath(out_path,params.save_prefix))_$(t)_$(i).$(ext)"
                FileIO.save(File(format"PNG", fn("true", "png", i)), grayimg(b[:,:,i]))
                FileIO.save(File(format"PNG", fn("data", "png", i)), grayimg(b_data[:,:,i]))
                FileIO.save(File(format"PNG", fn("reco", "png", i)), grayimg(opt_img[:,:,i]))
            end
            write(io,"\t\t\t\t\t $(mean_ssim/O)\t $(mean_psnr/O)\n")
        end
        p = PatchOp(x,b[:,:,1]) # Adjust parameter size
        x̄ = p(x)
        fn_par = (t, ext, i) -> "$(joinpath(out_path,params.save_prefix))_$(t)_$(i).$(ext)"
        adjust_histogram!(x̄,LinearStretching())
        FileIO.save(File(format"PNG", fn_par("par", "png", 1)), grayimg(x̄[:,:,1]))
        FileIO.save(File(format"PNG", fn_par("par", "png", 2)), grayimg(x̄[:,:,2]))
        FileIO.save(File(format"PNG", fn_par("par", "png", 3)), grayimg(x̄[:,:,3]))
    end
end

###########################
# Scalar Bilevel Experiment
###########################

const default_params = (
    verbose_iter = 1,
    maxiter = 20,
    save_results = true,
    dataset_name = "cameraman_128_5",
    save_iterations = false,
    tol = 1e-5,
    num_samples = 1
)

const bilevel_params = (
    η₁ = 0.25,
    η₂ = 0.75,
    β₁ = 0.25,
    β₂ = 1.9,
    Δ₀ = 0.1,
    α₀ = 0.1
)

function scalar_bilevel_tv_learn(;visualise=true, save_prefix=default_save_prefix, kwargs...)
    # Parameters for this experiment
    params = default_params ⬿ bilevel_params ⬿ kwargs
    params = params ⬿ (save_prefix = "tv_optimal_parameter_scalar_" * params.dataset_name,)
    # Load dataset
    b,b_noisy = testdataset(params.dataset_name)
    b = Float64.(Gray{Float64}.(b))[:,:,1:params.num_samples]
    b_noisy = Float64.(Gray{Float64}.(b_noisy))[:,:,1:params.num_samples]
    # Launch (background) visualiser
    st, iterate = initialise_bilevel_visualisation(visualise)
    # Run algorithm
    x, u, st = bilevel_learn((b,b_noisy),tv_op_learning_function; xinit=params.α₀,iterate=iterate, params=params)
    adjust_histogram!(u,LinearStretching())
    adjust_histogram!(b,LinearStretching())
    adjust_histogram!(b_noisy,LinearStretching())
    # Save results
    save_results(params, b, b_noisy, x, u, st)
    # Exit background visualiser
    finalise_visualisation(st)
end

###########################
# Patch Bilevel Experiment
###########################

const patch_bilevel_params = (
    η₁ = 0.25,
    η₂ = 0.75,
    β₁ = 0.25,
    β₂ = 1.9,
    Δ₀ = 0.0001,
    α₀ = 0.0001*ones(2,2)
)

function patch_bilevel_tv_learn(;visualise=true, save_prefix=default_save_prefix, kwargs...)
    # Parameters for this experiment
    params = default_params ⬿ patch_bilevel_params ⬿ kwargs
    params = params ⬿ (save_prefix = "tv_optimal_parameter_$(size(params.α₀))_" * params.dataset_name,)
    # Load dataset
    b,b_noisy = Datasets.testdataset(params.dataset_name)
    b = Float64.(Gray{Float64}.(b))[:,:,1:params.num_samples]
    b_noisy = Float64.(Gray{Float64}.(b_noisy))[:,:,1:params.num_samples]
    # Launch (background) visualiser
    st, iterate = initialise_bilevel_visualisation(visualise)
    # Run algorithm
    x, u, st = bilevel_learn((b,b_noisy),tv_op_learning_function; xinit=params.α₀,iterate=iterate, params=params)
    adjust_histogram!(u,LinearStretching())
    # Save results
    save_results(params, b, b_noisy, x, u, st)
    # Exit background visualiser
    finalise_visualisation(st)
end

###########################
# Validate learned tv parameter
###########################
function validate_tv_parameter(parameter; kwargs...)
    params = default_params ⬿ bilevel_params ⬿ kwargs
    params = params ⬿ (save_prefix = "val_tv_optimal_parameter_scalar_$(size(parameter))_" * params.dataset_name,)
    println(params)
    img,noisy = testdataset(params.dataset_name)
    u = TVDenoise(noisy,parameter)
    cost = L2CostFunction(u,img)
    @info "Denoising parameter $(parameter): cost = $(cost)"
    out_path = joinpath(default_save_prefix,params.dataset_name)
    if isdir(out_path) == false
        mkpath(out_path)
    end
    qualityfile = joinpath(out_path,params.save_prefix * "_quality.txt")
    open(qualityfile,"w") do io
        write(io,"img_num \t orig_ssim \t orig_psnr \t out_ssim \t out_psnr\n")
        M,N,O = size(img)
        ssim_mean = 0
        psnr_mean = 0
        for i = 1:O
            noisy_ssim = assess_ssim(img[:,:,i],noisy[:,:,i])
            noisy_psnr = assess_psnr(img[:,:,i],noisy[:,:,i])
            out_ssim = assess_ssim(img[:,:,i],u[:,:,i])
            out_psnr = assess_psnr(img[:,:,i],u[:,:,i])
            write(io,"$i\t $noisy_ssim \t $noisy_psnr \t $out_ssim \t $out_psnr\n")
            ssim_mean += out_ssim
            psnr_mean += out_psnr

            fn = (t, ext, i) -> "$(joinpath(out_path,params.save_prefix))_$(t)_$(i).$(ext)"
            FileIO.save(File(format"PNG", fn("true", "png", i)), grayimg(img[:,:,i]))
            FileIO.save(File(format"PNG", fn("data", "png", i)), grayimg(noisy[:,:,i]))
            FileIO.save(File(format"PNG", fn("reco", "png", i)), grayimg(u[:,:,i]))
        end
        write(io,"\t\t\t\t\t $(ssim_mean/O)\t $(psnr_mean/O)\n")
    end
end

###########################
# Scalar Sumregs Bilevel Experiment
###########################

const sumregs_bilevel_params = (
    η₁ = 0.25,
    η₂ = 0.75,
    β₁ = 0.25,
    β₂ = 1.9,
    Δ₀ = 0.01,
    α₀ = [0.001;0.001;0.001]
)

function scalar_bilevel_sumregs_learn(;visualise=true, save_prefix=default_save_prefix, kwargs...)
    # Parameters for this experiment
    params = default_params ⬿ sumregs_bilevel_params ⬿ kwargs
    params = params ⬿ (save_prefix = "sumregs_optimal_parameter_scalar_" * params.dataset_name,)
    # Load dataset
    b,b_noisy = testdataset(params.dataset_name)
    b = Float64.(Gray{Float64}.(b))[:,:,1:params.num_samples]
    b_noisy = Float64.(Gray{Float64}.(b_noisy))[:,:,1:params.num_samples]
    # Launch (background) visualiser
    st, iterate = initialise_bilevel_visualisation(visualise)
    # Run algorithm
    x, u, st = bilevel_learn((b,b_noisy),sumregs_learning_function; xinit=params.α₀,iterate=iterate, params=params)
    adjust_histogram!(u,LinearStretching())
    # Save results
    save_results(params, b, b_noisy, x, u, st)
    # Exit background visualiser
    finalise_visualisation(st)
end

###########################
# Patch Sumregs Bilevel Experiment
###########################

const patch_sumregs_bilevel_params = (
    η₁ = 0.25,
    η₂ = 0.75,
    β₁ = 0.25,
    β₂ = 1.5,
    Δ₀ = 0.1,
    α₀ = 0.001*ones(2,2,3)
)

function patch_bilevel_sumregs_learn(;visualise=true, save_prefix=default_save_prefix, kwargs...)
    # Parameters for this experiment
    params = default_params ⬿ patch_sumregs_bilevel_params ⬿ kwargs
    params = params ⬿ (save_prefix = "sumregs_optimal_parameter_patch_$(size(params.α₀))" * params.dataset_name,)
    # Load dataset
    b,b_noisy = testdataset(params.dataset_name)
    b = Float64.(Gray{Float64}.(b))[:,:,1:params.num_samples]
    b_noisy = Float64.(Gray{Float64}.(b_noisy))[:,:,1:params.num_samples]
    # Launch (background) visualiser
    st, iterate = initialise_bilevel_visualisation(visualise)
    # Run algorithm
    x, u, st = bilevel_learn((b,b_noisy),sumregs_learning_function; xinit=params.α₀,iterate=iterate, params=params)
    adjust_histogram!(u,LinearStretching())
    # Save results
    save_results(params, b, b_noisy, x, u, st)
    # Exit background visualiser
    finalise_visualisation(st)
end

function patch_bilevel_sumregs_learn(image_pair::AbstractArray{T,3},dataset_name;visualise=true, save_prefix=default_save_prefix, kwargs...) where T
    # Parameters for this experiment
    params = default_params ⬿ patch_sumregs_bilevel_params ⬿ kwargs
    params = params ⬿ (save_prefix = "sumregs_optimal_parameter_patch_$(size(params.α₀))" * dataset_name,)
    # Load dataset
    b = image_pair[:,:,1]
    b_noisy = image_pair[:,:,2]
    b = Float64.(Gray{Float64}.(b))[:,:,1:params.num_samples]
    b_noisy = Float64.(Gray{Float64}.(b_noisy))[:,:,1:params.num_samples]
    # Launch (background) visualiser
    st, iterate = initialise_bilevel_visualisation(visualise)
    # Run algorithm
    x, u, st = bilevel_learn((b,b_noisy),sumregs_learning_function; xinit=params.α₀,iterate=iterate, params=params)
    adjust_histogram!(u,LinearStretching())
    # Save results
    save_results(params, b, b_noisy, x, u, st)
    # Exit background visualiser
    finalise_visualisation(st)
end

###########################
# Validate learned tv parameter
###########################
function validate_sumregs_parameter(parameter; kwargs...)
    params = default_params ⬿ bilevel_params ⬿ kwargs
    params = params ⬿ (save_prefix = "val_sumregs_optimal_parameter_scalar_$(size(parameter))_" * params.dataset_name,)
    println(params)
    img,noisy = testdataset(params.dataset_name)
    u,cost,grad = sumregs_learning_function(parameter,noisy,0.1)
    @info "Denoising parameter $(parameter): cost = $(cost)"
    out_path = joinpath(default_save_prefix,params.dataset_name)
    if isdir(out_path) == false
        mkpath(out_path)
    end
    qualityfile = joinpath(out_path,params.save_prefix * "_quality.txt")
    open(qualityfile,"w") do io
        write(io,"img_num \t orig_ssim \t orig_psnr \t out_ssim \t out_psnr\n")
        M,N,O = size(img)
        ssim_mean = 0
        psnr_mean = 0
        for i = 1:O
            noisy_ssim = assess_ssim(img[:,:,i],noisy[:,:,i])
            noisy_psnr = assess_psnr(img[:,:,i],noisy[:,:,i])
            out_ssim = assess_ssim(img[:,:,i],u[:,:,i])
            out_psnr = assess_psnr(img[:,:,i],u[:,:,i])
            write(io,"$i\t $noisy_ssim \t $noisy_psnr \t $out_ssim \t $out_psnr\n")
            ssim_mean += out_ssim
            psnr_mean += out_psnr

            fn = (t, ext, i) -> "$(joinpath(out_path,params.save_prefix))_$(t)_$(i).$(ext)"
            FileIO.save(File(format"PNG", fn("true", "png", i)), grayimg(img[:,:,i]))
            FileIO.save(File(format"PNG", fn("data", "png", i)), grayimg(noisy[:,:,i]))
            FileIO.save(File(format"PNG", fn("reco", "png", i)), grayimg(u[:,:,i]))
        end
        write(io,"\t\t\t\t\t $(ssim_mean/O)\t $(psnr_mean/O)\n")
    end
end

end # Module

max_stars_repo_name: dvillacis/BPLDenoising

max_issues_repo_name: dvillacis/BPLDenoising

max_forks_repo_name: dvillacis/BPLDenoising


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "TetGen"
version = v"1.5"

# Collection of sources required to build tetgen
sources = [
    "https://github.com/ufz/tetgen.git" =>
    "9c4b70d479c6f67cb9e16dbe65f81ec3b821acad",
    "https://github.com/JuliaGeometry/TetgenBuilder.git" =>
    "f3fd2796459b3b4217ff4c240603e9193940b7ea"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd tetgen/
${CXX} -c -fPIC -std=c++11 -O3 predicates.cxx -o predicates.o
${CXX} -c -fPIC -std=c++11 -O3 -DTETLIBRARY -c tetgen.cxx -o tetgen.o
libdir="lib"
if [[ ${target} == *-mingw32 ]]; then     libdir="bin"; else     libdir="lib"; fi
mkdir ${prefix}/${libdir}
cp tetgen.h ../TetgenBuilder/tetgen.h
${CXX} $LDFLAGS -shared -fPIC tetgen.o predicates.o -DTETLIBRARY -o ${prefix}/${libdir}/libtet.${dlext} ../TetgenBuilder/cwrapper.cpp

exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = BinaryBuilder.supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libtet", :libtet)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl


text:
#   This file is part of AbstractTensors.jl. It is licensed under the GPL license
#   Grassmann Copyright (C) 2019 Michael Reed

module AbstractTensors

# universal root Tensor type

abstract type TensorAlgebra{V} <: Number end

# V, VectorSpace produced by DirectSum

import DirectSum: vectorspace, value, dual
import LinearAlgebra: dot, cross, norm, UniformScaling, I

# parameters accessible from anywhere

Base.@pure vectorspace(::T) where T<:TensorAlgebra{V} where V = V
Base.@pure Base.ndims(::T) where T<:TensorAlgebra{V} where V = ndims(V)

# universal vector space interopability

@inline interop(op::Function,a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{V}} where V = op(a,b)

# ^^ identity ^^ | vv union vv #

@inline function interop(op::Function,a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{W}} where {V,W}
    VW = V ∪ W
    return op(VW(a),VW(b))
end

# abstract tensor form evaluation

@inline interform(a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{V}} where V = a(b)
@inline function interform(a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{W}} where {V,W}
    VW = V ∪ W
    return VW(a)(VW(b))
end

# extended compatibility interface

export TensorAlgebra, interop, interform, scalar, involute, norm, unit, even, odd
export ⊖, ⊗, ⊛, ⊙, ⊠, ⨼, ⨽, ⋆, ∗, ⁻¹, ǂ, ₊, ₋, ˣ

# some shared presets

for op ∈ (:(Base.:+),:(Base.:-),:(Base.:*),:⊗,:⊛,:∗,:⨼,:⨽,:dot,:cross,:contraction,:(Base.:|),:(Base.:(==)),:(Base.:<),:(Base.:>),:(Base.:<<),:(Base.:>>),:(Base.:>>>),:(Base.div),:(Base.rem),:(Base.:&))
    @eval begin
        @inline $op(a::A,b::B) where {A<:TensorAlgebra,B<:TensorAlgebra} = interop($op,a,b)
        @inline $op(a::A,b::UniformScaling) where A<:TensorAlgebra{V} where V = $op(a,V(b))
        @inline $op(a::UniformScaling,b::B) where B<:TensorAlgebra{V} where V = $op(V(a),b)
    end
end

const ⊖ = *
@inline ⋆(t::UniformScaling{T}) where T = T<:Bool ? (t.λ ? 1 : -1) : t.λ
@inline ⊛(a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{V}} where V = scalar(contraction(a,b))
@inline ⨽(a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{V}} where V = contraction(a,b)
@inline ⨼(a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{V}} where V = contraction(b,a)
@inline Base.:<(a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{V}} where V = contraction(b,a)
@inline Base.:>(a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{V}} where V = contraction(a,b)
@inline Base.:|(a::A,b::B) where {A<:TensorAlgebra{V},B<:TensorAlgebra{V}} where V = contraction(a,b)

for op ∈ (:(Base.:+),:(Base.:*))
    @eval $op(t::T) where T<:TensorAlgebra = t
end
for op ∈ (:|,:!), T ∈ (TensorAlgebra,UniformScaling)
    @eval Base.$op(t::$T) = ⋆(t)
end
for op ∈ (:⊙,:⊠)
    @eval function $op end
end
for op ∈ (:scalar,:involute,:even)
    @eval $op(t::T) where T<:Real = t
end
odd(::T) where T<:Real = 0

@inline Base.exp(t::T) where T<:TensorAlgebra = 1+expm1(t)
@inline Base.log(b,t::T) where T<:TensorAlgebra = log(t)/log(b)
@inline Base.:^(b::S,t::T) where {S<:Number,T<:TensorAlgebra} = exp(t*log(b))
@inline Base.:^(a::A,b::UniformScaling) where A<:TensorAlgebra{V} where V = ^(a,V(b))
@inline Base.:^(a::UniformScaling,b::B) where B<:TensorAlgebra{V} where V = ^(V(a),b)

for base ∈ (2,10)
    fl,fe = (Symbol(:log,base),Symbol(:exp,base))
    @eval Base.$fl(t::T) where T<:TensorAlgebra = $fl(ℯ)*log(t)
    @eval Base.$fe(t::T) where T<:TensorAlgebra = exp(log($base)*t)
end

@inline Base.cos(t::T) where T<:TensorAlgebra{V} where V = cosh(V(I)*t)
@inline Base.sin(t::T) where T<:TensorAlgebra{V} where V = (i=V(I);sinh(i*t)/i)
@inline Base.tan(t::T) where T<:TensorAlgebra = sin(t)/cos(t)
@inline Base.cot(t::T) where T<:TensorAlgebra = cos(t)/sin(t)
@inline Base.sec(t::T) where T<:TensorAlgebra = inv(cos(t))
@inline Base.csc(t::T) where T<:TensorAlgebra = inv(sin(t))
@inline Base.asec(t::T) where T<:TensorAlgebra = acos(inv(t))
@inline Base.acsc(t::T) where T<:TensorAlgebra = asin(inv(t))
@inline Base.sech(t::T) where T<:TensorAlgebra = inv(cosh(t))
@inline Base.csch(t::T) where T<:TensorAlgebra = inv(sinh(t))
@inline Base.asech(t::T) where T<:TensorAlgebra = acosh(inv(t))
@inline Base.acsch(t::T) where T<:TensorAlgebra = asinh(inv(t))
@inline Base.tanh(t::T) where T<:TensorAlgebra = sinh(t)/cosh(t)
@inline Base.coth(t::T) where T<:TensorAlgebra = cosh(t)/sinh(t)
@inline Base.asinh(t::T) where T<:TensorAlgebra = log(t+sqrt(1+t^2))
@inline Base.acosh(t::T) where T<:TensorAlgebra = log(t+sqrt(t^2-1))
@inline Base.atanh(t::T) where T<:TensorAlgebra = (log(1+t)-log(1-t))/2
@inline Base.acoth(t::T) where T<:TensorAlgebra = (log(t+1)-log(t-1))/2
Base.asin(t::T) where T<:TensorAlgebra{V} where V =(i=V(I);-i*log(i*t+sqrt(1-t^2)))
Base.acos(t::T) where T<:TensorAlgebra{V} where V =(i=V(I);-i*log(t+i*sqrt(1-t^2)))
Base.atan(t::T) where T<:TensorAlgebra{V} where V =(i=V(I);(-i/2)*(log(1+i*t)-log(1-i*t)))
Base.acot(t::T) where T<:TensorAlgebra{V} where V =(i=V(I);(-i/2)*(log(t-i)-log(t+i)))
Base.sinc(t::T) where T<:TensorAlgebra{V} where V = iszero(t) ? one(V) : (x=(1*π)*t;sin(x)/x)
Base.cosc(t::T) where T<:TensorAlgebra{V} where V = iszero(t) ? zero(V) : (x=(1*π)*t; cos(x)/t - sin(x)/(x*t))

# absolute value norm

@inline Base.abs(t::T) where T<:TensorAlgebra{V} where V = sqrt(abs2(t))
@inline Base.abs2(t::T) where T<:TensorAlgebra = t∗t
@inline norm(t::T) where T<:TensorAlgebra = norm(value(t))
@inline unit(t::T) where T<:TensorAlgebra = t/abs(t)
@inline Base.iszero(t::T) where T<:TensorAlgebra = norm(t) ≈ 0
@inline Base.isone(t::T) where T<:TensorAlgebra = norm(t) ≈ value(scalar(t)) ≈ 1

# identity elements

for id ∈ (:zero,:one)
    @eval begin
        @inline Base.$id(::T) where T<:TensorAlgebra{V} where V = zero(V)
        @inline Base.$id(::Type{T}) where T<:TensorAlgebra{V} where V = zero(V)
    end
end

# postfix operators

struct Postfix{Op} end
@inline Base.:*(t,op::Postfix) = op(t)
for op ∈ (:⁻¹,:ǂ,:₊,:₋,:ˣ)
    @eval const $op = $(Postfix{op}())
end
@inline (::Postfix{:⁻¹})(t) = inv(t)
@inline (::Postfix{:ǂ})(t) = conj(t)
@inline (::Postfix{:₊})(t) = even(t)
@inline (::Postfix{:₋})(t) = odd(t)
@inline (::Postfix{:ˣ})(t) = involute(t)

end # module


text:
mutable struct History{I,V} <: UnivalueHistory{I}
    lastiter::I
    iterations::Vector{I}
    values::Vector{V}

    function History(::Type{V}, ::Type{I} = Int) where {I,V}
        new{I,V}(typemin(I), Array{I}(undef, 0), Array{V}(undef, 0))
    end
end

Base.length(history::History) = length(history.iterations)
Base.enumerate(history::History) = zip(history.iterations, history.values)
Base.first(history::History) = history.iterations[1], history.values[1]
Base.last(history::History) = history.iterations[end], history.values[end]
Base.get(history::History) = history.iterations, history.values

function Base.push!(
        history::History{I,V},
        iteration::I,
        value::V) where {I,V}
    lastiter = history.lastiter
    iteration > lastiter || throw(ArgumentError("Iterations must increase over time"))
    history.lastiter = iteration
    push!(history.iterations, iteration)
    push!(history.values, value)
    value
end

function Base.push!(
        history::History{I,V},
        value::V) where {I,V}
    lastiter = history.lastiter == typemin(I) ? zero(I) : history.lastiter
    iteration = lastiter + one(history.lastiter)
    history.lastiter = iteration
    push!(history.iterations, iteration)
    push!(history.values, value)
    value
end

Base.print(io::IO, history::History{I,V}) where {I,V} = print(io, "$(length(history)) elements {$I,$V}")

function Base.show(io::IO, history::History{I,V}) where {I,V}
    println(io, "History")
    println(io, "  * types: $I, $V")
    print(io,   "  * length: $(length(history))")
end

function Base.resize!(history::History, n)
    iter, vals = get(history)
    resize!(iter, n)
    resize!(vals, n)
    history.lastiter = last(iter)

    history
end

"""
    increment!(trace, iter, val)

Increments the value for a given iteration if it exists, otherwise adds the iteration with an ordinary push.
"""
function increment!(trace::History{I,V}, iter::Number, val)  where {I,V}
    if !isempty(trace.iterations)
        if trace.lastiter == iter # Check most common case to make it faster
            i = length(trace.iterations)
        else
            i = findfirst(isequal(iter), trace.iterations)
        end
        if i != nothing
            return (trace.values[i] += val)
        end
    end
    push!(trace, iter, val)
end


text:
# This file contains types and structures for economic DEA model
"""
    AbstractEconomicDEAlModel
An abstract type representing an economic DEA model.
"""
abstract type AbstractEconomicDEAModel <: AbstractTechnicalDEAModel end

"""
    AbstractCostDEAModel
An abstract type representing a cost DEA model.
"""
abstract type AbstractCostDEAModel <: AbstractEconomicDEAModel end

"""
    AbstractRevenueDEAModel
An abstract type representing a revenue DEA model.
"""
abstract type AbstractRevenueDEAModel <: AbstractEconomicDEAModel end

"""
    AbstractProfitDEAModel
An abstract type representing a revenue DEA model.
"""
abstract type AbstractProfitDEAModel <: AbstractEconomicDEAModel end

"""
    AbstractProfitabilityDEAModel
An abstract type representing a revenue DEA model.
"""
abstract type AbstractProfitabilityDEAModel <: AbstractEconomicDEAModel end

"""
    efficiency(model::AbstractEconomicDEAModel)
Return efficiency scores of an economic DEA model.

# Optional Arguments
- `type=Economic`: type of efficiency scores to return.

Type specification:
- `:Economic`: returns economic efficiency of the model.
- `:CRS`: returns technical efficiency under constant returns to scale.
- `:VRS`: returns technical efficiency under variable returns to scale.
- `:Scale`: returns scale efficiency.
- `:Allocative`: returns allocative efficiency.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];
julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];
julia> W = [2 1; 2 1; 2 1; 2 1; 2 1.0];
julia> P = [3 2; 3 2; 3 2; 3 2; 3 2.0];
julia> profitbl = deaprofitability(X, Y, W, P)
julia> efficiency(profitbl)
5-element Array{Float64,1}:
 0.38795983677810825
 0.9999999082180037
 0.7652173234322985
 0.249999984621284
 0.15879016408241559

julia> efficiency(profitbl, :Allocative)
5-element Array{Float64,1}:
 0.6096511887087087
 0.9999999387832732
 0.765217344390068
 0.9999999256504388
 0.6086956149040872
```
"""
function efficiency(model::AbstractEconomicDEAModel, type::Symbol = :Economic)::Vector

    if type == :Economic
        return model.eff
    end

    if type == :Technical
        if isdefined(model, :techeff)
            return model.techeff
        end
    end

    if type == :CRS
        if isdefined(model, :crseff)
            return model.crseff
        end
    end

    if type == :VRS
        if isdefined(model, :vrseff)
            return model.vrseff
        end
    end

    if type == :Scale
        if isdefined(model, :scaleff)
            return model.scaleff
        end
    end

    if type == :Allocative
        if isdefined(model, :alloceff)
            return model.alloceff
        end
    end

    error(typeof(model), " has no efficiency type ", string(type))

end


text:
# Run tests locally
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
# Pkg.test("AlgamesDriving")

Pkg.activate(joinpath(@__DIR__, "../test"))
Pkg.activate(joinpath(@__DIR__, "../docs"))


text:
"""
    CosineKernel

The cosine kernel is a stationary kernel for a sinusoidal given by
```
    κ(x,y) = cos( π * (x-y) )
```

"""
struct CosineKernel <: BaseKernel end

kappa(κ::CosineKernel, d::Real) = cospi(d)

metric(::CosineKernel) = Euclidean()


text:
using ToyFHE
using OffsetArrays
using Test
using Primes
using GaloisFields
using LinearAlgebra
using Random

N = 2^5
q₀ = nextprime(Int128(2)^40 + 1, 1; interval=2N)
q₁ = nextprime(q₀ + 2N, 1; interval=2N)
ps = nextprime(q₁ + 2N, 1; interval=2N)

ℛ = let CT = ToyFHE.CRTEncoded{3, Tuple{GaloisField.((q₀,q₁,ps))...}}
    ζ₂n = GaloisFields.minimal_primitive_root(CT, 2N)
    ToyFHE.NegacyclicRing{CT, N}(ζ₂n)
end

ckks_params = ModulusRaised(CKKSParams(ℛ, 0, 3.2))
kp = keygen(ckks_params)

scale = 2^40
Tscale = FixedRational{scale}

plain = CKKSEncoding{Tscale}(zero(ℛ_cipher(ckks_params)))
plain .= OffsetArray(1:div(N,2), 0:div(N,2)-1)

c = encrypt(kp, plain)
@test decrypt(kp, keyswitch(ToyFHE.make_eval_key(Random.GLOBAL_RNG, kp.priv.secret=>kp.priv), c)) ≈ plain atol=10^-8


text:
filepath = joinpath(Pkg.dir("Survival", "examples"), "rossi.csv")
rossi = readtable(filepath)
rossi[:event] = Event.(rossi[:week],rossi[:arrest] .== 0)

outcome = coxph( @formula(event ~ fin+age+race+wexp+mar+paro+prio) ,rossi ; tol = 1e-8)

filepath_coefs = joinpath(Pkg.dir("Survival", "test"), "expected_coefmat.jld")
expected_coefmat = JLD.load(filepath_coefs, "expected_coefmat")
@test outcome.coefmat.cols[1:3] ≈ expected_coefmat.cols[1:3] atol = 1e-6
println("Cox regression is fine")
println("$(outcome.coefmat) \n ≈ \n $expected_coefmat")


text:
"""
    Mode(k0, frequency, data, weighting, basis, label)

Eigenmode of a photonic crystal expressed on a plane-wave `basis` with a weighted inner product.
"""
struct Mode
    k0::Vector{Float64}
    frequency::ComplexF64
    data::Vector{ComplexF64}
    weighting::Matrix{ComplexF64}
    basis::PlaneWaveBasis
    label::String
end

"""
    HilbertSpace(k0, data, weighting, basis)

Hilbert space spanned by the eigenvectors in each column of `data`.

The eigenvectors will be orthonormalised using Gram-Schmidt orthonormalisation,
see [`orthonormalise`](@ref).
"""
struct HilbertSpace
    k0::Vector{Float64}
    data::Matrix{ComplexF64}
    weighting::Matrix{ComplexF64}
    basis::PlaneWaveBasis
    function HilbertSpace(k0::Vector{Float64}, data::Matrix{ComplexF64},
                            weighting::Matrix{ComplexF64}, basis::PlaneWaveBasis)
        # Inner constructor guarantees data will be orthonormalised
        data = orthonormalise(data, weighting=weighting)
        return new(k0, data, weighting, basis)
    end
end

"""
    HilbertSpace(modes::Array{Mode,1})

Returns the Hilbert space spanned by the `modes`.

The `data` of the Hilbert space is guaranteed to be orthonormal
under the weighting of the `modes`.
"""
function HilbertSpace(modes::Array{Mode,1})
    k0 = modes[1].k0
    data = zeros(ComplexF64, length(modes[1].data), length(modes))
    weighting = modes[1].weighting
    basis = modes[1].basis
    for (col,mode) in enumerate(modes)
        @assert mode.k0 == k0
        @assert mode.weighting == weighting
        @assert mode.basis == basis
        data[:,col] = mode.data
    end
    return HilbertSpace(k0, data, weighting, basis)
end

"""
    shift_k0(space::HilbertSpace, dp::Int, dq::Int)

Shift the basis of the Hilbert space by `dp*b1 + dq*b2`, where `b1` and `b2`
are reciprocal lattice vectors.

This is required when we need the overlaps of modes that are at the same k-point
but in different Brillouin zones.
"""
function shift_k0(space::HilbertSpace, dp::Int, dq::Int)
    data = space.data
    ps = space.basis.ps
    qs = space.basis.qs
    data_new = zeros(ComplexF64, size(data))
    for row in 1:size(data,1)
        # Identify the corresponding (p_new,q_new) indices
        # in the plane wave basis
        p_new = ps[row] - dp
        q_new = qs[row] - dq
        row_new = findfirst(x->x==(p_new,q_new), collect(zip(ps,qs)))
        if !isnothing(row_new)
            data_new[row_new,:] = data[row,:]
        end
    end
    k0_new = space.k0 + dp*space.basis.b1 + dq*space.basis.b2
    return HilbertSpace(k0_new, data_new, space.weighting, space.basis)
end

"""
    normalise(data; weighting=I)

Normalisation of vectors with a weighted inner product.
"""
function normalise(data; weighting=I)
    out = copy(data)
    for n in 1:size(out,2)
        out[:,n] /= sqrt(abs(dot(out[:,n],weighting*out[:,n])))
    end
    return out
end

"""
    orthonormalise(data; weighting=I)

Gram-Schmidt orthonormalisation of vectors with a weighted inner product.
"""
function orthonormalise(data; weighting=I)
    out = copy(data)
    out[:,1] = normalise(out[:,1], weighting=weighting)
    # Gram-Schmidt orthogonalisation
    projector = zeros(size(out,1),size(out,1))
    for n in 2:size(out,2)
        projector += out[:,n-1]*out[:,n-1]' * weighting
        out[:,n] = out[:,n] - projector*out[:,n]
        out[:,n] = normalise(out[:,n], weighting=weighting)
    end
    return out
end


text:
module TreeShakeTest

using FileIO

function foo()
    f, io = mktemp()
    FileIO.save("$(f).png", rand(UInt8, 100, 100))
end

end # module


text:
@testset "794.valid-tic-tac-toe-state.jl" begin
    @test validTicTacToe([
        'O' '_' '_'
        '_' '_' '_'
        '_' '_' '_'
    ]) == false

    @test validTicTacToe([
        'X' 'O' 'X'
        '_' 'X' '_'
        '_' '_' '_'
    ]) == false

    @test validTicTacToe([
        'X' 'X' 'X'
        '_' '_' '_'
        'O' 'O' 'O'
    ]) == false

    @test validTicTacToe([
        'X' 'O' 'X'
        'O' '_' 'O'
        'X' 'O' 'X'
    ]) == true
end


text:
function mock_spectrum(n::Int = Int(1e3); use_units::Bool = false)
    wave = range(1e4, 3e4, length = n)
    sigma = 0.1 .* sin.(wave)
    T = 6700
    flux = @. 1e14 / (wave^5 * (exp(1 / (wave * T)) - 1)) ± sigma
    if use_units
        wave *= u"angstrom"
        flux *= u"erg/s/cm^2/angstrom"
    end
    spectrum(wave, flux, name="Test Spectrum")
end

include("redden.jl")


text:
function write_nl_file(f::IO, m::AmplNLMathProgModel)
    write_nl_header(f, m)

    if m.ncon > 0
        write_nl_c_blocks(f, m)
    end

    if has_objective(m)
        write_nl_o_block(f, m)
    end

    if m.ncon > 0
        write_nl_d_block(f, m)
    end

    write_nl_x_block(f, m)

    if m.ncon > 0
        write_nl_r_block(f, m)
    end

    write_nl_b_block(f, m)

    if m.ncon > 0
        write_nl_k_block(f, m)
        write_nl_j_blocks(f, m)
    end

    if has_objective(m)
        write_nl_g_block(f, m)
    end
end

function has_objective(m::AmplNLMathProgModel)
    return !isempty(m.lin_obj) || (isa(m.obj, Expr))
end

function write_nl_header(f, m::AmplNLMathProgModel)
    # Line 1: Always the same
    println(f, "g3 1 1 0")
    # Line 2: vars, constraints, objectives, ranges, eqns, logical constraints
    n_ranges = sum(m.r_codes .== 0)
    n_eqns = sum(m.r_codes .== 4)
    nobj = has_objective(m) ? 1 : 0
    println(f, " $(m.nvar) $(m.ncon) $nobj $n_ranges $n_eqns 0")
    # Line 3: nonlinear constraints, objectives
    nlc = sum(m.conlinearities .== :Nonlin)
    nlo = m.objlinearity == :Nonlin ? 1 : 0
    println(f, " $nlc $nlo")
    # Line 4: network constraints: nonlinear, linear
    println(f, " 0 0")
    # Line 5: nonlinear vars in constraints, objectives, both
    nonlinear_obj = m.varlinearities_obj .== :Nonlin
    nonlinear_con = m.varlinearities_con .== :Nonlin
    nonlinear = (nonlinear_con + nonlinear_obj) .> 0
    nonlinear_both = (nonlinear_con + nonlinear_obj) .> 1
    nlvc = sum(nonlinear_con .> 0)
    nlvo = sum(nonlinear_obj .> 0)
    nlvb = sum(nonlinear_both)
    println(f, " $nlvc $nlvo $nlvb")
    # Line 6: linear network variables; functions; arith, flags
    println(f, " 0 0 0 1")  # flags set to 1 to get suffixes in .sol file
    # Line 7: discrete variables: binary, integer, nonlinear (b,c,o)
    binary = m.vartypes .== :Bin
    integer = m.vartypes .== :Int
    discrete = binary + integer .> 0
    # Julia 0.6 syntax
    nbv = sum(binary + .!nonlinear .> 1)
    niv = sum(integer + .!nonlinear .> 1)
    nlvbi = sum(nonlinear_both + discrete .> 1)
    nlvci = sum(nonlinear_con - nonlinear_obj + discrete .> 1)
    nlvoi = sum(nonlinear_obj - nonlinear_con + discrete .> 1)
    println(f, " $nbv $niv $nlvbi $nlvci $nlvoi")
    # Line 8: nonzeros in Jacobian, gradients
    nzc = sum(m.j_counts)
    nzo = length(m.lin_obj)
    println(f, " $nzc $nzo")
    # Line 9: max name lengths: constraints, variables
    println(f, " 0 0")
    # Line 10: common exprs: b,c,o,c1,o1
    println(f, " 0 0 0 0 0")
end

# Nonlinear constraint trees
function write_nl_c_blocks(f, m::AmplNLMathProgModel)
    for index in 0:(m.ncon - 1)
        i = m.c_index_map_rev[index]
        println(f, "C$index")
        write_nl_expr(f, m, m.constrs[i])
    end
end

# Nonlinear objective tree
function write_nl_o_block(f, m::AmplNLMathProgModel)
    println(f, string("O0 ", sense_to_nl[m.sense]))
    write_nl_expr(f, m, m.obj)
end

# Initial dual guesses - unused
function write_nl_d_block(f, m::AmplNLMathProgModel)
    println(f, "d$(m.ncon)")
    for index in 0:(m.ncon - 1)
        i = m.c_index_map_rev[index]
        println(f, "$index 0")
    end
end

# Initial primal guesses
function write_nl_x_block(f, m::AmplNLMathProgModel)
    println(f, "x$(m.nvar)")
    for index in 0:(m.nvar - 1)
        i = m.v_index_map_rev[index]
        println(f, "$index $(m.x_0[i])")
    end
end

# Constraint bounds
function write_nl_r_block(f, m::AmplNLMathProgModel)
    println(f, "r")
    for index in 0:(m.ncon - 1)
        i = m.c_index_map_rev[index]
        lower = m.g_l[i]
        upper = m.g_u[i]
        rel = m.r_codes[i]
        if rel == 0
            println(f, "$rel $lower $upper")
        elseif rel == 1
            println(f, "$rel $upper")
        elseif rel == 2
            println(f, "$rel $lower")
        elseif rel == 3
            println(f, "$rel")
        elseif rel == 4
            println(f, "$rel $lower")
        end
    end
end

# Variable bounds
function write_nl_b_block(f, m::AmplNLMathProgModel)
    println(f, "b")
    for index in 0:(m.nvar - 1)
        i = m.v_index_map_rev[index]
        lower = m.x_l[i]
        upper = m.x_u[i]
        if lower == -Inf
            if upper == Inf
                println(f, "3")
            else
                println(f, "1 $upper")
            end
        else
            if lower == upper
                println(f, "4 $lower")
            elseif upper == Inf
                println(f, "2 $lower")
            else
                println(f, "0 $lower $upper")
            end
        end
    end
end

# Jacobian counts
function write_nl_k_block(f, m::AmplNLMathProgModel)
    println(f, "k$(m.nvar - 1)")
    total = 0
    for index = 0:(m.nvar - 2)
        i = m.v_index_map_rev[index]
        total += m.j_counts[i]
        println(f, total)
    end
end

# Linear constraint expressions
function write_nl_j_blocks(f, m::AmplNLMathProgModel)
    for index in 0:(m.ncon - 1)
        i = m.c_index_map_rev[index]
        num_vars = length(m.lin_constrs[i])
        # Only print linear blocks with vars, .nl file is malformed otherwise
        if num_vars > 0
            println(f, "J$index $num_vars")

            # We need to output .nl index and constraint coeff, ordered by .nl
            # index. `lin_constrs` contains our variable index as the key and
            # constraint coeff as the value

            # Assemble tuples of (.nl index, constraint value)
            output = collect(zip(
                (m.v_index_map[j] for j in keys(m.lin_constrs[i])),
                values(m.lin_constrs[i]))
            )

            # Loop through output in .nl index order
            for (index2, value) in sort(output)
                println(f, "$index2 $value")
            end
        end
    end
end

# Linear objective expression
function write_nl_g_block(f, m::AmplNLMathProgModel)
    println(f, string("G0 ", length(m.lin_obj)))
    for index in 0:(m.nvar - 1)
        i = m.v_index_map_rev[index]
        if i in keys(m.lin_obj)
            println(f, "$index $(m.lin_obj[i])")
        end
    end
end

# Convert an expression tree (with .nl formulae only) to .nl format
write_nl_expr(f, m, c) = println(f, string(c))
# Handle numerical constants e.g. pi
write_nl_expr(f, m, c::Symbol) =  write_nl_expr(f, m, float(eval(c)))
function write_nl_expr(f, m, c::Real)
    println(f, nl_number(c == round(Integer, c) ? round(Integer, c) : c))
end
write_nl_expr(f, m, c::LinearityExpr) = write_nl_expr(f, m, c.c)
function write_nl_expr(f, m, c::Expr)
    if c.head == :ref
        # Output variable as `v$index`
        if c.args[1] == :x
            @assert isa(c.args[2], Int)
            println(f, nl_variable(m.v_index_map[c.args[2]]))
        else
            error("Unrecognized reference expression $c")
        end
    elseif c.head == :call
        # Output function as `o$opcode`
        println(f, nl_operator(c.args[1]))
        if c.args[1] in nary_functions
            # Output nargs on subsequent line if n-ary function
            println(f, (string(length(c.args) - 1)))
        end
        map(arg -> write_nl_expr(f, m, arg), c.args[2:end])

    elseif c.head == :comparison
        # .nl only handles binary comparison
        @assert length(c.args) == 3
        # Output comparison type first, followed by args
        println(f, nl_operator(c.args[2]))
        map(arg -> write_nl_expr(f, m, arg), c.args[1:2:end])

    elseif c.head in [:&&, :||]
        # Only support binary and/or for now
        @assert length(c.args) == 2
        println(f, nl_operator(c.head))
        map(arg -> write_nl_expr(f, m, arg), c.args)
    else
        error("Unrecognized expression $c")
    end
end

nl_variable(index::Integer) = "v$index"
nl_number(value::Real) = "n$value"

function nl_operator(operator::Symbol)
    if !haskey(func_to_nl, operator)
        error("translation of the function \"$operator\" to NL is not defined")
    end
    return "o$(func_to_nl[operator])"
end


text:
module Spinners

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export spinner

register_normal_element("q__spinner", context = @__MODULE__)

for spinner in ["audio", "ball", "bars", "box", "clock", "comment", "cube", "dots", "facebook", "gears", "grid",
                "hearts", "hourglass", "infinity", "ios", "orbit", "oval", "pie", "puff", "radio", "rings", "tail"]
  register_normal_element("q__spinner__$spinner", context = @__MODULE__)
end

function spinner(spinner_type::Union{String,Symbol} = "",
                  args...;
                  wrap::Function = StippleUI.DEFAULT_WRAPPER,
                  kwargs...)

  wrap() do
    getfield(@__MODULE__, Symbol("q__spinner$(isempty(string(spinner_type)) ? "" : "__")$spinner_type"))(args...; kwargs...)
  end
end

end


text:
import PathWeightSampling
using Test

system = PathWeightSampling.gene_expression_system(dtimes=0:0.2:1.0)
conf = PathWeightSampling.generate_configuration(system)
ensemble = PathWeightSampling.MarginalEnsemble(system)

# Test the shooting moves
jump_problem = ensemble.jump_problem
original_traj = conf.s_traj
new_traj = copy(conf.s_traj)
@test original_traj == new_traj
PathWeightSampling.shoot_forward!(new_traj, original_traj, jump_problem, 0.5)
@test original_traj != new_traj
@test issorted(new_traj.t) && allunique(new_traj.t)
@test length(new_traj.t) == length(new_traj.u) == length(new_traj.i)

new_traj = copy(conf.s_traj)
@test original_traj == new_traj
PathWeightSampling.shoot_backward!(new_traj, original_traj, jump_problem, 0.5)
@test original_traj != new_traj
@test issorted(new_traj.t) && allunique(new_traj.t)
@test length(new_traj.t) == length(new_traj.u) == length(new_traj.i)

mchain = PathWeightSampling.chain(ensemble; θ=1.0)
@test PathWeightSampling.energy(conf, mchain) < 0

# Test Thermodynamic Integration

alg = PathWeightSampling.TIEstimate(100, 16, 1000)
result = PathWeightSampling.simulate(alg, conf, ensemble)
@test all(PathWeightSampling.log_marginal(result) .>= 0)


text:
## Lebesgue Measure - Volume of box
## ================================
measure(b::Boxes) = prod([measure(b[dim]) for dim in dims(b)])
logmeasure(b::Boxes) = sum([log(measure(b[dim])) for dim in dims(b)])
measures{B<:Boxes}(bs::Vector{B}) = [measure(b) for b in bs]


text:
function main_recursion(
    Ns::Int64, N::Int64, expβϵ::Array{T,1}, FLAG_OccNum::Bool
) where {T<:FloatType}

    # recursion setup
    expβμ = (abs(expβϵ[Ns - N + 1]) +
            abs(expβϵ[Ns - N])) / 2       # chemical potential
    expβϵμ = expβϵ / expβμ                # shifted value

    # recursions for PFs
    Z̃ = APF_recursion(Ns, expβϵμ)
    # use the Z_1 instead of Z_0 as Z̃_1 is more stable numerically
    Z₁ = real(sum(expβϵ))
    Z = Z̃[N + 1] / Z̃[2] * (expβμ) ^ (N - 1) * Z₁

    # skip the occupation number recursion
    if FLAG_OccNum == false
        return Z
    end

    # num. of energy levels below the Fermi level
    # use this formula to ensure complex conjugate pairs are in the same section
    N_below = sum(abs.(expβϵμ) .> 1)

    # separate recursions for occupancies of energy levels above/below the Fermi level
    n̄ = complex(zeros(Ns))
    n₁ = expβϵ / Z₁
    n̄_above = APF_occupancy_recursion(N, expβϵμ[1 : Ns - N_below], Z̃, n₁)
    n̄_below = APF_occupancy_recursion_reverse(Ns, N, expβϵμ[Ns - N_below + 1 : Ns], Z̃)
    # then concatenate
    n̄ = vcat(n̄_above, n̄_below)

    Z, n̄

end

function APF_recursion(
    Ns::Int64, expβϵμ::Array{T,1}
) where {T<:FloatType}

    ν1 = expβϵμ ./ (1 .+ expβϵμ)      # occupation probability in GCE
    ν2 = 1 ./ (1 .+ expβϵμ)           # hole probability in GCE
    Z̃ = complex(zeros(Ns + 1, Ns))    # APF

    # Initialization
    Z̃[1, 1] = ν2[1]
    Z̃[2, 1] = ν1[1]
    for i = 2 : Ns                    # iteration over energy levels
        Z̃[1, i] = ν2[i] * Z̃[1, i - 1]
        for j = 2 : i + 1             # iteration over particles
            Z̃[j, i] = ν2[i] * Z̃[j, i - 1] + ν1[i] * Z̃[j - 1, i - 1]
        end
    end

    real(Z̃[:, Ns])                    # PF can be proven to be strictly real

end

# auxiliary partition function recursion
# considering complex conjugate pairs explicitly
function APF_recursionCompConj(
    Ns::Int64, expβϵμ::Array{T,1}
) where {T<:FloatType}

    # find all complex conj. pairs (ccp) first
    expβϵμ_ccp = [x for x in expβϵμ if imag(x) > 1e-15]
    expβϵμ_real = real([x for x in expβϵμ if abs(imag(x)) <= 1e-15])
    num_ccp = length(expβϵμ_ccp)

    Z̃ = zeros(Ns + 1, Ns - num_ccp)         # APF

    # do the real spectra first
    ν = expβϵμ_real ./ (1 .+ expβϵμ_real)
    Z̃[1, 1] = 1 - ν[1]
    Z̃[2, 1] = ν[1]
    for i = 2 : Ns - 2 * num_ccp            # iteration over real energy levels
        Z̃[1, i] = (1 - ν[i]) * Z̃[1, i - 1]
        for j = 2 : i + 1                   # iteration over particles
            Z̃[j, i] = (1 - ν[i]) * Z̃[j, i - 1] + ν[i] * Z̃[j - 1, i - 1]
        end
    end

    if num_ccp == 0
        return Z̃[:, Ns]
    end

    # then the complex conjugate pairs part
    ν_denom = 1 .+ 2 * real(expβϵμ_ccp) .+ real(expβϵμ_ccp).^2 .+ imag(expβϵμ_ccp).^2
    ν1 = 2 * real(expβϵμ_ccp) ./ ν_denom
    ν2 = (real(expβϵμ_ccp).^2 .+ imag(expβϵμ_ccp).^2) ./ ν_denom
    Nsc = Ns - 2 * num_ccp                  # num. of spectra that have been counted
    for i = 1 : num_ccp                     # iteration over ccp energy levels
        Z̃[1, Nsc + i] = (1 - ν1[i] - ν2[i]) * Z̃[1, Nsc + i - 1]
        Z̃[2, Nsc + i] = (1 - ν1[i] - ν2[i]) * Z̃[2, Nsc + i - 1] +
        ν1[i] * Z̃[1, Nsc + i - 1]
        for j = 3 : Nsc + 2 * i + 1         # iteration over particles
            Z̃[j, Nsc + i] = (1 - ν1[i] - ν2[i]) * Z̃[j, Nsc + i - 1] +
            ν1[i] * Z̃[j - 1, Nsc + i - 1] + ν2[i] * Z̃[j - 2, Nsc + i - 1]
        end
    end

    Z̃[:, Ns - num_ccp]

end

# Occupation number recursion (above Fermi level) using rescaled spectrum and APF
function APF_occupancy_recursion(
    N::Int64, expβϵμ::Array{T1,1}, Z̃::Array{T2,1}, n₁::Array{T3,1}
) where {T1<:FloatType, T2<:FloatType, T3<:FloatType}

    # Initialization
    Ñs = length(expβϵμ)
    n̄ = complex(zeros(N + 1, Ñs))
    # use n₁ as a correction
    n̄[2, :] = n₁[1 : Ñs]

    for i = 2 : N
        n̄[i + 1, :] = real(Z̃[i] / Z̃[i + 1]) * expβϵμ .* (1 .- n̄[i, :])
        # Truncate the values that are smaller than 10^-10
        n̄[i + 1, :] = n̄[i + 1, :] .* (abs.(real(n̄[i + 1, :])) .> 1e-10)
    end

    n̄[N + 1, :]

end

# Occupation number reverse recursion (below Fermi level) using rescaled spectrum and APF
function APF_occupancy_recursion_reverse(
    Ns::Int64, N::Int64, expβϵμ::Array{T1,1}, Z̃::Array{T2,1}
) where {T1<:FloatType, T2<:FloatType}

    # Initialization
    Ñs = length(expβϵμ)
    N_rev = Ns - N                  # num. of reverse recursions
    n̄ = complex(ones(N_rev, Ñs))

    for i = 1 : N_rev - 1
        n̄[i + 1, :] = real(Z̃[Ns - i + 1] / Z̃[Ns - i]) * n̄[i, :] ./ expβϵμ
        # Truncate the values that are smaller than 10^-10
        n̄[i + 1, :] = n̄[i + 1, :] .* (abs.(real(n̄[i + 1, :])) .> 1e-10)
        n̄[i + 1, :] = 1 .- n̄[i + 1, :]
    end

    n̄[N_rev, :]

end

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
#=##############################################################################
# DESCRIPTION
    Types defining maneuvers of flight vehicles.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

################################################################################
# ABSTRACT MANEUVER TYPE
################################################################################
"""
    `AbstractManeuver{N, M}`

`N` indicates the number of tilting systems in this maneuver, while and `M`
indicates the number of rotor systems.

Every implementation of `AbstractManeuver` must have the properties:

 * `angle::NTuple{N, Function}` where `angle[i](t)` returns the angle of the
        i-th tilting system at time `t` (t is nondimensionalized by the total
        time of the maneuver, from 0 to 1, beginning to end).
 * `RPM::NTuple{M, Function}` where `RPM[i](t)` returns the normalized RPM of
        the i-th rotor system at time `t`. This RPM values are normalized by the
        an arbitrary RPM value (usually RPM in hover or cruise).
"""
abstract type AbstractManeuver{N, M} end

##### FUNCTIONS REQUIRED IN IMPLEMENTATIONS ####################################
"""
    `calc_dV(maneuver::AbstractManeuver, vehicle::Vehicle, t, dt, ttot, Vref)`

Returns the change in velocity `dV=[dVx, dVy, dVz]` (m/s) of `vehicle` performing
`maneuver` at time `t` (s) after a time step `dt` (s).  `Vref` and `tot` are the
reference velocity and the total time at which this maneuver is being performed,
respectively. `dV` is in the global reference system.
"""
function calc_dV(self::AbstractManeuver, vehicle::AbstractVehicle, t::Real,
                                            dt::Real, ttot::Real, Vref::Real)
    error("$(typeof(self)) has no implementation yet!")
end

"""
    `calc_dw(maneuver::AbstractManeuver, vehicle::Vehicle, t, dt, Vref, ttot)`

Returns the change in angular velocity `dW=[dWx, dWy, dWz]` (about global
axes, in radians) of `vehicle` performing `maneuver` at time `t` (s) after a
time step `dt` (s). `ttot` is the total time at which this maneuver is to be
performed.
"""
function calc_dW(self::AbstractManeuver, vehicle::AbstractVehicle, t::Real,
                                                        dt::Real, ttot::Real)
    error("$(typeof(self)) has no implementation yet!")
end

##### COMMON FUNCTIONS  ########################################################
"""
    `get_ntltsys(self::AbstractManeuver)`
Return number of tilting systems.
"""
get_ntltsys(self::AbstractManeuver) = typeof(self).parameters[1]

"""
    `get_nrtrsys(self::AbstractManeuver)`
Return number of rotor systems.
"""
get_nrtrsys(self::AbstractManeuver) = typeof(self).parameters[2]

"""
    `get_angle(maneuver::AbstractManeuver, i::Int, t::Real)`

Returns the angle (in degrees) of the i-th tilting system at the non-dimensional
time t.
"""
function get_angle(self::AbstractManeuver, i::Int, t::Real)
    if i<=0 || i>get_ntltsys(self)
        error("Invalid tilting system #$i (max is $(get_ntltsys(self))).")
    end
    if t<0 || t>1
        @warn("Got non-dimensionalized time $(t).")
    end
    return self.angle[i](t)
end

"""
    `get_angles(maneuver::AbstractManeuver, t::Real)`

Returns the angle (in degrees) of every tilting systems at the non-dimensional
time t.
"""
get_angles(self::AbstractManeuver, t::Real) = Tuple(a(t) for a in self.angle)

"""
    `get_RPM(maneuver::AbstractManeuver, i::Int, t::Real)`

Returns the normalized RPM of the i-th rotor system at the non-dimensional time
t.
"""
function get_RPM(self::AbstractManeuver, i::Int, t::Real)
    if i<=0 || i>get_nrtrsys(self)
        error("Invalid rotor system #$i (max is $(get_nrtrsys(self))).")
    end
    if t<0 || t>1
        @warn("Got non-dimensionalized time $(t).")
    end
    return self.RPM[i](t)
end

"""
    `get_RPMs(maneuver::AbstractManeuver, t::Real)`

Returns the normalized RPM of every rotor systems at the non-dimensional time
t.
"""
get_RPMs(self::AbstractManeuver, t::Real) = Tuple(rpm(t) for rpm in self.RPM)

##### COMMON INTERNAL FUNCTIONS  ###############################################

##### END OF ABSTRACT MANEUVER #################################################

################################################################################
# KINEMATIC MANEUVER TYPE
################################################################################
"""
    `KinematicManeuver(angle, RPM, Vvehicle, anglevehicle)`

A vehicle maneuver where the kinematic are prescribed.

# ARGUMENTS
* `angle::NTuple{N, Function}` where `angle[i](t)` returns the angles
        `[Ax, Ay, Az]` (in degrees)of the i-th tilting system at time `t` (t is
        nondimensionalized by the total time of the maneuver, from 0 to 1,
        beginning to end).
* `RPM::NTuple{M, Function}` where `RPM[i](t)` returns the normalized RPM of
        the i-th rotor system at time `t`. This RPM values are normalized by the
        an arbitrary RPM value (usually RPM in hover or cruise).
* `Vvehicle::Function` where `Vvehicle(t)` returns the normalized vehicle
        velocity `[Vx, Vy, Vz]` at the normalized time `t`. Velocity is
        normalized by a reference velocity (typically, cruise velocity).
* `anglevehicle::Function` where `anglevehicle(t)` returns the angles
        `[Ax, Ay, Az]` (in degrees) of the vehicle relative to the global
        coordinate system at the normalized time `t`.
"""
struct KinematicManeuver{N, M} <: AbstractManeuver{N, M}
    angle::NTuple{N, Function}
    RPM::NTuple{M, Function}
    Vvehicle::Function
    anglevehicle::Function
end

# # Implicit N and M constructor
# KinematicManeuver(a::NTuple{N, Function}, b::NTuple{M, Function},
#                     c::Function, d::Function
#                  ) where {N, M} = KinematicManeuver{N, M}(a, b, c, d)

##### FUNCTIONS  ###############################################################
function calc_dV(self::KinematicManeuver, vehicle::AbstractVehicle, t::Real,
                                            dt::Real, ttot::Real, Vref::Real)
    return Vref * (self.Vvehicle((t+dt)/ttot) - self.Vvehicle(t/ttot))
end

function calc_dW(self::KinematicManeuver, vehicle::AbstractVehicle, t::Real,
                                                         dt::Real, ttot::Real)
    prev_W = (self.anglevehicle(t/ttot) - self.anglevehicle((t-dt)/ttot)) / dt
    cur_W = (self.anglevehicle((t+dt)/ttot) - self.anglevehicle(t/ttot)) / dt
    return pi/180 * (cur_W - prev_W)
end

##### INTERNAL FUNCTIONS  ######################################################

##### END OF KINEMATICMANEUVER  ################################################

"""
    MassProp(m, Ixx, Iyy, Izz, Ixz, Ixy, Iyz)

Mass and moments of inertia in the body frame.
Ixx = int(y^2 + z^2, dm)
Ixz = int(xz, dm)

We can model our quadcopter as two thin uniform rods crossed at the origin with a point mass
(motor) at the end of each. With this in mind, it’s clear that the symmetries result in a diagonal
inertia matrix
"""
struct MassProp{TF}
    m::TF
    Ixx::TF
    Iyy::TF
    Izz::TF
end

"""
    ConstantAtmosphere(Wi, Wb, rho, asound, g)

Constant atmospheric properties.
"""
struct Atmosphere{TF,TV}
    Wi::TV
    Wb::TV
    rho::TF
    asound::TF
    g::TF
end

################################################################################
# KINEMATIC MANEUVER TYPE
################################################################################
mutable struct DynamicManeuver{N, M} <: AbstractManeuver{N, M}
    ttot::Real
    tinit::Real
    angle::NTuple{N, Function}
    RPM::NTuple{M, Function}
    RPMref::Real
    mp::MassProp{Float64}
    atm::Atmosphere{Float64,Array{Float64,1}}
    MovementMatrix::Matrix{Float64}
    k::Real
    b::Real
    kd::Real
    l::Real
    dt::Real
    rotors::Int
end

# # Implicit N and M constructor
# DynamicManeuver(a::NTuple{N, Function}, b::NTuple{M, Function}
#                                     ) where {N, M} = DynamicManeuver{N, M}(a, b)

##### FUNCTIONS  ###############################################################
function calc_dV(self::DynamicManeuver, t::Real)
a = Array{Float64}(undef, 3)
thetaprev = Array{Float64}(undef, 3)
xdotprev = Array{Float64}(undef, 3)
omegaprev = Array{Float64}(undef, 3)
for i in 1:4
rpm[i] = self.RPMref*self.RPM[i](t)
end
for i in 1:3
xprev[i] = self.MovementMatix[3+i,selfi.nt]
xdotprev[i] = self.MovementMatrix[6+i,selfi.nt]
thetaprev[i] = self.MovementMatrix[9+i,selfi.nt]
omegaprev[i] = self.MovementMatix[12+i,selfi.nt]
end
a = acceleration(rpm, thetaprev, xdotprev, self.atm, self.k, self.mp)
for i in 1:3
self.MovementMatix[i,selfi.nt+1] = a
end
omegadot = Array{Float64}(undef, 3)
I = [self.mp.Ixx   0.0   0.0;
     0.0     self.mp.Iyy 0.0;
     0.0      0.0  self.mp.Izz]
tau = torque(rpm, self.l, self.k, self.b)
omegadot = inv(I) * (tau - cross(omegaprev, I * omegaprev))
omega = Array{Float64}(undef, 3)
omega = omegaprev + dt * omegadot
for i in 1:3
self.MovementMatix[12+i,selfi.nt+1] = omega
end
thetadot = Array{Float64}(undef, 3)
thetadot = omega2thetadot(thetaprev, omega)
theta = Array{Float64}(undef, 3)
theta = thetaprev + dt * thetadot
for i in 1:3
self.MovementMatix[9+i,selfi.nt+1] = theta
end
xdot = Array{Float64}(undef, 3)
xdot = xdotprev + dt * a
x = Array{Float64}(undef, 3)
x = xprev + dt * xdot
for i in 1:3
self.MovementMatix[3+i,selfi.nt+1] = x
end
dv = xdot - xdotprev
return dv
end

function calc_dw(self::DynamicManeuver, t, nt)
a = Array{Float64}(undef, 3)
thetaprev = Array{Float64}(undef, 3)
xdotprev = Array{Float64}(undef, 3)
omegaprev = Array{Float64}(undef, 3)
for i in 1:4
rpm[i] = self.RPMref*self.RPM[i](t)
end
for i in 1:3
xprev[i] = self.MovementMatix[3+i,nt]
xdotprev[i] = self.MovementMatrix[6+i,nt]
thetaprev[i] = self.MovementMatrix[9+i,nt]
omegaprev[i] = self.MovementMatix[12+i,nt]
end
a = acceleration(rpm, thetaprev, xdotprev, self.atm, self.k, self.mp)
for i in 1:3
self.MovementMatix[i,nt+1] = a
end
omegadot = Array{Float64}(undef, 3)
I = [self.mp.Ixx   0.0   0.0;
     0.0     self.mp.Iyy 0.0;
     0.0      0.0  self.mp.Izz]
tau = torque(rpm, self.l, self.k, self.b)
omegadot = inv(I) * (tau - cross(omegaprev, I * omegaprev))
omega = Array{Float64}(undef, 3)
omega = omegaprev + dt * omegadot
for i in 1:3
self.MovementMatix[12+i,nt+1] = omega
end
thetadot = Array{Float64}(undef, 3)
thetadot = omega2thetadot(thetaprev, omega)
theta = Array{Float64}(undef, 3)
theta = thetaprev + dt * thetadot
for i in 1:3
self.MovementMatix[9+i,nt+1] = theta
end
xdot = Array{Float64}(undef, 3)
xdot = xdotprev + dt * a
x = Array{Float64}(undef, 3)
x = xprev + dt * xdot
for i in 1:3
self.MovementMatix[3+i,nt+1] = x
end
return  omega - omegaprev
end

##### INTERNAL FUNCTIONS  ######################################################
function bodytoinertial(theta)

    Ri = Array{eltype(theta)}(undef,3, 3)

    cphi, ctht, cpsi = cos.([theta[1], theta[2], theta[3]])
    sphi, stht, spsi = sin.([theta[1], theta[2], theta[3]])

    Ri[1, 1] = cphi*cpsi - ctht*sphi*spsi
    Ri[1, 2] = -cpsi*sphi - cphi*ctht*spsi
    Ri[1, 3] = stht*spsi

    Ri[2, 1] = ctht*sphi*cpsi + cphi*spsi
    Ri[2, 2] = cphi*ctht*cpsi - sphi*spsi
    Ri[2, 3] = -cpsi*stht

    Ri[3, 1] = sphi*stht
    Ri[3, 2] = cphi*stht
    Ri[3, 3] = ctht

    return Ri

end

"""
 Total thrust, giving the square rotation velocity and the coefficient k
"""
function thrust(rpm, k)
sum=0
for i in 1:self.rotors
    sum  = sum + rpm[i]
end
    T = [0.0; 0.0; k*sum]
    return T
end

"""
Torque, giving the square rotation velocity, the distance of rotors from the center
and the coefficients k and b
"""
function torque(rpm, l, k, b)
    tau = [l*k*(rpm[1]-rpm[3]); l*k*(rpm[2]-rpm[4]); b*(rpm[1]-rpm[2]+rpm[3]-rpm[4])]
    return tau
end

"""
Computation of the linear acceleration
"""
function acceleration(rpm, theta, xdot, atm, k, mp)
    gravity= [0.0; 0.0; -atm.g]
    R = bodytoinertial(theta)
    Tb = thrust(rpm, k)
    Ti = R * Tb
    Fd = -kd * xdot
    a = gravity + 1/mp.m  * Ti + 1/mp.m  * Fd
return a
end

"""
    UniformGravitationalField()

Assumes center of mass and center of gravity are coincident.
"""
#struct UniformGravitationalField <: AbstractInertialModel end

function omega2thetadot(theta,omega)
    ct, cp = cos.([theta[2], theta[1]])
    st, sp = sin.([theta[2], theta[1]])

    C = [1.0  0.0  -st;
        0.0   cp  ct*sp
        0.0  -sp  ct*cp]

    thetadot = inv(C)*omega
    return thetadot
end
##### END OF KINEMATICMANEUVER  ################################################

#


text:
using GeometryTypes, ColorTypes
using Test
using Test: @inferred

# 0.7 Base still contains a (deprecated) contains; resolve ambiguity
if VERSION < v"1.0-"
    using GeometryTypes: contains
end

@testset "GeometryTypes" begin
    include("baseutils.jl")
    include("convexhulls.jl")
    include("polygons.jl")
    include("hyperrectangles.jl")
    include("faces.jl")
    include("meshes.jl")
    include("distancefields.jl")
    include("primitives.jl")
    include("decompose.jl")
    include("simplerectangle.jl")
    include("hypersphere.jl")
    include("typeutils.jl")
    include("simplices.jl")
    include("gjk.jl")
    include("lines.jl")
    include("polygons.jl")
    include("cylinder.jl")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c83a4764-1d65-11eb-0211-3f30ca230c32
using Pkg, DrWatson

# ╔═╡ cf3e1fd6-1d65-11eb-043f-b1e0c0d3620e
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Roots
end

# ╔═╡ 8be3ec46-1d49-11eb-0cb0-bb04f5925358
md"## Listing 1.6"

# ╔═╡ 837f96be-1d3e-11eb-1df7-878544276c0a
function bubbleSort!(a)
    n = length(a)
    for i in 1:n-1
        for j in 1:n-i
            if a[j] > a[j+1]
                a[j], a[j+1] = a[j+1], a[j]
            end
        end
    end
    return a
end

# ╔═╡ b14809a4-1d3f-11eb-1567-0364629fc97d
data = [65, 51, 32, 12, 23, 84, 68, 1]

# ╔═╡ b1485852-1d3f-11eb-1ebe-09648d32a6f8
bubbleSort!(data)

# ╔═╡ 3e65dfd4-1d66-11eb-275a-8747968227c7
md"## Listing 1.7"

# ╔═╡ 13f2150a-1d4e-11eb-30bc-c1bf235f675b
function polynomialGenerator(a...)
    n = length(a)-1
    poly =  function(x)
                return sum([a[i+1]*x^i for i in 0:n])
            end
    return poly
end

# ╔═╡ 2a3e53a6-1d66-11eb-0b8f-03b03d500b4f
polynomial = polynomialGenerator(1, 3, -10)

# ╔═╡ 2a3ea644-1d66-11eb-22d0-254c46fbb49d
zeroVals = find_zeros(polynomial, -10 ,10)

# ╔═╡ 2a3f3f3c-1d66-11eb-0ee9-ff23b3c7804f
"Zeros of the function f(x): ", zeroVals

# ╔═╡ Cell order:
# ╟─8be3ec46-1d49-11eb-0cb0-bb04f5925358
# ╠═c83a4764-1d65-11eb-0211-3f30ca230c32
# ╠═cf3e1fd6-1d65-11eb-043f-b1e0c0d3620e
# ╠═837f96be-1d3e-11eb-1df7-878544276c0a
# ╠═b14809a4-1d3f-11eb-1567-0364629fc97d
# ╠═b1485852-1d3f-11eb-1ebe-09648d32a6f8
# ╟─3e65dfd4-1d66-11eb-275a-8747968227c7
# ╠═13f2150a-1d4e-11eb-30bc-c1bf235f675b
# ╠═2a3e53a6-1d66-11eb-0b8f-03b03d500b4f
# ╠═2a3ea644-1d66-11eb-22d0-254c46fbb49d
# ╠═2a3f3f3c-1d66-11eb-0ee9-ff23b3c7804f


text:
# ------------------------------------------------------------------------------------------
# ## Intro to neurons
#
# At this point, we know how to build models and have a computer automatically learn how to
# match the model to data. This is the core of how any machine learning method works.
#
# Now, let's narrow our focus and look at **neural networks**. Neural networks (or "neural
# nets", for short) are a specific choice of a **model**. It's a network made up of
# **neurons**; this, in turn, leads to the question, "what is a neuron?"
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Models with multiple inputs
#
# So far, we have been using the sigmoid function as our model. One of the forms of the
# sigmoid function we've used is
#
# $$\sigma_{w, b}(x) = \frac{1}{1 + \exp(-wx + b)},$$
#
# where `x` and `w` are both single numbers. We have been using this function to model how
# the amount of the color green in an image (`x`) can be used to determine if an image shows
# an apple or a banana.
#
# But what if we had multiple data features we wanted to fit?
#
# We could then extend our model to include multiple features like
#
# $$\sigma_{\mathbf{w},b}(\mathbf{x}) = \frac{1}{1 + \exp(-w_1 x_1 - w_2 x_2 - \cdots - w_n
# x_n + b)}$$
#
# Note that now $\mathbf{x}$ and $\mathbf{w}$ are vectors with many components, rather than
# a single number.
#
# For example, $x_1$ might be the amount of the color green in an image, $x_2$ could be the
# height of the object in the picture, $x_3$ could be the width, and so forth. We can add
# information for as many features as we have! Our model now has more parameters, but the
# same idea of gradient descent ("rolling the ball downhill") will still work to train our
# model.
#
# This version of the sigmoid model that takes multiple inputs is an example of a
# **neuron**.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In the video, we see that one huge class of learning techniques is based around neurons,
# that is, *artificial neurons*. These are caricatures of real, biological neurons. Both
# *artificial* and *biological* neurons have several inputs $x_1, \ldots, x_n$, and a single
# output, $y$. Schematically they look like this:
# ------------------------------------------------------------------------------------------

include("draw_neural_net.jl")

number_inputs, number_neurons = 4, 1

draw_network([number_inputs, number_neurons])

# ------------------------------------------------------------------------------------------
# We should read this as showing how information flows from left to right:
# - 4 pieces of input information arrive (shown in green on the left);
#
# - the neuron (shown in red on the right) receives all the inputs, processes them, and
# outputs a single number to the right.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In other words, a neuron is just a type of function that takes multiple inputs and returns
# a single output.
#
# The simplest interesting case that we will look at in this notebook is when there are just
# two pieces of input data:
# ------------------------------------------------------------------------------------------

draw_network([2, 1])

# ------------------------------------------------------------------------------------------
# Each link between circles above represents a **weight** $w$ that can be modified to allow
# the neuron to learn, so in this case the neuron has two weights, $w_1$ and $w_2$.
#
# The neuron also has a single bias $b$, and an **activation function**, which we will take
# to be the $\sigma$ sigmoidal function that we have been using. (Note that other activation
# functions can be used!)
#
# Let's call our neuron $f_{w_1,w_2, b}(x_1, x_2)$, where
#
# $$f_{w_1,w_2, b}(x_1, x_2) := \sigma(w_1 x_1 + w_2 x_2 + b).$$
#
# Note that $f_{w_1,w_2, b}(x_1, x_2)$ has 3 parameters: two weights and a bias.
#
# To simplify the notation, and to prepare for more complicated scenarios later, we put the
# two weights into a vector, and the two data values into another vector:
#
# $$
# \mathbf{w} = \begin{pmatrix} w_1 \\ w_2 \end{pmatrix};
# \qquad
# \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}.
# $$
#
# We thus have
#
# $$f_{\mathbf{w}, b}(\mathbf{x}) = \sigma(\mathbf{w} \cdot \mathbf{x} + b),$$
#
# where the dot product $\mathbf{w} \cdot \mathbf{x}$ is an abbreviated syntax for $w_1 x_1
# + w_2 x_2$.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Declare the function `f(x, w, b)` in Julia. `f` should take vectors `x` and `w` as vectors
# and `b` as a scalar. Furthermore `f` should call
#
# ```julia
# σ(x) = 1 / (1 + exp(-x))
# ```
#
# What output do you get for
#
# ```julia
# f(3, 4, 5)
# ```
# ?
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/11. ML - Intro to neurons.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using LinearAlgebra
using Test

import MathOptInterface
const MOI = MathOptInterface

import Linda
import GLPK

const testdir = dirname(@__FILE__)

@testset "Linda" begin
    include(joinpath(testdir, "env.jl"))
    include(joinpath(testdir, "master.jl"))
    include(joinpath(testdir, "colgen.jl"))
end


text:
# ================
# Homogeneous case
# ================

# homogeneous case without invariant, non-recursive implementation
function reach_homog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
                          Ω0::Hyperrectangle{N, VNC, VNR},
                          Φ::AbstractMatrix,
                          NSTEPS::Integer,
                          δ::N,
                          X::Universe,
                          recursive::Val{false},
                          Δt0::TN) where {N, TN, VNC, VNR}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    # preallocations
    n = dim(Ω0)
    c = [VNC(undef, n) for _ in 1:NSTEPS]
    r = [VNR(undef, n) for _ in 1:NSTEPS]

    copy!(c[1], Ω0.center)
    copy!(r[1], Ω0.radius)

    # cache for powers of Φ
    Φ_power_k = similar(Φ)
    copyto!(Φ_power_k, Φ)
    Φ_power_k_cache = similar(Φ_power_k)
    Φ_power_k_abs = similar(Φ_power_k)

    k = 2
    @inbounds while k <= NSTEPS
        # compute the overapproximation of Φ^k * Ω0 with a hyperrectangle
        mul!(c[k], Φ_power_k, c[1])
        Φ_power_k_abs .= abs.(Φ_power_k)
        mul!(r[k], Φ_power_k_abs, r[1])
        Hₖ = Hyperrectangle(c[k], r[k], check_bounds=false)

        Δt += δ
        F[k] = ReachSet(Hₖ, Δt)

        mul!(Φ_power_k_cache, Φ_power_k, Φ)
        copyto!(Φ_power_k, Φ_power_k_cache)
        k += 1
    end
    return F
end

# homogeneous case without invariant, recursive implementation
# note that this implementation suffers from wrapping
function reach_homog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
                          Ω0::Hyperrectangle{N, VNC, VNR},
                          Φ::AbstractMatrix,
                          NSTEPS::Integer,
                          δ::Float64,
                          X::Universe,
                          recursive::Val{true},
                          Δt0::TN) where {N, TN, VNC, VNR}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    k = 2
    @inbounds while k <= NSTEPS
        Rₖ = overapproximate(Φ * set(F[k-1]), Hyperrectangle)
        Δt += δ
        F[k] = ReachSet(Rₖ, Δt)
        k += 1
    end
    return F
end

# homogeneous case with an invariant, non-recursive implementation
function reach_homog_BOX!(F::Vector{ReachSet{N, Hyperrectangle{N, VNC, VNR}}},
                          Ω0::Hyperrectangle{N, VNC, VNR},
                          Φ::AbstractMatrix,
                          NSTEPS::Integer,
                          δ::N,
                          X::LazySet,
                          recursive::Val{false},
                          Δt0::TN) where {N, TN, VNC, VNR}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    # preallocations
    n = dim(Ω0)
    r = [VNR(undef, n) for _ in 1:NSTEPS]
    c = [VNC(undef, n) for _ in 1:NSTEPS]

    copy!(c[1], Ω0.center)
    copy!(r[1], Ω0.radius)

    # cache for powers of Φ
    Φ_power_k = similar(Φ)
    copyto!(Φ_power_k, Φ)
    Φ_power_k_cache = similar(Φ_power_k)
    Φ_power_k_abs = similar(Φ_power_k)

    k = 2
    @inbounds while k <= NSTEPS
      # compute the overapproximation of Φ^k * Ω0 with a hyperrectangle
      mul!(c[k], Φ_power_k, c[1])
      Φ_power_k_abs .= abs.(Φ_power_k)
      mul!(r[k], Φ_power_k_abs, r[1])
      Hₖ = Hyperrectangle(c[k], r[k], check_bounds=false)

      _is_intersection_empty(X, Hₖ) && break
      Δt += δ
      F[k] = ReachSet(Hₖ, Δt)

      mul!(Φ_power_k_cache, Φ_power_k, Φ)
      copyto!(Φ_power_k, Φ_power_k_cache)
      k += 1
    end
    if k < NSTEPS + 1
        resize!(F, k-1)
    end
    return F
end


text:
# Author: AJ Tyler
# Date: 24/01/22
module Comparison # This module contains functions which compare configurations using their CNA profiles.
include("io.jl") # This file contains the IO module.
using .IO # This module contains functions that read in and write out files.
using Statistics # This module allows easy mean and standard deviation.
export compareConfigs,groupConfigs,configStats # These functions ared used to compare configurations.

"""
Finds the most similar configuration to every configuration over all rCut values based on their CNA profile.

Inputs:
rCutRange: (Vector{Float64}) A range of cut-off radii to test, in units of equilibriumBondLength.
energies: (Vector{Float64}) The energy of each configuration in each the file.
totalProfiles: Array{Dict{String,Int}} 2D array containg the totalCNA profiles of each configuration, for each rCut value.
atomicProfiles: Array{Dict{String,Int}} 3D array containing the CNA profile of each configuration, for each rCut value, for each atom.
configDir: (String) The path to (or name of if in same location as this program) the directory containing the .xyz files.
filename: (String) The prefix of the file.
similarityMeasure: (String) The method of comparing configurations, either "atomic" or "total" CNA.
M: (Int64) The number of rCut values to test.
L: (Int64) The number of configurations in the file.
N: (Int64) The number of atoms in each configuration.

Outputs:
similarConfigs: Vector{Dict{Int64,Vector{Float64}}} A vector of dictionaries which map the most similar configurations to the rCut values
													for which it is maximised, for each configuraiton.
maxSims: Vector{Float64} A Vector of the maximum similarity score for each configuration when compared to every other configuration
						 over all rCut values.
"""
function compareConfigs(rCutRange,energies,totalProfiles,atomicProfiles,configDir,filename,similarityMeasure,M,L,N)
	# Initialise vector to contain most similar configurations for each configuration
	similarConfigs = Vector{Dict{Int64,Vector{Float64}}}(undef,L)
	maxSims = zeros(Float64,1,L) # Initialise the maximum similarity score for each configuration
	for i in 1:M # For all rCut values
		writeProfiles(rCutRange,configDir,filename,energies,totalProfiles,i,L) # Write out the total CNA profiles
		if L == 1
			continue
		end
		atomicProfilesDict = convertAtomicProfiles(atomicProfiles,i,L) # Convert atomic CNA profile format
		# For each pair of configurations
		for configA in 1:L-1
			for configB in configA+1:L
				# Compute the similarity for the cluster pair
				similarity = similarityScore(configA, configB, totalProfiles,atomicProfilesDict,similarityMeasure,N,i)
				# Update the most similar configurations to configA
				if (similarity > maxSims[configA]) # If new similarity score is greater than previous highest
					similarConfigs[configA] = Dict{Int64,Vector{Float64}}(configB=>Vector{Float64}([rCutRange[i]])) # Create new vector, containing configB
					maxSims[configA] = similarity # Update maximum similarity score
				elseif (similarity == maxSims[configA]) # If the new similarity is the same as the previous highest
					if (haskey(similarConfigs[configA],configB)) # If the same configuration pair is already most similar
						# Add the current rCut value to the dictionary
						push!(similarConfigs[configA], configB => push!(get!(similarConfigs[configA],configB,0),rCutRange[i]))
					else
						push!(similarConfigs[configA], configB => Vector{Float64}([rCutRange[i]])) # Add the configuration to the dictionary
					end
				end

				# Same for configB
				if (similarity > maxSims[configB]) # If new similarity score is greater than previous highest
					similarConfigs[configB] = Dict{Int64,Vector{Float64}}(configA=>Vector{Float64}([rCutRange[i]])) # Create new vector, containing configB
					maxSims[configB] = similarity # Update maximum similarity score
				elseif (similarity == maxSims[configB]) # If the new similarity is the same as the previous highest
					if (haskey(similarConfigs[configB],configA))
						push!(similarConfigs[configB], configA => push!(get!(similarConfigs[configB],configA,0),rCutRange[i]))
					else
						push!(similarConfigs[configB], configA => Vector{Float64}([rCutRange[i]]))
					end
				end
			end
		end
	end
	return similarConfigs, maxSims # Return outputs
end

"""
Converts a 2D array of containing the CNA profile of each configuration, for each atom, to a vector of dictionaries,
mapping atom CNA profiles to their frequency, for each configuration.

Inputs:
atomicProfiles: Array{Dict{String,Int}} 3D array containing the CNA profile of each configuration, for each rCut value, for each atom.
i: (Int64) The index of the current rCut value.
L: (Int64) The number of configurations in the file.

Outputs:
atomicProfileDict: Vector{Dict{String,Int64}} ector of dictionaries, mapping atom CNA profiles to their frequency, for each configuration.
"""
function convertAtomicProfiles(atomicProfiles,i,L)
	# Convert atomicProfiles into format better for computing similarity
	atomicProfilesDict = [Dict{String,Int64}() for j in 1:L] # Initialise new format
	for j in 1:L # For all configurations
		atomicProfile = atomicProfiles[i,j,:] # Get atom profile
		for atomProfile in atomicProfile # For each atom
			# Increments atomic CNA profile frequency by 1
			push!(atomicProfilesDict[j], "$atomProfile" => get!(atomicProfilesDict[j],"$atomProfile",0)+1)
		end
	end

	return atomicProfilesDict # Return output
end

"""
This function returns the similarity score of two configurations based on their CNA profiles.

Inputs:
configA: (Int64) The configuration number of one of the configurations to compare.
configB: (Int64) The configuration number of the other of the configurations to compare.
totalProfiles: Array{Dict{String,Int}} 2D array containg the totalCNA profiles of each configuration, for each rCut value.
atomicProfilesDict: Vector{Dict{String,Int64}} ector of dictionaries, mapping atom CNA profiles to their frequency, for each configuration.
similarityMeasure: (String) The method of comparing configurations, either "atomic" or "total" CNA.
N: (Int64) The number of atoms in each configuration.
i: (Int64) The index of the current rCut value.

Outputs:
similarity: (Float64) The similarity score between the two configurations.
"""
function similarityScore(configA, configB, totalProfiles,atomicProfilesDict,similarityMeasure,N,i)
	similarity = 0.0 # Initialise similarity score
	if (similarityMeasure == "atomic") # If want to compare using atomic CNA profile
		# Compute similarity of atomic CNA profiles
		for key in keys(atomicProfilesDict[configA]) # For each atom CNA profile
			# Add the number of atoms with the same atom CNA profile
			similarity += min(get(atomicProfilesDict[configA],key,0),get(atomicProfilesDict[configB],key,0))
		end
		similarity /= N # Normalise the similarity score

	elseif (similarityMeasure == "total") # If want to compare using total CNA profile
		# Compute Jacard index of their total CNA profiles
		similarity = JacardIndex(totalProfiles[i,configA],totalProfiles[i,configB])
	else # If none of the valid similarity measures were used
		throw(DomainError(similarityMeasure, "Invalid similarity measure")) # Throw an error
	end

	return round(similarity; digits = 3) # Return rounded similarity score
end

"""
This function calculates the Jacard Index for two dictionaries,
with their keys being the unique elements of a multiset and the values are their frequencies

Inputs:
A: (Dictionary) Multiset A
B: (Dictionary) Multiset B

Outputs:
similarity: (Float64) The Jacard Index of multsets A & B
"""
function JacardIndex(A,B)
	intersection = 0.0 # Initialise the number of elements in the intersection of A & B
	union = 0.0 # Initialise the number of elements in the union of A & B
	for key in keys(A) # For all the unique elements of A
		intersection += min(get(A,key,0),get(B,key,0)) # Add the minimum number of occurances of the key across both A & B
		union += get(A,key,0) # Add all the elements of A
	end

	for key in keys(B) # For all the unique elements of B
		union += get(B,key,0) # Add all the elements of B
	end

	union -= intersection # union = A + B - intersection
	similarity = intersection/union # Comput the Jacard index
	return similarity # Return the Jacard Index
end

"""
This function groups configurations into distinguisable sets based on their similarity scores.

Inputs:
energies: (Vector{Float64}) The energy of each configuration in each the file.
sortingArray: Vector{Int64} A vector of configuration numbers, sorted by maximum similarity score.
maxSims: Vector{Float64} A Vector of the maximum similarity score for each configuration when compared to every other configuration
						 over all rCut values, in decreasing order.
similarConfigs: Vector{Dict{Int64,Vector{Float64}}} A vector of dictionaries which map the most similar configurations to the rCut values
						 							for which it is maximised, for each configuraiton, sorted by similarity.
similarityThreshold: (Float64) Similarity threshold above which two configurations are considered identical.
M: (Int64) The number of rCut values tested.

Keyword Arguments:
rCutThreshold: (Float64) Proportion of rCut values that have to meet the similarity threshold for two configs to be considered identical.

Outputs:
uniqueConfigs: Vector{Vector{Int64}} Vector of configuration numbers of the distinguisable sets of configurations.
uniqueEnergies: Vector{Vector{Float64}} Vector of configuration energies of the distinguisable sets of configurations.
uniqueSimilarities: Vector{Float64} Vector of similarity scores of each of the distinguiable sets.
"""
function groupConfigs(energies,sortingArray,maxSims,similarConfigs,similarityThreshold,M;rCutThreshold=0.9)
	# Initialise arrays
	uniqueConfigs = Vector{Vector{Int64}}()
	uniqueEnergies = Vector{Vector{Float64}}()
	uniqueSimilarities = Vector{Float64}()
	classified = Vector{Int64}() # Vector of configuration numbers that have already been grouped
	for (j,config) in enumerate(sortingArray) # For all configurations
		if (!(config in classified)) # If the configuration has not already been grouped
			same = Vector{Int64}([config])
			if (maxSims[j] > similarityThreshold) # If the maximum similarity of the configuration exceeds a threshold
				for key in keys(similarConfigs[j])
					if (length(get!(similarConfigs[j],key,0)) > rCutThreshold*M && !(key in classified))
						# Get vector of configuration numbers of configurations that are similar
						push!(same,key)
					end
				end
			end
			push!(uniqueConfigs,same) # Add the set of configurations to the list
			push!(uniqueEnergies,energies[same]) # Add the energies of the configurations to the list
			push!(uniqueSimilarities,maxSims[j]) # Add the similarity score to the list
			classified = vcat(classified,same) # Add the configurations to the list of already classified
		end
	end

	return uniqueConfigs, uniqueEnergies, uniqueSimilarities # Return outputs
end

"""
This function returns some basic statistics of the energies of each of the distinguisbale configuration sets.

Inputs:
uniqueEnergies: (Vector{Vector{Float64}}) Vector of configuration energies of the distinguisable sets of configurations.
numUnique (Int64) Number of distinguiable sets of configurations.

Outputs:
uniqueEnergyStats: Array{Float64} 2D array of energy statistics. First row is the mean, second row is standard deviation.
"""
function configStats(uniqueEnergies,numUnique)
	uniqueEnergyStats = Array{Float64}(undef,2,numUnique) # Initialse array
	for i in 1:numUnique # For each grouped set of configurations
		uniqueEnergyStats[1,i] = Statistics.mean(uniqueEnergies[i]) # Compute mean
		uniqueEnergyStats[2,i] = Statistics.std(uniqueEnergies[i]) # Compute standard deviation
	end

	return uniqueEnergyStats # Return statistics
end

end # End of Module