text:
#=
  GLM Link Functions
=#

using Distributions

abstract type AbstractMatrixType end
struct RegularData <: AbstractMatrixType end
struct Block1D <: AbstractMatrixType end
struct Block1DParallel <: AbstractMatrixType end

# Defining the standard normal for further use
const standardNormal = Normal()

# Standard Normal pdf
function pdfSN(eta::T)::T where {T <: AbstractFloat}
  return exp(-(eta^2)/2)/sqrt(2*pi)
end

# Compliment of epsilon 1 - eps(Float)
function ceps(T::Type{<: AbstractFloat})
  return 1 - eps(T)
end

#=
  Types for Link Functions
=#
abstract type AbstractLink end
struct IdentityLink <: AbstractLink end
struct LogLink <: AbstractLink end
struct InverseLink <: AbstractLink end
struct NegativeBinomialLink{T <: AbstractFloat} <: AbstractLink
  Alpha::T
end

# abstract type AbstractProbabilityLink <: AbstractLink end
struct LogitLink <: AbstractLink end
struct ProbitLink <: AbstractLink end
struct CauchitLink <: AbstractLink end
struct OddsPowerLink{T <: AbstractFloat} <: AbstractLink
  Alpha::T
end
struct LogComplimentLink <: AbstractLink end
struct LogLogLink <: AbstractLink end
struct ComplementaryLogLogLink <: AbstractLink end
struct PowerLink{T <: AbstractFloat} <: AbstractLink
  Alpha::T
end

#====================================================================
  Copy constructors for link types
====================================================================#
import Base.copy
copy(::IdentityLink) = IdentityLink()
copy(::LogLink) = LogLink()
copy(::InverseLink) = InverseLink()
function copy(distrib::NegativeBinomialLink{T}) where {T}
  return NegativeBinomialLink(distrib.Alpha)
end
copy(::LogitLink) = LogitLink()
copy(::ProbitLink) = ProbitLink()
copy(::CauchitLink) = CauchitLink()
function copy(distrib::OddsPowerLink{T}) where {T}
  return OddsPowerLink(distrib.Alpha)
end
copy(::LogComplimentLink) = LogComplimentLink()
copy(::LogLogLink) = LogLogLink()
copy(::ComplementaryLogLogLink) = ComplementaryLogLogLink()
function copy(distrib::PowerLink{T}) where {T}
  return PowerLink(distrib.Alpha)
end
#===================================================================#

#=
  Implementation of the groups of link functions
=#

# IdentityLink Functions
@inline function linkfun(::IdentityLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return copy(mu)
end
@inline function deta_dmu(::IdentityLink, mu::Array{T, 1}, eta::Array{T})::Array{T, 1} where {T <: AbstractFloat}
  return fill(T(1.0), size(mu))
end
@inline function linkinv(::IdentityLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return copy(eta)
end
# LogLink Functions
@inline function linkfun(::LogLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return log.(mu)
end
@inline function deta_dmu(::LogLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return mu.^(-1)
end
@inline function linkinv(::LogLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return exp.(eta)
end
# InverseLink Functions
@inline function linkfun(::InverseLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return mu.^(-1)
end
@inline function deta_dmu(::InverseLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return -(mu.^-2)
end
@inline function linkinv(::InverseLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return eta.^(-1)
end
# LogitLink Functions
@inline function linkfun(::LogitLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return log.(mu ./(1 .- mu))
end
@inline function deta_dmu(::LogitLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return (mu .* (1 .- mu)).^(-1)
end
@inline function linkinv(::LogitLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  tmp = exp.(eta)
  return tmp ./ (1 .+ tmp)
end
# CauchitLink Functions
@inline function linkfun(::CauchitLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return tan.(pi .* (mu .- (1/2)))
end
@inline function deta_dmu(::CauchitLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return pi .* (cos.(pi .* (mu .- 0.5)).^(-2))
end
@inline function linkinv(::CauchitLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return (atan.(eta)./pi) .+ 0.5
end
# ProbitLink Functions
@inline function linkfun(::ProbitLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return quantile.(standardNormal, mu)
end
@inline function deta_dmu(::ProbitLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return pdf.(standardNormal, eta).^-1
end
@inline function linkinv(::ProbitLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return cdf.(standardNormal, eta)
end
# PowerLink{Alpha} Functions
@inline function linkfun(link::PowerLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  if link.Alpha == T(0)
    return log.(mu)
  end
  return mu .^link.Alpha
end
@inline function deta_dmu(link::PowerLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  if link.Alpha == T(0)
    return mu.^(-1)
  end
  return link.Alpha .* (mu .^(link.Alpha - 1))
end
@inline function linkinv(link::PowerLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  if link.Alpha == T(0)
    return exp.(eta)
  end
  return eta .^(link.Alpha^-1)
end
# OddsPowerLink Functions
#=
  Reference from:
  1. Stata Documentation "rglm - Robust variance estimates for generalized linear models"
  2. Multivariate Methods in Epidemiology by Theodore Holford, Appendix 5 , Table A5-1
=#
# OddsPowerLink Functions
@inline function linkfun(link::OddsPowerLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  if link.Alpha == 0
    return linkfun(LogitLink(), mu)
  end
  return ((mu ./ (1 .- mu)).^link.Alpha .- 1)./link.Alpha
end
@inline function deta_dmu(link::OddsPowerLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  if link.Alpha == 0
    return deta_dmu(LogitLink(), mu, eta)
  end
  return (mu .^ (link.Alpha - 1))./((1 .- mu).^(link.Alpha + 1))
end
@inline function linkinv(link::OddsPowerLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  if link.Alpha == 0
    return linkinv(LogitLink(), eta)
  end
  tmp = ((eta .* link.Alpha .+ 1).^(1/link.Alpha))
  return min.(max.(tmp ./ (1 .+ tmp), eps(T)), ceps(T))
end
# LogComplement Link Functions - needs further work
@inline function linkfun(::LogComplimentLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return log.(max.(1 .- mu, eps(T)))
end
@inline function deta_dmu(::LogComplimentLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return -(max.(1 .- mu, eps(T))).^(-1)
end
@inline function linkinv(::LogComplimentLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return min.(max.(-expm1.(eta), eps(T)), ceps(T))
end
# LogLogLink Link Functions - needs further work
@inline function linkfun(::LogLogLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return -log.(-log.(mu))
end
@inline function deta_dmu(::LogLogLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return -1 ./(mu .* log.(mu))
end
@inline function linkinv(::LogLogLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return exp.(-exp.(-eta))
end
# ComplementaryLogLogLink Functions
@inline function linkfun(::ComplementaryLogLogLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return log.(-log.(1 .- mu))
end
@inline function deta_dmu(::ComplementaryLogLogLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return ((mu .- 1) .* log.(1 .- mu)).^-1
end
@inline function linkinv(::ComplementaryLogLogLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return 1 .- exp.(-exp.(eta))
end
# NegativeBinomialLink Functions - rewrite with 1/Alpha
@inline function linkfun(link::NegativeBinomialLink, mu::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  iAlpha = T(1)/link.Alpha
  return log.(mu ./(mu .+ iAlpha))
end
@inline function deta_dmu(link::NegativeBinomialLink, mu::Array{T, 1}, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  return 1 ./ (mu .+ link.Alpha .* mu .^2)
end
@inline function linkinv(link::NegativeBinomialLink, eta::Array{T, 1})::Array{T, 1} where {T <: AbstractFloat}
  tmp = exp.(eta)
  return tmp ./(link.Alpha .* (1 .- tmp))
end

# For block algorithm
function linkfun(link::AbstractLink, mu::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [linkfun(link, mu[i]) for i in 1:nBlocks]
end
function deta_dmu(link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [deta_dmu(link, mu[i], eta[i]) for i in 1:nBlocks]
end
function linkinv(link::AbstractLink, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(eta)
  return [linkinv(link, eta[i]) for i in 1:nBlocks]
end

# For parallel block algorithm
function linkfun(::Block1DParallel, link::AbstractLink, mu::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = linkfun(link, mu[i])
  end
  return ret
end
function deta_dmu(::Block1DParallel, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = deta_dmu(link, mu[i], eta[i])
  end
  return ret
end
function linkinv(::Block1DParallel, link::AbstractLink, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(eta)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = linkinv(link, eta[i])
  end
  return ret
end

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
# This file was generated, do not modify it. # hide
LinMod = @pipeline(FeatureSelector(features=[:hp1]),
                   LinearRegressor());

ext: jl


text:
# we should find sole marker on infinite field without borders and so

# count value ++ when its bigger geworden,  we move in digit horizonside
function find_marker(r::Robot)
    count = 1
    while !ismarker(r)
        for i = 0:1
            move_side(r, count, HorizonSide(i))
        end
        count += 1
        for i = 2:3
            move_side(r, count, HorizonSide(i))
        end
        count += 1;
    end
end

# move in side count times
function move_side(r::Robot, count::Int, side::HorizonSide)
    for i = 1:count
        if ismarker(r)
            break
        end
        move!(r, side)
    end
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
#
#  The following functions allow us to give printed feedback to the user
#
#	pretty_print_list(myList)
# 	warning(msg::String)
#

function pretty_print_list(myList; orient="vertical", digits=3, enum=false)
	#
	if orient=="horizontal"
		for ii = 1:length(myList)
			i = myList[ii]
			if typeof(i)<:Number
				i = round(i, digits=digits)
			end
			if enum
				print(ii, ". ")
			end
	    	print(i)
	    	print("  ")
	    end
	else
	    for ii = 1:length(myList)
	    	i = myList[ii]
	    	if typeof(i)<:Number
				i = round(i, digits=digits)
			end
			if enum
				print(ii, ". ")
			end
	    	println(i)
	    end
    end
end;
function warning(msg::String)
    red = "\033[1m\033[31m"
    println("********************************************************************************************")
    println(string("     WARNING: ", msg))
    println("********************************************************************************************")
    return
end
function progressbar(iter,total)
    done = ["=", "=", "=", "=", "=", "=", "=", "=", "=", "="]
    incomplete = ["-", "-", "-", "-", "-", "-", "-", "-", "-", "-"]
    if mod(iter,total*0.1) == 0
        ndone = round(Int,iter/total * 10);
        nincomp = round(Int, (1 - iter/total) * 10);
        println("   *", join(done[1:ndone]), join(incomplete[1:nincomp]), " (", iter, "/", total, ") ",Dates.format(now(), "mm/dd/yy HH:MM:ss") )
    elseif iter==1
        ndone = 0;
        nincomp = 10;
        println("   *", join(done[1:ndone]), join(incomplete[1:nincomp]), " (", iter, "/", total, ") ",Dates.format(now(), "mm/dd/yy HH:MM:ss") )
    end
end;


text:
using DataFrames, Gadfly
import Compose: w, h

set_default_plot_size(12inch,3.3inch)

# Issue #1357

df = DataFrame(x=[0.6w, 0.6w], y=[0.6h, 0.7h], label=["Text 1", "Text 2"])
lyr1 = layer(x=1.0:5, y=[0.21, 0.81, 0.68, 0.76, 0.18], Geom.point)
lyr2 = layer(df, x=:x, y=:y, label=:label, Geom.label(position=:right))
lyr3 = layer(xintercept=[0.5inch], Geom.vline(color="orange"))

p1 = plot(xmin=[.25], xmax=[.75], Geom.band, color=[colorant"red"])
p2 = plot(xmin=[.25], xmax=[.75], Geom.band, color=[colorant"red"],
    layer(Geom.rect, xmin=[0], xmax=[1], ymin=[0], ymax=[1], color=[colorant"blue"]))
p3 = plot(lyr1, lyr2, lyr3)
p4 = plot(lyr1, lyr2, lyr3, Scale.x_discrete)

hstack(p1, p2, p3, p4)


text:
module MicroCoverage

include("types.jl")

include("public_interface.jl")

include("assert.jl")
include("clean.jl")
include("coverage_compute.jl")
include("coverage_write.jl")
include("instrument.jl")
include("preview_coverage.jl")
include("read_write_julia_files.jl")
include("restore.jl")
include("start_tracking_file.jl")
include("tracked_files.jl")
include("tracker.jl")
include("utils.jl")

end # end module MicroCoverage


text:
using SQLite

dir_elements = joinpath(Pkg.dir("Alpine"), "src/elements")

const db = SQLite.DB(joinpath(dir_elements, "elements.db"))

function element_properties()
    df = SQLite.columns(db, "elements")
    keys = SQLite.columns(db, "elements")[:name]
    vals = SQLite.columns(db, "elements")[:type]
    Dict(zip(keys, vals))
end

function element_property(symbol::Symbol, property::Symbol)
    df = SQLite.query(db, "SELECT * FROM elements WHERE symbol == ?;", values = [string(symbol)])
    df[property][1]
end

max_stars_repo_name: NickMcNutt/Alpine

max_issues_repo_name: NickMcNutt/Alpine

max_forks_repo_name: NickMcNutt/Alpine


text:
####################
####### DEMO #######
####################

## import required packages ##
using Distributions
import LightGraphs
import Gadfly
import Plots
Plots.gr()
using LaTeXStrings

include("rw_smc.jl")
include("utils.jl")
include("rand_rw.jl")
include("gibbs_updates.jl")

# set seed
srand(0)

########################
######### DATA #########
########################

######### SYNTHETIC #########
    syn_data = true
    # set synthetic data parameters
    n_edges_data = 25 # number of edges
    α = 0.25 # parameter for new vertex probability
    λ = 4.0 # paramter for random walk length distribution
    ld = Poisson(λ) # random walk length distribution (automatically shifted by +1 if minimum is zero)
    sb = true # degree-biased distribution for first vertex at each step (false for uniform distribution)

    # sample simple graph (returns edge list)
    g = randomWalkSimpleGraph(n_edges=n_edges_data,alpha_prob=α,length_distribution=ld,sizeBias=sb)

    # convert edge list to adjacency matrix
    A = edgelist2adj(g)
    Gadfly.spy(A)

######### REAL DATA #########
    syn_data = false
    # SJS
    # g = readcsv("data/sjsAdj.csv",Int64)
    # NIPS
    g = readcsv("data/nipsAdj_02_03.csv",Int64)
    # PPI
    # g = readcsv("data/interactomeAdj.csv",Int64)

    n_edges_data = size(g,1)
    # convert edge list to adjacency matrix
    A = edgelist2adj(g)
    Gadfly.spy(A)

######################################
# INTERFACE WITH LIGHTGRAPHS PACKAGE #
######################################

# create LightGraphs.Graph object for use with LightGraphs functions
    lg = LightGraphs.Graph(A)

    # Examples:
    lg_conn = LightGraphs.is_connected(lg)
    lg_diam = LightGraphs.diameter(lg)

#######################################
#### RUN SMC WITH FIXED PARAMETERS ####
#######################################

n_particles = 100

α_fixed = α
λ_fixed = λ

particle_paths,log_marginal = SMC(g,n_particles,α_fixed,λ_fixed,sb)

############################
#### RUN PARTICLE GIBBS ####
############################

# set sampler parameters
n_particles = 100
n_mcmc_iter = 500 # total number of iterations; includes burn-in
n_burn = 100 # burn-in
n_collect = 10 # collect a sample every n_collect iterations
k_trunc = 10*lg_diam # truncation for support of latent K distribution

# set starting values for sampling (can be anything)
α_start = convert(Float64,0.5*(maximum(g)-1)/(n_edges_data-1))
λ_start = convert(Float64,sqrt(lg_diam))

# set prior parameters
a_α = 1.0
b_α = 1.0
a_λ = 1.0
b_λ = 0.25

n_print = 1 # print progress updates every n_print iterations

samples = ParticleGibbs(g,n_particles,
                  a_α,b_α,a_λ,b_λ,
                  α_start,λ_start,sb,k_trunc,
                  n_mcmc_iter,n_burn,n_collect,n_print)

# plot samples
Plots.scatter(samples["lambda"],samples["alpha"],legend=false,
          ylims=[0,1],xlims=[0,15],
          markershape=:circle,markersize=5.5,markeralpha=0.7,markercolor=:grey,markerstrokewidth=0.0,
          xlabel=L"\lambda",ylabel=L"\alpha")
if syn_data
    Plots.scatter!([λ],[α],legend=false,marker=(:star4,12.0,1.0,:black))
end

# plot edge sequence
Plots.plot(samples["edge_sequence"]',legend=false,xlabel="Sampled arrival order",ylabel="Edge index")


text:
tocursors(range, start=1) =
  cursor(range[1]-start+1, range[2]+1),
  cursor(range[1]-start+1, range[3]+1)

function fillranges(code, values, ranges, line = 1)
  reader = LineNumberingReader(code)
  result = IOBuffer()
  for i = 1:length(ranges)
    start, stop = tocursors(ranges[i], line)
    while cursor(reader) < start
      print(result, read(reader, Char))
    end
    while cursor(reader) < stop
      read(reader, Char)
    end
    print(result, values[i])
  end
  while !eof(reader)
    print(result, read(reader, Char))
  end
  return takebuf_string(result)
end


text:
## fix for heatmap visualization
function spatial_array(x::AbstractArray{T,2}) where T <: Real
    flipud(transpose(x))
end

flipud(x::AbstractArray{T, 2}) where T <: Real = x[end:-1:1, :]
flipud(x::AbstractArray{T, 3}) where T <: Real = x[end:-1:1, :, :]

fliplr(x::AbstractArray{T, 2}) where T <: Real = x[:, end:-1:1]
fliplr(x::AbstractArray{T, 3}) where T <: Real = x[:, end:-1:1, :]

"""Function to flip GeoArray upside down to adjust to GDAL ecosystem."""
function flipud!(ga::GeoArray)
    # Flip data upside down
    ga.A = reverse(ga.A, dims=2)

    # Find new corner coordinates
    ux, uy = ga.f(SVector{2}([0, size(ga)[2]]))

    # Define y mirror and compose
    lm = LinearMap(SMatrix{2,2}([1.0 0.0; 0.0 -1.0]))
    am = compose(ga.f, lm)  # AffineMap
    translate = SVector{2}([ux, uy])
    f = AffineMap(am.linear, translate)

    ga.f = f
    ga
end

export fliplr, flipud, flipud!


text:
module TestLearningCurves

using Test
using Distributed

@everywhere begin
    using MLJBase
    using MLJTuning
    using ..Models
    import ComputationalResources: CPU1, CPUProcesses, CPUThreads
    #using Distributed
end
using Random
Random.seed!(1234*myid())
using ..TestUtilities

x1 = rand(100);
x2 = rand(100);
x3 = rand(100);
X = (x1=x1, x2=x2, x3=x3);
y = 2*x1 .+ 5*x2 .- 3*x3 .+ 0.2*rand(100);

@testset_accelerated "learning curves" accel begin
    atom = FooBarRegressor()
    ensemble = EnsembleModel(atom=atom, n=50)
    mach = machine(ensemble, X, y)
    r_lambda = range(ensemble, :(atom.lambda),
                     lower=0.0001, upper=0.1, scale=:log10)
    printstyled("\n Testing progressmeter rngs option with $(accel) and CPU1 grid \n", color=:bold)
    if accel == CPU1()
        curve = @test_logs((:info, r"No measure"),
                           (:info, r"Training"),
                           (:info, r"Attempting"),
                           learning_curve(mach; range=r_lambda,
                                          acceleration=accel))
    else
        curve = learning_curve(mach; range=r_lambda,
                               acceleration=accel)
    end
    @test curve isa NamedTuple{(:parameter_name,
                                :parameter_scale,
                                :parameter_values,
                                :measurements)}
    @test length(curve.parameter_values) == length(curve.measurements)
    atom.lambda=0.3
    r_n = range(ensemble, :n, lower=10, upper=100)
    curves = learning_curve(mach; range=r_n, resolution=7,
                             acceleration=accel,
                             rngs = MersenneTwister.(1:3),
                             rng_name=:rng, verbosity = 1)
    @test size(curves.measurements) == (length(curves.parameter_values), 3)
    @test length(curves.parameter_values) == 7

    # individual curves are different:
    @test !(curves.measurements[1,1] ≈ curves.measurements[1,2])
    @test !(curves.measurements[1,1] ≈ curves.measurements[1,3])

    # reproducibility:
    curves2 = learning_curve(mach; range=r_n, resolution=7,
                             acceleration=accel,
                             rngs = 3,
                             rng_name=:rng, verbosity=0)
    @test curves2.measurements ≈ curves.measurements

    # alternative signature:
    curves3 = learning_curve(ensemble, X, y; range=r_n, resolution=7,
                             acceleration=accel,
                             rngs = 3,
                             rng_name=:rng, verbosity=0)

    @test curves2.measurements ≈ curves3.measurements

   # restricting rows gives different answer:
    curves4 = learning_curve(ensemble, X, y; range=r_n, resolution=7,
                             rows = 1:60,
                             acceleration=accel,
                             rngs = 3,
                             rng_name=:rng, verbosity=0)

    @test !(curves4.measurements[1] ≈ curves2.measurements[1])

end

@static if VERSION >= v"1.3.0-DEV.573"

@testset_accelerated "learning curves (CPUThreads grid) " accel begin
    atom = FooBarRegressor()
    ensemble = EnsembleModel(atom=atom, n=50)
    mach = machine(ensemble, X, y)
    r_lambda = range(ensemble, :(atom.lambda),
                     lower=0.0001, upper=0.1, scale=:log10)
    printstyled("\n Testing progressmeter rngs option with $(accel) and CPUThreads grid \n", color=:bold)
    if accel == CPU1()
        curve = @test_logs((:info, r"No measure"),
                           (:info, r"Training"),
                           (:info, r"Attempting"),
                           learning_curve(mach; range=r_lambda,
                                          acceleration=accel, acceleration_grid = CPUThreads()))
    else
        curve = learning_curve(mach; range=r_lambda,
                               acceleration=accel, acceleration_grid = CPUThreads())
    end
    @test curve isa NamedTuple{(:parameter_name,
                                :parameter_scale,
                                :parameter_values,
                                :measurements)}
    @test length(curve.parameter_values) == length(curve.measurements)
    atom.lambda=0.3
    r_n = range(ensemble, :n, lower=10, upper=100)
    curves = learning_curve(mach; range=r_n, resolution=7,
                             acceleration=accel, acceleration_grid = CPUThreads(),
                             rngs = MersenneTwister.(1:3),
                             rng_name=:rng, verbosity = 1)
    @test size(curves.measurements) == (length(curves.parameter_values), 3)
    @test length(curves.parameter_values) == 7

    # individual curves are different:
    @test !(curves.measurements[1,1] ≈ curves.measurements[1,2])
    @test !(curves.measurements[1,1] ≈ curves.measurements[1,3])

    # reproducibility:
    curves2 = learning_curve(mach; range=r_n, resolution=7,
                             acceleration=accel, acceleration_grid = CPUThreads(),
                             rngs = 3,
                             rng_name=:rng, verbosity=0)
    @test curves2.measurements ≈ curves.measurements

    # alternative signature:
    curves3 = learning_curve(ensemble, X, y; range=r_n, resolution=7,
                             acceleration=accel, acceleration_grid = CPUThreads(),
                             rngs = 3,
                             rng_name=:rng, verbosity=0)

    @test curves2.measurements ≈ curves3.measurements

    # restricting rows gives different answer:
    curves4 = learning_curve(ensemble, X, y; range=r_n, resolution=7,
                             acceleration=accel, acceleration_grid = CPUThreads(),
                             rows = 1:60,
                             rngs = 3,
                             rng_name=:rng, verbosity=0)

    @test !(curves4.measurements[1] ≈ curves2.measurements[1])

end
end

@testset_accelerated "learning curves (CPUProcesses grid) " accel begin
    atom = FooBarRegressor()
    ensemble = EnsembleModel(atom=atom, n=50)
    mach = machine(ensemble, X, y)
    r_lambda = range(ensemble, :(atom.lambda),
                     lower=0.0001, upper=0.1, scale=:log10)
    printstyled("\n Testing progressmeter rngs option with $(accel) and CPUProcesses grid \n", color=:bold)
    if accel == CPU1()
        curve = @test_logs((:info, r"No measure"),
                           (:info, r"Training"),
                           (:info, r"Attempting"),
                           learning_curve(mach; range=r_lambda,
                                          acceleration=accel, acceleration_grid = CPUProcesses()))
    else
        curve = learning_curve(mach; range=r_lambda,
                               acceleration=accel, acceleration_grid = CPUProcesses())
    end
    @test curve isa NamedTuple{(:parameter_name,
                                :parameter_scale,
                                :parameter_values,
                                :measurements)}
    @test length(curve.parameter_values) == length(curve.measurements)
    atom.lambda=0.3
    r_n = range(ensemble, :n, lower=10, upper=100)
    curves = learning_curve(mach; range=r_n, resolution=7,
                             acceleration=accel, acceleration_grid = CPUProcesses(),
                             rngs = MersenneTwister.(1:3),
                             rng_name=:rng, verbosity = 1)
    @test size(curves.measurements) == (length(curves.parameter_values), 3)
    @test length(curves.parameter_values) == 7

    # individual curves are different:
    @test !(curves.measurements[1,1] ≈ curves.measurements[1,2])
    @test !(curves.measurements[1,1] ≈ curves.measurements[1,3])

    # reproducibility:
    curves2 = learning_curve(mach; range=r_n, resolution=7,
                             acceleration=accel, acceleration_grid = CPUProcesses(),
                             rngs = 3,
                             rng_name=:rng, verbosity=0)
    @test curves2.measurements ≈ curves.measurements

    # alternative signature:
    curves3 = learning_curve(ensemble, X, y; range=r_n, resolution=7,
                             acceleration=accel, acceleration_grid = CPUProcesses(),
                             rngs = 3,
                             rng_name=:rng, verbosity=0)

    @test curves2.measurements ≈ curves3.measurements

    # restricting rows gives different answer:
    curves4 = learning_curve(ensemble, X, y; range=r_n, resolution=7,
                             acceleration=accel, acceleration_grid = CPUProcesses(),
                             rows=1:60,
                             rngs = 3,
                             rng_name=:rng, verbosity=0)

    @test !(curves4.measurements[1] ≈ curves2.measurements[1])

end

end # module
true


text:
"""
A julia wrapper and interface for both XFOIL and airfoilprep.py, and other
functions for airfoil processing.
"""
module AirfoilPrep

# ------------ GENERIC MODULES -------------------------------------------------
using PyCall
using PyPlot
using JLD
using Dierckx
using Interpolations
using Roots
using LaTeXStrings
using Statistics: mean
import CSV
import DataFrames

# ------------ FLOW CODES ------------------------------------------------------
# Xfoil from https://github.com/byuflowlab/Xfoil.jl
import Xfoil

# ------------ GLOBAL VARIABLES ------------------------------------------------
const module_path = splitdir(@__FILE__)[1]          # Path to this module

# ------------ LOAD airfoilprep.py ---------------------------------------------
path_airfoilpreppy = module_path                    # Path to airfoilprep.py
prepy = PyNULL()                                    # airfoilpreppy module

function __init__()
    imp = pyimport("imp")
    (file, filename, data) = imp.find_module("airfoilprep", [path_airfoilpreppy])
    copy!(prepy, imp.load_module("airfoilprep", file, filename, data))
end

# ------------ HEADERS ---------------------------------------------------------
for header_name in ["pywrapper", "xfoil", "ndtools", "liftprops", "misc"]
    include("AirfoilPrep_"*header_name*".jl")
end

end # END OF MODULE


text:
@using LinearAlgebraicRepresentation as LAR
#@using LinearAlgebra as LA
@using ViewerGL as GL
using SparseArrays
include("../../src/Lar.jl")
LarA = Lar.Arrangement

function build_copFE(FV::Lar.Cells, EV::Lar.Cells)
	copFE = LAR.u_coboundary_1(FV, EV) # unsigned
	faceedges = [findnz(copFE[f,:])[1] for f=1:size(copFE,1)]

	f_edgepairs = Array{Array{Int64,1}}[]
	for f=1:size(copFE,1)
		edgepairs = Array{Int64,1}[]
		for v in FV[f]
			push!(edgepairs, [e for e in faceedges[f] if v in EV[e]])
		end
		push!(f_edgepairs, edgepairs)
	end
	copFE.=0
	for f=1:size(copFE,1)
		facepairs = sort(sort.(f_edgepairs[f]))
		# setting a reference
		copFE[f, facepairs[1][1]] = 1
		for (e1,e2) in facepairs
			@assert copFE[f, e1] != 0
			if copFE[f, e2] == 0
				copFE[f, e2] =
					-((EV[e1][1]==EV[e2][1])||(EV[e1][2]==EV[e2][2])) * copFE[f,e1]+
					((EV[e1][1]==EV[e2][2])||(EV[e1][2]==EV[e2][1])) * copFE[f,e1]
			end
		end
	end
	return copFE
end

function getGFVModel(V::Lar.Points, FV::Array{Array{Int,1}})::Lar.Model
	function getFaceEdges(fV::Array{Int,1})
		return [fV[1], fV[2]], [fV[1], fV[3]], [fV[2], fV[3]]
	end
	model = Lar.Model(V);
	EV = unique([(map(f -> getFaceEdges(f), FV)...)...])
	Lar.addModelCells!(model, 1, LAR.build_copEV(EV))
	#FE = [[3x+1, 3x+2, 3x+3] for x = 0 : length(FV)-1]
	Lar.addModelCells!(model, 2, build_copFE(FV, EV))
	return Lar.mergeModelVertices(model)
end

scaling = 1.2
no_sphere = 40
expl = 2.0

carry = Array{Tuple{Array{Float64,2},Array{Array{Int64,1},1}}}(undef, no_sphere)
for i = 1 : no_sphere
	carry[i] = LAR.struct2lar(LAR.Struct([
		LAR.t(rand(0.5:1e-2:5, 3)...),
		LAR.sphere(rand(1:1e-2:3) * scaling)([2,3])
	]))
end
V, FV = LAR.struct2lar(LAR.Struct(carry))
FV = sort(sort.(FV))
GL.VIEW([GL.GLGrid(V, FV)])

m = getGFVModel(V, FV)
Lar.viewModel(m, 1)
Lar.viewModel(m, 2)
GL.VIEW(GL.GLExplode(m.G,[[s] for s in FV],1.2,1.2,1.2,99,1));

sp_idx = LarA.spaceIndex(m, 2)

de_models = [
	LarA.face_decomposition(m, face_idx, sp_idx[face_idx])
	for face_idx = 1 : size(m, 2, 1)
]

mfrag = Lar.mergeMultipleModels(de_models)

triangulated_faces = LAR.triangulate2D(convert(Lar.Points, mfrag.G'), mfrag.T[1:2])
FVs = convert(Array{Lar.Cells}, triangulated_faces);
GL.VIEW( GL.GLExplode(mfrag.G,FVs,expl,expl,expl,99,1) );

mArranged = Lar.deepcopy(mfrag)
Lar.addModelCells!(mArranged, 3, LAR.Arrangement.minimal_3cycles(
		convert(Lar.Points, mfrag.G'), mfrag.T[1], mfrag.T[2])
	)

V,CVs,FVs,EVs = LAR.pols2tria(mArranged.G, mArranged.T[1], mArranged.T[2], mArranged.T[3])

GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));
meshes = GL.GLExplode(V,CVs[1:end],8,4,6,99,0.5);
push!( meshes, GL.GLFrame)
GL.VIEW( meshes );

#-------------------------------------------------------------------------------
#
# V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(
# 	 W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp);
#
#
# triangulated_faces = Lar.triangulate2D(V, [copEV, copFE]);
# FVs = convert(Array{Lar.Cells}, triangulated_faces);
# V = convert(Lar.Points, V');
# GL.VIEW( GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1) );
#
#
# EVs = Lar.FV2EVs(copEV, copFE); # polygonal face fragments
# GL.VIEW( GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1) );


text:
function remove_forecasts_by_type!(system::PowerSystems.System, ::Type{T}) where T <: PowerSystems.Component

    forecasts_to_delete = Vector()
    for first_step in PowerSystems.get_forecast_initial_times(uc_system)
        for f in PowerSystems.get_component_forecasts(T, uc_system, first_step)
            push!(forecasts_to_delete, f)
        end
    end

    for f in forecasts_to_delete
        PowerSystems.remove_forecast!(uc_system, f)
    end

    return
end

function explore_cost_functions(sys::PSY.System)
    p = plot(legend = :bottomright)
    for t in PowerSystems.get_components(PowerSystems.ThermalStandard, sys)
        plot!([x[2] for x in t.op_cost.variable.cost], [x[1] for x in t.op_cost.variable.cost], label = t.name)
    end
    return(p)
end


text:
# Line geometry connects (x, y) coordinates with lines.
immutable LineGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement

    # Do not reorder points along the x-axis.
    preserve_order::Bool

    function LineGeometry(default_statistic=Gadfly.Stat.identity();
                          preserve_order=false)
        new(default_statistic, preserve_order)
    end
end

const line = LineGeometry

# Only allowing identity statistic in paths b/c I don't think any
# any of the others will work with `preserve_order=true` right now
function path()
    return LineGeometry(preserve_order=true)
end

function density()
    return LineGeometry(Gadfly.Stat.density())
end

function smooth(; method::Symbol=:loess, smoothing::Float64=0.75)
    return LineGeometry(Gadfly.Stat.smooth(method=method, smoothing=smoothing))
end

function step(; direction::Symbol=:hv)
    return LineGeometry(Gadfly.Stat.step(direction=direction))
end

function default_statistic(geom::LineGeometry)
    return geom.default_statistic
end

function element_aesthetics(::LineGeometry)
    return [:x, :y, :color]
end

# Render line geometry.
#
# Args:
#   geom: line geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::LineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.line", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.line", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    aes = inherit(aes, default_aes)

    if length(aes.color) == 1 &&
        !(isa(aes.color, PooledDataArray) && length(levels(aes.color)) > 1)
        points = {(x, y) for (x, y) in zip(aes.x, aes.y)}
        if !geom.preserve_order
            sort!(points, by=first)
        end
        form = compose(lines(points...),
                       stroke(aes.color[1]),
                       svgclass("geometry"))
    else
        # group points by color
        points = Dict{ColorValue, Array{(Any, Any),1}}()
        for (x, y, c) in zip(aes.x, aes.y, cycle(aes.color))
            if !haskey(points, c)
                points[c] = Array((Any, Any),0)
            end
            push!(points[c], (x, y))
        end

        forms = Array(Any, length(points))
        for (i, (c, c_points)) in enumerate(points)
            if !geom.preserve_order
                sort!(c_points, by=first)
            end
            forms[i] =
                compose(lines({(x, y) for (x, y) in c_points}...),
                        stroke(c),
                        svgclass(@sprintf("geometry color_%s",
                                          escape_id(aes.color_label([c])[1]))))
        end
        form = combine(forms...)
    end

    compose(
        canvas(units_inherited=true, order=2),
        form, fill(nothing), linewidth(theme.line_width))
end


text:
include("mean0.jl")

function sd0(a)
    # create variable to store ouput
    out = 0.0

    #mean
    m = mean0(a)

    # sum of square
    sq = 0.0
    for i in a
        sq += (i - m)^2
    end

    # n-1
    n = length(a) - 1

    # calculation
    out = sqrt(sq/n)

    #return
    return out
end


text:
"""
Parameters for the AST version of the `MCTS.DPWSolver` solver.
"""
function MCTSPWSolver(; kwargs...)
    return MCTS.DPWSolver(; estimate_value=AST.rollout, # required.
                            enable_state_pw=false, # required.
                            reset_callback=AST.go_to_state, # Custom fork of MCTS.jl
                            tree_in_info=true,
                            show_progress=true,
                            kwargs...)
end

"""
Return optimal action path from MCTS tree (using `info[:tree]` from `(, info) = action_info(...)`).
"""
function get_optimal_path(mdp, tree, snode::Int, actions::Vector{ASTAction}; verbose::Bool=false)
    best_Q = -Inf
    sanode = 0
    for child in tree.children[snode]
        if tree.q[child] > best_Q
            best_Q = tree.q[child]
            sanode = child
        end
    end

    if verbose
        print("State = 0x", string(tree.s_labels[snode].hash, base=16), "\t:\t")
    end
    if sanode != 0
        if verbose
            print("Q = ", tree.q[sanode], "\t:\t")
            println("Action = ", string(tree.a_labels[sanode]))
        end
        push!(actions, tree.a_labels[sanode])

        # Find subsequent maximizing state node
        best_Q_a = -Inf
        snode2 = 0
        for tran in tree.transitions[sanode]
            if tran[2] > best_Q_a
                best_Q_a = tran[2]
                snode2 = tran[1]
            end
        end

        if snode2 != 0
            get_optimal_path(mdp, tree, snode2, actions, verbose=verbose)
        end
    else
        AST.go_to_state(mdp, tree.s_labels[snode])
        if verbose
            if BlackBox.isevent(mdp.sim)
                println("Event.")
            else
                println("End of tree.")
            end
        end
    end

    return actions::Vector{ASTAction}
end
get_optimal_path(mdp, tree, state, actions::Vector{ASTAction}=ASTAction[]; kwargs...) = get_optimal_path(mdp, tree, tree.s_lookup[state], actions; kwargs...)

"""
    AST.search!(planner::DPWPlanner; return_tree::Bool=false, verbose=false)

Search for failure events using the planner of the `MCTSPWSolver`.

This is the main entry function to get a failure trajectories from the planner.
"""
function AST.search!(planner::DPWPlanner; return_tree::Bool=false, verbose=false)
    mdp::ASTMDP = planner.mdp
    Random.seed!(mdp.params.seed) # Determinism

    initstate = AST.initialstate(mdp)
    tree = MCTS.action_info(planner, initstate, tree_in_info=true)[2][:tree] # this runs MCTS
    action_path::Vector = get_optimal_path(mdp, tree, initstate, verbose=verbose)

    return return_tree ? tree : action_path
end


text:
#########
# TRACE #
#########

function update!(
    tr::OptimizationTrace,
    state::S,
    iteration::Integer,
    population,
    val,
    minimizer,
    dt::Dict,
    store_trace::Bool,
    show_trace::Bool,
    show_every::Int = 1,
    callback = nothing,
) where {S<:AbstractOptimizerState,T,O}
    os = OptimizationTraceRecord(
        iteration,
        copy(population),
        val,
        minimizer,
        dt,
    )
    store_trace && push!(tr, os)
    if show_trace
        if iteration % show_every == 0
            show(os)
            print("\n")
            flush(stdout)
        end
    end
    if callback !== nothing && (iteration % show_every == 0)
        if store_trace
            stopped = callback(tr)
        else
            stopped = callback(os)
        end
    else
        stopped = false
    end
    stopped
end

function trace!(
    tr,
    iteration,
    objfun,
    state,
    method,
    options,
    curr_time = time(),
)
    dt = Dict{String,Any}()
    dt["time"] = curr_time
    # set additional trace value
    trace!(dt, objfun, state, method, options)
    update!(
        tr,
        state,
        iteration,
        state.fittest,
        value(state),
        minimizer(state),
        dt,
        options.store_trace,
        options.show_trace,
        options.show_every,
        options.callback,
    )
end

"""
    trace!(record::Dict{String,Any}, objfun, state, population, method, options)

Update the trace `record`. This function allows to supplement an additional information into the optimization algorithm trace by modifying a trace `record`. It can be overridden by specifying particular parameter types.
"""
trace!(record::Dict{String,Any}, objfun, state, method, options) = ()

########
# MISC #
########

default_values(x::AbstractArray{T}) where {T} = fill!(similar(x), zero(T))
default_values(x::AbstractArray{T}) where {T<:AbstractFloat} = fill!(similar(x), T(NaN))

function funargnum(f)
    fobj = try
        first(methods(f))
    catch
        @error "Cannot find function $f"
    end
    fobj.nargs
end

function vswap!(v1::T, v2::T, idx::Int) where {T<:AbstractVector}
    val = v1[idx]
    v1[idx] = v2[idx]
    v2[idx] = val
end

function dump_result(result::OptimizationResult, directory)
    mkpath(directory)
    result_file = joinpath(directory, "result.toml")
    save(joinpath(directory, "trace.jld"), "data", result.trace)
    if length(result.metrics) > 0
        save(joinpath(directory, "convergence.jld"), "data", result.metrics)
    end
    open(result_file, "w") do io
        data = Dict(
            "algorithm" => string(result.method),
            "iterations" => result.iterations,
            "converged" => result.converged,
            "trace" => "trace.jld",
            "time_run" => time_run(result),
            "time_limit" => time_limit(result),
            "f_calls" => f_calls(result),
            "minimum" => minimum(result),
            "minimizer" => minimizer(result),
        )
        if length(result.metrics) > 0
            data["convergence"] = "convergence.jld"
        end
        TOML.print(io, data)
    end
end

function draw_convergence(result::OptimizationResult, directory)
    images_directory = joinpath(directory, "convergence")
    mkpath(images_directory)

end

function draw_generations(result::OptimizationResult, directory=pwd())
    images_directory = joinpath(directory, "tmp", "generations")
    mkpath(images_directory)
    gen = 0
    anim = @animate for state in result.trace
        gen += 1
        println("Plotting generation $(gen) / $(length(result.trace))")
        for i in state.population
            if length(variables(i)) > 3 || length(variables(i)) < 2
                error("Can't draw individual with dimensions of $(length(variables(i)))")
            end
        end
        fig = Plots.plot(
            [variables(i)[1] for i in state.population],
            [variables(i)[2] for i in state.population],
            seriestype = :scatter,
            title = "Generation $(gen)",
        )
        xlims!((-3, 3))
        ylims!((-3, 3))
        # contour!(-3:0.1:3, -3:0.1:3, (x, y) -> x^2+y^2, legend = false, colorbar = false)
        savefig(fig, joinpath(images_directory, "$(gen).png"))
    end
    gif(anim, joinpath(directory, "generations.gif"), fps = 15, loop = 1)
end

function draw_pareto_fronts(result::OptimizationResult, directory=pwd())
    images_directory = joinpath(directory, "tmp", "objectives")
    mkpath(images_directory)
    gen = 0
    anim = @animate for state in result.trace[2:end]
        gen += 1
        # println("Plotting generation $(gen) / $(length(result.trace))")
        mF = maxF(Population(state.population))
        F = calcF(Population(state.population))
        cols = distinguishable_colors(mF, [RGB(1,1,1), RGB(0,0,0)], dropseed=false)
        fig = nothing
        for i = 1:length(F)
            f = F[i]
            if length(f) == 0
                continue
            end
            if i == 1
                fig = Plots.plot(
                    objectives(Population(state.population[f]))[:, 1],
                    objectives(Population(state.population[f]))[:, 2],
                    seriestype = :scatter,
                    color = cols[i],
                    title = "Generation $(gen)",
                    legend = false,
                )
            else
                Plots.plot!(
                    objectives(Population(state.population[f]))[:, 1],
                    objectives(Population(state.population[f]))[:, 2],
                    seriestype = :scatter,
                    color = cols[i]
                )
            end
        end
        # xlims!((0, 5))
        # ylims!((0, 5))
        savefig(fig, joinpath(images_directory, "$(gen).png"))
    end
    gif(anim, joinpath(images_directory, "objectives.gif"), fps = 20, loop = 1)
end

function draw_objectives(result::OptimizationResult, directory=pwd())
    images_directory = joinpath(directory, "tmp", "objectives")
    mkpath(images_directory)
    gen = 0
    anim = @animate for state in result.trace[2:end]
        gen += 1
        # println("Plotting generation $(gen) / $(length(result.trace))")
        fig = Plots.plot(
            objectives(Population(pfront(state)))[:, 1],
            objectives(Population(pfront(state)))[:, 2],
            seriestype = :scatter,
            title = "Generation $(gen)",
            legend = false,
        )
        savefig(fig, joinpath(images_directory, "$(gen).png"))
    end
    gif(anim, joinpath(images_directory, "objectives.gif"), fps = 20, loop = 1)
end

function draw_metrics(result::OptimizationResult)
    plot(
        1:length(result.trace),
        [state.metrics[2].value for state in result.trace],
        legend = false,
    )
end

function draw_igd(result, truepf)
    Plots.plot(
        1:length(result.trace),
        [state_igd(state, truepf) for state in result.trace],
        legend = false,
    )
end

function state_igd(state, truepf)
    igd(pfront(state), truepf)
end


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 1ef22dae-5c88-11eb-1276-231730d144f6
begin
	using Images
	using Random
	using PlutoUI
	using Statistics
	using Rotations
	using Serialization

	gd = "C:/Users/dusty/Documents/PlaymatProjects/PlaymatGames/Mage"

	USER_SETTINGS = deserialize("$gd/tmp/user_selection.jls")

	DECKS_DIR = "C:/Users/dusty/Documents/PlaymatProjects/PlaymatGames/Mage/EDH/decks"

	md"""
	# Carve custom faces for your cards!
	"""
end

# ╔═╡ 9d2e6708-88ea-11eb-3852-5b09b6960506
md"""
Mage.jl Deck directory path: $(@bind DECK_DIR Select(DECKS_DIR |> readdir))
"""

# ╔═╡ 9d987b58-88e9-11eb-1b8b-e733e3b39391
DECK_NAME = split(DECK_DIR, '/')[end]

# ╔═╡ c0a76bea-88e9-11eb-144e-65dd27acc30b
begin
	custom_imgs_dir = "$DECKS_DIR/$DECK_NAME/custom_images/"
	custom_img_names = [ fn for fn in readdir(custom_imgs_dir) if occursin("png", fn) ]
	push!(custom_img_names, "Backside.png")
	[ try load("$custom_imgs_dir/$img_fn") catch end for img_fn in custom_img_names ]
end

# ╔═╡ 8823075c-88e9-11eb-3aa1-e9a9af738720
deck = deserialize("$DECKS_DIR/$DECK_DIR/$DECK_NAME.jls")

# ╔═╡ b978978c-5ee0-11eb-1443-33825d99e295
md"""
#### $(@bind carve_custom_card CheckBox()) Carve cards in custom images dir $(@bind custom_img_select Slider(1:length(custom_img_names), show_value=true))
"""

# ╔═╡ c4bb34e4-5ede-11eb-0871-e16410f0c048
if carve_custom_card
	md"""Shrink custom card $(@bind custom_img_ratio Slider(0:0.01:1, default=0.5, show_value=true)) by percentage or $(@bind specify_custom_size CheckBox(default=true)) specify card width $(@bind custom_width NumberField(200:10:300, default=270))"""
end

# ╔═╡ b92971e6-5d70-11eb-1713-7327a0bca668
md"""
#### $(@bind carve_deck_card CheckBox()) Carve cards in deck $(@bind deck_img_select Slider(1:length(deck[:CARD_FACES]), show_value=true))
"""

# ╔═╡ 6ff626d4-5d72-11eb-110e-49b0464e95a7
if carve_deck_card
	deck_img = deck[:CARD_FACES][deck_img_select][end][end]
	save("blot_img.png", deck_img)
	deck_img
end

# ╔═╡ 501cb384-5d73-11eb-20cb-d3b96fc930d0
begin
	custom_img = load("$custom_imgs_dir/$(custom_img_names[custom_img_select])")
	imgs = carve_deck_card ? [ deck_img ] : [ custom_img ]
end;

# ╔═╡ 4c43b0ba-5ed7-11eb-3fa3-7d9b1ba5e434
if carve_custom_card
	card_ratio = height(custom_img) / width(custom_img)
	imgs[end] = specify_custom_size ? imresize(custom_img, ceil(Int, custom_width * card_ratio), custom_width) : imresize(custom_img, ratio=custom_img_ratio)
end

# ╔═╡ 4ecd1258-5d11-11eb-1b53-71740a55b68d
md"""
##### Vertically shrink image up to : $(@bind vert_shrink_by NumberField(0:200, default=160)) pixels $(@bind shrink_greedy_v CheckBox())
"""

# ╔═╡ 72c7f7b6-5c8e-11eb-26a8-61e826cea075
md"""
##### Horizontally shrink image up to : $(@bind hor_shrink_by NumberField(0:40, default=20)) pixels $(@bind shrink_greedy_h CheckBox())
"""

# ╔═╡ 4a156ba2-5d1b-11eb-07a5-730b3f65033f
md"""
 **$(@bind export_img CheckBox()) Export card image to $DECK_DIR?**
"""

# ╔═╡ 148ab572-5d15-11eb-2a54-89eaf1d95ca4
if export_img
	save("$DECK_DIR/carved_$(randstring(5)).png", imgs[end])
	if carve_deck_card
		deck[:CARD_FACES][deck_img_select] = deck[:CARD_FACES][deck_img_select][begin] => [ imgs[end] ]
	end
end

# ╔═╡ 26ba7e48-5ecb-11eb-22f1-b5013af4ec24
md"""
**Functions**
"""

# ╔═╡ deffb654-5c8d-11eb-3138-a5c7634793bc
begin
	brightness(c::RGB) = mean((c.r, c.g, c.b))
	brightness(c::RGBA) = mean((c.r, c.g, c.b))
end

# ╔═╡ e5dad786-5c8d-11eb-1da6-e5fd55842485
convolve(img, k) = imfilter(img, reflect(k))

# ╔═╡ f334877c-5ce4-11eb-1d51-83452cede3d8
energy(∇x, ∇y) = sqrt.(∇x.^2 .+ ∇y.^2)

# ╔═╡ ed9f0b4c-5c8d-11eb-0bda-4d2fe4942325
function energy(img)
	∇y = convolve(brightness.(img), Kernel.sobel()[1])
	∇x = convolve(brightness.(img), Kernel.sobel()[2])
	energy(∇x, ∇y)
end

# ╔═╡ 88d7b806-5c88-11eb-220a-9bae6e26693b
function mark_path(img, path)
	img′ = copy(img)
	m = size(img, 2)
	for (i, j) in enumerate(path)
		for j′ in j-1:j+1
			img′[i, clamp(j′, 1, m)] = RGB(1,0,1)
		end
	end
	img′
end

# ╔═╡ 93f5643a-5ba3-11eb-3bfa-3f1e453ec23b
function remove_in_each_row_views(img, column_numbers)
	@assert size(img, 1) == length(column_numbers)
	m, n = size(img)
	local img′ = similar(img, m, n-1)

	for (i, j) in enumerate(column_numbers)
		img′[i, 1:j-1] .= @view img[i, 1:j-1]
		img′[i, j:end] .= @view img[i, j+1:end]
	end

	img′
end

# ╔═╡ a1b7250c-5ba3-11eb-0e28-2b8e62c2007a
function shrink_n(img, n, min_seam, imgs=[]; show_lightning=false)
	n==0 && return push!(imgs, img)
	e = energy(img)

	seam_energy(seam) = sum(e[i, seam[i]] for i in 1:size(img, 1))
	_, min_j = findmin(map(j->seam_energy(min_seam(e, j)), 1:size(e, 2)))
	min_seam_vec = min_seam(e, min_j)
	img′ = remove_in_each_row_views(img, min_seam_vec)

	if show_lightning
		push!(imgs, mark_path(img, min_seam_vec))
	else
		push!(imgs, img′)
	end

	shrink_n(img′, n-1, min_seam, imgs)
end

# ╔═╡ 26500512-5d00-11eb-37a8-458570a946fd
function shrink!(img, min_seam, cols=true)
	img = cols ? img : rotr90(img)
	e = energy(img)

	seam_energy(seam) = sum(e[i, seam[i]] for i in 1:size(img, 1))
	_, min_j = findmin(map(j->seam_energy(min_seam(e, j)), 1:size(e, 2)))
	min_seam_vec = min_seam(e, min_j)

	img = remove_in_each_row_views(img, min_seam_vec)
	img = cols ? img : rotl90(img)

	return img
end

# ╔═╡ a7b310bc-5ba3-11eb-38e5-3d025712d57e
function greedy_seam(energies, starting_pixel::Int)
	is = [ starting_pixel ]
	m, n = size(energies)

	for k in 2:m
		es = energies[k, clamp(is[end]-1,1,n):clamp(is[end]+1,1,n)]
		push!(is, clamp(is[end] + argmin(es) - clamp(is[end],0,2), 1, n))
	end

	return is
end

# ╔═╡ d1dd1162-5d10-11eb-3d7f-8f10be30efa7
if shrink_greedy_v
	greedy_carved_v = shrink_n(rotr90(imgs[end]), vert_shrink_by, greedy_seam, show_lightning=false)
	md"""Shrink by: $(@bind greedy_v Slider(1:vert_shrink_by, show_value=true)) Save changes to img? $(@bind save_v CheckBox())"""
end

# ╔═╡ 38bc0512-5ee4-11eb-1576-1fac86b5e112
if shrink_greedy_v
	size(rotl90(greedy_carved_v[greedy_v]))
end

# ╔═╡ e061aca2-5d12-11eb-3661-59430c9f0ce4
if shrink_greedy_v && (@isdefined greedy_carved_v)
	rotl90(greedy_carved_v[greedy_v])
end

# ╔═╡ 0d2b79d4-5d13-11eb-03c5-f90edc2114e0
if (@isdefined save_v) && save_v
	imgs[end] = deepcopy(rotl90(greedy_carved_v[greedy_v]))
	"Changes saved!"
end

# ╔═╡ 3fe076c8-5c8e-11eb-3222-ad52ddd23256
if shrink_greedy_h
	greedy_carved_h = shrink_n(imgs[end], hor_shrink_by, greedy_seam, show_lightning=false)

	md"""
	Shrink by: $(@bind greedy_h Slider(1:hor_shrink_by, show_value=true))
	Save changes to img? $(@bind save_h CheckBox())
	"""
end

# ╔═╡ e1c05af0-5ee4-11eb-153f-af5abf6bd5a9
if shrink_greedy_h
	size(greedy_carved_h[greedy_h])
end

# ╔═╡ a3c7f5f8-5c8e-11eb-3807-290228dc3ea8
if shrink_greedy_h && (@isdefined greedy_carved_h)
	greedy_carved_h[greedy_h]
end

# ╔═╡ b4d3252a-5d10-11eb-16c2-1330d81f5f7c
if (@isdefined save_h) && save_h
	imgs[end]=deepcopy(greedy_carved_h[greedy_h])
	"Changes saved!"
end

# ╔═╡ 2c06e12e-5d7d-11eb-0da4-c71de86166e8
if ((@isdefined save_v) && save_v) || ((@isdefined save_h) && save_h)
	serialize("$DECKS_DIR/$DECK_DIR/$DECK_NAME.jls", deck)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
Images = "~0.24.1"
PlutoUI = "~0.7.9"
Rotations = "~1.0.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[ArrayInterface]]
deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "045ff5e1bc8c6fb1ecb28694abba0a0d55b5f4f5"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.17"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[AxisAlgorithms]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9"
uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
version = "1.0.0"

[[AxisArrays]]
deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"]
git-tree-sha1 = "f31f50712cbdf40ee8287f0443b57503e34122ef"
uuid = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
version = "0.4.3"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[CatIndices]]
deps = ["CustomUnitRanges", "OffsetArrays"]
git-tree-sha1 = "a0f80a09780eed9b1d106a1bf62041c2efc995bc"
uuid = "aafaddc9-749c-510e-ac4f-586e18779b91"
version = "0.2.2"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "dcc25ff085cf548bc8befad5ce048391a7c07d40"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "0.10.11"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "42a9b08d3f2f951c9b283ea427d96ed9f1f30343"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.5"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.31.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[ComputationalResources]]
git-tree-sha1 = "52cb3ec90e8a8bea0e62e275ba577ad0f74821f7"
uuid = "ed09eef8-17a6-5b46-8889-db040fac31e3"
version = "0.3.2"

[[CoordinateTransformations]]
deps = ["LinearAlgebra", "StaticArrays"]
git-tree-sha1 = "6d1c23e740a586955645500bbec662476204a52c"
uuid = "150eb455-5306-5404-9cee-2592286d6298"
version = "0.6.1"

[[CustomUnitRanges]]
git-tree-sha1 = "537c988076d001469093945f3bd0b300b8d3a7f3"
uuid = "dc8bdbbb-1ca9-579f-8c36-e416f6a65cce"
version = "1.0.1"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distances]]
deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
git-tree-sha1 = "abe4ad222b26af3337262b8afb28fab8d215e9f8"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.3"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EllipsisNotation]]
deps = ["ArrayInterface"]
git-tree-sha1 = "8041575f021cba5a099a456b4163c9a08b566a02"
uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
version = "1.1.0"

[[FFTViews]]
deps = ["CustomUnitRanges", "FFTW"]
git-tree-sha1 = "70a0cfd9b1c86b0209e38fbfe6d8231fd606eeaf"
uuid = "4f61f5a4-77b1-5117-aa51-3ab5ef4ef0cd"
version = "0.3.1"

[[FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "f985af3b9f4e278b1d24434cbb546d6092fca661"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.4.3"

[[FFTW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3676abafff7e4ff07bbd2c42b3d8201f31653dcc"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.9+8"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.10.1"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
git-tree-sha1 = "2c1cf4df419938ece72de17f368a021ee162762e"
uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
version = "1.1.0"

[[IdentityRanges]]
deps = ["OffsetArrays"]
git-tree-sha1 = "be8fcd695c4da16a1d6d0cd213cb88090a150e3b"
uuid = "bbac6d45-d8f3-5730-bfe4-7a449cd117ca"
version = "0.3.1"

[[IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.0"

[[ImageAxes]]
deps = ["AxisArrays", "ImageCore", "Reexport", "SimpleTraits"]
git-tree-sha1 = "794ad1d922c432082bc1aaa9fa8ffbd1fe74e621"
uuid = "2803e5a7-5153-5ecf-9a86-9b4c37f5f5ac"
version = "0.6.9"

[[ImageContrastAdjustment]]
deps = ["ColorVectorSpace", "ImageCore", "ImageTransformations", "Parameters"]
git-tree-sha1 = "2e6084db6cccab11fe0bc3e4130bd3d117092ed9"
uuid = "f332f351-ec65-5f6a-b3d1-319c6670881a"
version = "0.3.7"

[[ImageCore]]
deps = ["AbstractFFTs", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
git-tree-sha1 = "db645f20b59f060d8cfae696bc9538d13fd86416"
uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
version = "0.8.22"

[[ImageDistances]]
deps = ["ColorVectorSpace", "Distances", "ImageCore", "ImageMorphology", "LinearAlgebra", "Statistics"]
git-tree-sha1 = "6378c34a3c3a216235210d19b9f495ecfff2f85f"
uuid = "51556ac3-7006-55f5-8cb3-34580c88182d"
version = "0.2.13"

[[ImageFiltering]]
deps = ["CatIndices", "ColorVectorSpace", "ComputationalResources", "DataStructures", "FFTViews", "FFTW", "ImageCore", "LinearAlgebra", "OffsetArrays", "Requires", "SparseArrays", "StaticArrays", "Statistics", "TiledIteration"]
git-tree-sha1 = "bf96839133212d3eff4a1c3a80c57abc7cfbf0ce"
uuid = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
version = "0.6.21"

[[ImageIO]]
deps = ["FileIO", "Netpbm", "PNGFiles", "TiffImages", "UUIDs"]
git-tree-sha1 = "d067570b4d4870a942b19d9ceacaea4fb39b69a1"
uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
version = "0.5.6"

[[ImageMagick]]
deps = ["FileIO", "ImageCore", "ImageMagick_jll", "InteractiveUtils", "Libdl", "Pkg", "Random"]
git-tree-sha1 = "5bc1cb62e0c5f1005868358db0692c994c3a13c6"
uuid = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
version = "1.2.1"

[[ImageMagick_jll]]
deps = ["JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "1c0a2295cca535fabaf2029062912591e9b61987"
uuid = "c73af94c-d91f-53ed-93a7-00f77d67a9d7"
version = "6.9.10-12+3"

[[ImageMetadata]]
deps = ["AxisArrays", "ColorVectorSpace", "ImageAxes", "ImageCore", "IndirectArrays"]
git-tree-sha1 = "ae76038347dc4edcdb06b541595268fca65b6a42"
uuid = "bc367c6b-8a6b-528e-b4bd-a4b897500b49"
version = "0.9.5"

[[ImageMorphology]]
deps = ["ColorVectorSpace", "ImageCore", "LinearAlgebra", "TiledIteration"]
git-tree-sha1 = "68e7cbcd7dfaa3c2f74b0a8ab3066f5de8f2b71d"
uuid = "787d08f9-d448-5407-9aad-5290dd7ab264"
version = "0.2.11"

[[ImageQualityIndexes]]
deps = ["ColorVectorSpace", "ImageCore", "ImageDistances", "ImageFiltering", "OffsetArrays", "Statistics"]
git-tree-sha1 = "1198f85fa2481a3bb94bf937495ba1916f12b533"
uuid = "2996bd0c-7a13-11e9-2da2-2f5ce47296a9"
version = "0.2.2"

[[ImageShow]]
deps = ["Base64", "FileIO", "ImageCore", "OffsetArrays", "Requires", "StackViews"]
git-tree-sha1 = "832abfd709fa436a562db47fd8e81377f72b01f9"
uuid = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
version = "0.3.1"

[[ImageTransformations]]
deps = ["AxisAlgorithms", "ColorVectorSpace", "CoordinateTransformations", "IdentityRanges", "ImageCore", "Interpolations", "OffsetArrays", "Rotations", "StaticArrays"]
git-tree-sha1 = "d966631de06f36c8cd4bec4bb2e8fa731db16ed9"
uuid = "02fcd773-0e25-5acc-982a-7f6622650795"
version = "0.8.12"

[[Images]]
deps = ["AxisArrays", "Base64", "ColorVectorSpace", "FileIO", "Graphics", "ImageAxes", "ImageContrastAdjustment", "ImageCore", "ImageDistances", "ImageFiltering", "ImageIO", "ImageMagick", "ImageMetadata", "ImageMorphology", "ImageQualityIndexes", "ImageShow", "ImageTransformations", "IndirectArrays", "OffsetArrays", "Random", "Reexport", "SparseArrays", "StaticArrays", "Statistics", "StatsBase", "TiledIteration"]
git-tree-sha1 = "8b714d5e11c91a0d945717430ec20f9251af4bd2"
uuid = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
version = "0.24.1"

[[IndirectArrays]]
git-tree-sha1 = "c2a145a145dc03a7620af1444e0264ef907bd44f"
uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
version = "0.5.1"

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[IntelOpenMP_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2018.0.3+2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[Interpolations]]
deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "1470c80592cf1f0a35566ee5e93c5f8221ebc33a"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.13.3"

[[IntervalSets]]
deps = ["Dates", "EllipsisNotation", "Statistics"]
git-tree-sha1 = "3cc368af3f110a767ac786560045dceddfc16758"
uuid = "8197267c-284f-5f27-9208-e0e47529a953"
version = "0.5.3"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "LinearAlgebra"]
git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.2.5"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "c253236b0ed414624b083e6b72bfe891fbd2c7af"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2021.1.1+1"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "6a8a2a625ab0dea913aba95c11370589e0239ff0"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.6"

[[MappedArrays]]
git-tree-sha1 = "18d3584eebc861e311a552cbb67723af8edff5de"
uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
version = "0.4.0"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MosaicViews]]
deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d"
uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
version = "0.3.3"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[Netpbm]]
deps = ["ColorVectorSpace", "FileIO", "ImageCore"]
git-tree-sha1 = "09589171688f0039f13ebe0fdcc7288f50228b52"
uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
version = "1.0.1"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "2bf78c5fd7fa56d2bbf1efbadd45c1b8789e6f57"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.2"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PNGFiles]]
deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
git-tree-sha1 = "520e28d4026d16dcf7b8c8140a3041f0e20a9ca8"
uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
version = "0.3.7"

[[PaddedViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "0fa5e78929aebc3f6b56e1a88cf505bb00a354c4"
uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
version = "0.5.8"

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.2"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RangeArrays]]
git-tree-sha1 = "b9039e93773ddcfc828f12aadf7115b4b4d225f5"
uuid = "b3c3ace0-ae52-54e7-9d0b-2c1406fd6b9d"
version = "0.3.2"

[[Ratios]]
git-tree-sha1 = "37d210f612d70f3f7d57d488cb3b6eff56ad4e41"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.0"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[Rotations]]
deps = ["LinearAlgebra", "StaticArrays", "Statistics"]
git-tree-sha1 = "2ed8d8a16d703f900168822d83699b8c3c1a5cd8"
uuid = "6038ab10-8711-5258-84ad-4b1120ba62dc"
version = "1.0.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a50550fa3164a8c46747e62063b4d774ac1bcf49"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.5.1"

[[StackViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
version = "0.1.1"

[[Static]]
deps = ["IfElse"]
git-tree-sha1 = "2740ea27b66a41f9d213561a04573da5d3823d4b"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.2.5"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "a43a7b58a6e7dc933b2fa2e0ca653ccf8bb8fd0e"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.6"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2f6792d523d7448bbe2fec99eca9218f06cc746d"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.8"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[TensorCore]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
version = "0.1.1"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TiffImages]]
deps = ["ColorTypes", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "OffsetArrays", "OrderedCollections", "PkgVersion", "ProgressMeter"]
git-tree-sha1 = "03fb246ac6e6b7cb7abac3b3302447d55b43270e"
uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
version = "0.4.1"

[[TiledIteration]]
deps = ["OffsetArrays"]
git-tree-sha1 = "52c5f816857bfb3291c7d25420b1f4aca0a74d18"
uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
version = "0.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[WoodburyMatrices]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "59e2ad8fd1591ea019a5259bd012d7aee15f995c"
uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
version = "0.5.3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─1ef22dae-5c88-11eb-1276-231730d144f6
# ╟─9d2e6708-88ea-11eb-3852-5b09b6960506
# ╟─9d987b58-88e9-11eb-1b8b-e733e3b39391
# ╟─c0a76bea-88e9-11eb-144e-65dd27acc30b
# ╟─8823075c-88e9-11eb-3aa1-e9a9af738720
# ╟─501cb384-5d73-11eb-20cb-d3b96fc930d0
# ╟─b978978c-5ee0-11eb-1443-33825d99e295
# ╟─c4bb34e4-5ede-11eb-0871-e16410f0c048
# ╟─4c43b0ba-5ed7-11eb-3fa3-7d9b1ba5e434
# ╟─b92971e6-5d70-11eb-1713-7327a0bca668
# ╟─6ff626d4-5d72-11eb-110e-49b0464e95a7
# ╟─4ecd1258-5d11-11eb-1b53-71740a55b68d
# ╟─38bc0512-5ee4-11eb-1576-1fac86b5e112
# ╟─d1dd1162-5d10-11eb-3d7f-8f10be30efa7
# ╟─e061aca2-5d12-11eb-3661-59430c9f0ce4
# ╟─0d2b79d4-5d13-11eb-03c5-f90edc2114e0
# ╟─72c7f7b6-5c8e-11eb-26a8-61e826cea075
# ╟─3fe076c8-5c8e-11eb-3222-ad52ddd23256
# ╟─e1c05af0-5ee4-11eb-153f-af5abf6bd5a9
# ╟─a3c7f5f8-5c8e-11eb-3807-290228dc3ea8
# ╟─b4d3252a-5d10-11eb-16c2-1330d81f5f7c
# ╟─4a156ba2-5d1b-11eb-07a5-730b3f65033f
# ╟─148ab572-5d15-11eb-2a54-89eaf1d95ca4
# ╟─2c06e12e-5d7d-11eb-0da4-c71de86166e8
# ╟─26ba7e48-5ecb-11eb-22f1-b5013af4ec24
# ╟─deffb654-5c8d-11eb-3138-a5c7634793bc
# ╟─e5dad786-5c8d-11eb-1da6-e5fd55842485
# ╟─f334877c-5ce4-11eb-1d51-83452cede3d8
# ╟─ed9f0b4c-5c8d-11eb-0bda-4d2fe4942325
# ╟─88d7b806-5c88-11eb-220a-9bae6e26693b
# ╟─93f5643a-5ba3-11eb-3bfa-3f1e453ec23b
# ╟─a1b7250c-5ba3-11eb-0e28-2b8e62c2007a
# ╟─26500512-5d00-11eb-37a8-458570a946fd
# ╟─a7b310bc-5ba3-11eb-38e5-3d025712d57e
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: dustyirwin/PlaymatSimulator

max_issues_repo_name: dustyirwin/PlaymatSimulator

max_forks_repo_name: dustyirwin/PlaymatSimulator


text:
include("magpoly.jl")

using .magpoly
const mp = magpoly
using DelimitedFiles, Random, Distributions, LinearAlgebra

function generate_saw(n_samples::Int, samples_lag::Int)
    io = open("features.txt", "r")
    features = readdlm(io,Float64)
    close(io)
    n_mono = size(features, 1)

    burn_in = ceil(Int, 20*n_mono^(1.11))
    #burn_in = 3
    pol = mp.Magnetic_polymer(n_mono, 1.0, 0.0, 1.0)
    mp.initialize_poly!(pol)

    for i_mono in 1:n_mono
        for j in 1:3
            pol.trial_coord[j,i_mono] = pol.coord[j,i_mono]
        end
    end

    hash_base = 2*n_mono +1

    # equilibration of the chain
    for i in 1:burn_in
        pivot = rand(2:n_mono-1)
        p_move = rand(1:47)
        mp.try_pivot!(pivot, p_move, pol.coord, pol.trial_coord, pol.n_mono)
        still_saw = mp.checksaw!(pivot, pol, hash_base)
        #println(still_saw)
        if still_saw
            for i_mono in 1:n_mono
                for j in 1:3
                    pol.coord[j,i_mono] = pol.trial_coord[j,i_mono]
                end
            end
        end
        empty!(pol.hash_saw)
    end
    println("End of burn in!")
    n_acc = 0

    n_steps = n_samples*samples_lag

    saw_configs = zeros(Int, n_samples, 3, n_mono)

    for i in 1:n_steps
        pivot = rand(2:n_mono-1)
        p_move = rand(1:47)
        if i%10000 == 0
            println("Pivot: ", pivot)
            println("Attmpted move: ", p_move)
            println("step number: ", i, "\n")
        end
        mp.try_pivot!(pivot, p_move, pol.coord, pol.trial_coord, pol.n_mono)
        still_saw = mp.checksaw!(pivot, pol, hash_base)
        if still_saw
            for i_mono in 1:n_mono
                for j in 1:3
                    pol.coord[j,i_mono] = pol.trial_coord[j,i_mono]
                end
            end
            n_acc += 1
        end

        if i%samples_lag == 0
            i_sample = div(i,samples_lag)
            for i_mono in 1:n_mono
                for j in 1:3
                    saw_configs[i_sample,j,i_mono] = pol.coord[j,i_mono]
                end
            end
        end
        empty!(pol.hash_saw)
    end
    println("Fraction of accepted pivots: ", n_acc/n_steps)

    !isdir("SAW_data") && mkdir("SAW_data")
    io = open("SAW_data/saw_configs.txt", "w")
    for i_sample in 1:n_samples
        writedlm(io, saw_configs[i_sample,:,:])
    end
    close(io)
end

# This function generates samples from the posterior distribution of the weights
# Assuming that the data were generated with zero spins coupling i.e. a logistic regression problem
# The sampling is done with a simple metroplis hastings sampler.
function vanilla_sampler(n_samples::Int, delta_w::Float64, initial_weights::Array{Float64}, features_file::String, data_file::String)
    io = open(features_file, "r")
    features = readdlm(io,Float64)
    close(io)
    n_mono = size(features, 1)
    n_feats = size(features, 2)
    println("Number of features: " ,n_feats)

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

    weights_series = zeros(n_feats, n_samples)
    energies = zeros(n_samples)

    weights = zeros(n_feats)
    trial_weights = zeros(n_feats)
    weights .= initial_weights
    weights_series[:,1] .= weights

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

    energy = 0.0
    for i_data in 1:n_data
        for i_mono in 1:n_mono
            energy += data_spins[i_data, i_mono]*fields[i_mono] - log(1+exp(fields[i_mono]))
        end
    end
    energies[1] = energy

    w_acceptance = 0.0
    n_acc = 0

    for i_sample in 2:n_samples
        trial_weights .= weights .+ randn(n_feats) .* delta_w
        fields = zeros(n_mono)
        for i in 1:n_feats
            fields .+= features[:,i] .* trial_weights[i]
        end

        #trial_energy = -0.5*(trial_weights[1]^2 + trial_weights[2]^2)/(w_std^2)
        trial_energy = 0.0
        for i_data in 1:n_data
            for i_mono in 1:n_mono
                trial_energy += data_spins[i_data, i_mono]*fields[i_mono] - log(1+exp(fields[i_mono]))
            end
        end

        delta_energy = trial_energy - energy
        delta_energy >=0 ? w_acceptance=1 : w_acceptance=exp(delta_energy)

        if i_sample%1000 == 0
            println(i_sample)
            for i_feat in 1:n_feats
                println("w$(i_feat): ",weights[i_feat]," ---> ",trial_weights[i_feat])
            end
            println("delta_energy: ", delta_energy)
            println("acceptance: ", w_acceptance)
        end

        if rand() <= w_acceptance
            energy = trial_energy
            weights .= trial_weights
            n_acc += 1
        end

        energies[i_sample] = energy
        weights_series[:, i_sample] .= weights
    end

    println("Acceptance fraction: ", n_acc/(n_samples-1))

    !isdir("IS_data") && mkdir("IS_data")
    open("IS_data/weights$(n_data)_$(delta_w).txt", "w") do io
        writedlm(io, weights_series)
    end
    open("IS_data/neg_energies$(n_data)_$(delta_w).txt", "w") do io
        writedlm(io, energies)
    end
end

function generate_spin_confs(features_file::String, n_data::Int, weights::Array{Float64})
    io = open(features_file, "r")
    features = readdlm(io,Float64)
    close(io)
    n_mono = size(features, 1)
    n_feats = size(features, 2)
    fields = zeros(n_mono)
    for i in 1:n_feats
        fields .+= features[:,i] .* weights[i]
    end
    proba = zeros(n_mono)
    for i in 1:n_mono
        proba[i] = exp(fields[i]) / (1 + exp(fields[i]))
    end
    spins_conf = zeros(Int, n_mono, n_data)
    for i_data in 1:n_data
        for i_mono in 1:n_mono
            if rand() <= proba[i_mono]
                spins_conf[i_mono, i_data] = 1
            end
        end
    end

    open("LR_saw_conf_data_1.txt", "w") do io
        writedlm(io,transpose([weights]))
        writedlm(io,transpose(spins_conf))
    end
end


text:
# make this into test script, spline proxy vs lorentzian.

include("../src/NMRSpectraSimulator.jl")
import .NMRSpectraSimulator

using LinearAlgebra
using FFTW
import PyPlot
import JSON

#import Clustering
import Statistics

import Random
Random.seed!(25)

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

### user inputs.

tol_coherence = 1e-2 # resonances are pairs of eigenvalues of the Hamiltonian that have quantum numbers that differ by -1. This is the tolerance away from -1 that is allowed.
α_relative_threshold = 0.05 # resonances with relative amplitude less than this factor compared to the maximum resonance in the spin group will be removed. Set to 0.0 to see every single resonance component.
Δc_partition_radius = 0.3 # determines how many resonances get grouped together. Larger number means less groups and thus more resonances per group.
λ0 = 3.4

Δr_default = 1.0 # the samples used to build the surrogate is taken every `Δr` radian on the frequency axis. Decrease for improved accuracy at the expense of computation resources.
Δκ_λ_default = 0.05 # the samples used to build thes urrogate for κ_λ are taken at this sampling spacing. Decrease for improved accuracy at the expense of computation resources.
Δcs_max_scalar_default = 0.2 # In units of ppm. interpolation border that is added to the lowest and highest resonance frequency component of the mixture being simulated.
κ_λ_lb_default = 0.5 # interpolation lower limit for κ_λ.
κ_λ_ub_default = 2.5 # interpolation upper limit for κ_λ.

SH_config_path = "/home/roy/Documents/repo/NMRData/input/SH_configs/select_compounds_SH_configs.json"
surrogate_config_path = "/home/roy/Documents/repo/NMRData/input/surrogate_configs/select_compounds_SH_configs.json"

molecule_names = ["L-Serine"; "L-Phenylalanine"; "DSS"; "Ethanol"; "L-Isoleucine"]
#molecule_names = ["D-(+)-Glucose"; "DSS"]

# machine values taken from the BMRB 700 MHz 20 mM glucose experiment.
fs = 14005.602240896402
SW = 20.0041938620844
ν_0ppm = 10656.011933076665

# # machine values for the BMRB 500 MHz glucose experiment.
# ν_0ppm = 6752.490995937095
# SW = 16.0196917451925
# fs = 9615.38461538462

# path to the json file that provides the mapping from a compound name to its spin system info file name.
H_params_path = "/home/roy/Documents/repo/NMRData/input/coupling_info"
dict_compound_to_filename = JSON.parsefile("/home/roy/Documents/repo/NMRData/input/compound_mapping/select_compounds.json")

### end inputs.

hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)

println("Timing: mag equivalence")
@time MEs = NMRSpectraSimulator.getmageqinfomixture(molecule_names,
    H_params_path,
    dict_compound_to_filename;
    unique_cs_atol = 1e-6)

println("Timing: setupmixtureproxies()")
@time mixture_params = NMRSpectraSimulator.setupmixtureSH(molecule_names,
    H_params_path, dict_compound_to_filename, fs, SW,
    ν_0ppm;
    MEs = MEs,
    config_path = SH_config_path,
    tol_coherence = tol_coherence,
    α_relative_threshold = α_relative_threshold,
    Δc_partition_radius = Δc_partition_radius)
As = mixture_params

dummy_SSFID = NMRSpectraSimulator.SpinSysParamsType1(0.0)
u_min = ppm2hzfunc(-0.5)
u_max = ppm2hzfunc(4.0)

Bs = NMRSpectraSimulator.fitproxies(As, dummy_SSFID, λ0;
    names = molecule_names,
    config_path = surrogate_config_path,
    Δcs_max_scalar_default = Δcs_max_scalar_default,
    κ_λ_lb_default = κ_λ_lb_default,
    κ_λ_ub_default = κ_λ_ub_default,
    u_min = u_min,
    u_max = u_max,
    Δr_default = Δr_default,
    Δκ_λ_default = Δκ_λ_default)

### plot.

# purposely distort the spectra by assigning random values to model parameters.
B = Bs[1]
B.ss_params.d[:] = rand(length(B.ss_params.d))
B.ss_params.κs_λ[:] = rand(length(B.ss_params.κs_λ)) .+ 1
B.ss_params.κs_β[:] = collect( rand(length(B.ss_params.κs_β[i])) .* (2*π) for i = 1:length(B.ss_params.κs_β) )

f = uu->NMRSpectraSimulator.evalmixture(uu, mixture_params, Bs)

# test params.
#ΩS_ppm = collect( hz2ppmfunc.( NMRSpectraSimulator.combinevectors(A.Ωs) ./ (2*π) ) for A in mixture_params )
#ΩS_ppm_flat = NMRSpectraSimulator.combinevectors(ΩS_ppm)

P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
U = ppm2hzfunc.(P)
U_rad = U .* (2*π)

## parameters that affect qs.
# A.d, A.κs_λ, A.κs_β
# A.d_singlets, A.αs_singlets, A.Ωs_singlets, A.β_singlets, A.λ0, A.κs_λ_singlets
q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, mixture_params, Bs)

Es = collect( NMRSpectraSimulator.καFIDModelType(Bs[i]) for i = 1:length(Bs) )
q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, mixture_params, Es)

f_U = f.(U_rad)
q_U = q.(U_rad)

discrepancy = abs.(f_U-q_U)
max_val, ind = findmax(discrepancy)
println("relative discrepancy = ", norm(discrepancy)/norm(f_U))
println("max discrepancy: ", max_val)
println()

# ## remove areas with low signal from plotting to reduce plot size.
# reduction_factor = 100
# threshold_factor =  α_relative_threshold/10
# inds, keep_inds, inds_for_reducing = NMRSpectraSimulator.prunelowsignalentries(q_U, threshold_factor, reduction_factor)
#
# q_U_display = q_U[inds]
# f_U_display = f_U[inds]
# P_display = P[inds]

q_U_display = q_U
f_U_display = f_U
P_display = P

## visualize.
PyPlot.figure(fig_num)
fig_num += 1

# PyPlot.plot(P, real.(f_U), label = "f")
# PyPlot.plot(P, real.(q_U), label = "q")
PyPlot.plot(P_display, real.(f_U_display), label = "f")
PyPlot.plot(P_display, real.(q_U_display), label = "q")
PyPlot.plot(P_display, real.(q_U_display), "x")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("f vs q")

## visualize. zoom in.

inds = findall(xx->(2.5<xx<3.9), P_display)

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P_display[inds], real.(f_U_display[inds]), label = "f")
PyPlot.plot(P_display[inds], real.(q_U_display[inds]), label = "q")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("f vs q")

# using BenchmarkTools
#
# m = 1
# A = As[1];
#
# println("qs[i][k], gs eval:")
# r0 = 2*π*U[m] - A.ss_params.d[1]
# @btime A.qs[1][1](r0, 1.0)
# @btime A.gs[1][1](r0, 1.0)
#
# println("q eval.")
# @btime q.(U_rad[m]);
#
# println("q_U eval")
# @btime q.(U_rad);
# println()


text:
"""A helper function that computes the centroid of _data_"""
function centroid{T}(data::AbstractArray{T, 1})
  centroid = zero(T)
  sum = zero(T)
  dataLen = size(data, 1)
  i = 0

  for harmonic in data
    centroid += i*harmonic
    sum += harmonic
    i += 1
  end

  if abs(centroid) < one(T)
    return T(dataLen/2)
  end
  return centroid / sum
end

"""
Compute the centroid frequency of a signal as a function of time. This is done by taking the STFT and finding the centroid of the FFT of each time slice
Parameters:
    -- `data`: The signal to find the STFT of
    -- `windowFunc`: see STFT.stft_processed
    -- `windowStep`: The step ( in samples) to move the window between FFTs
    -- `low`: The lowest harmonic to include in the output
    -- `high`: The highest harmonic to include in the output
    -- `transformFunction`: A function to apply to the windowed, FFT'd pieces of _data_
"""
function centroidfrequencies{dataType,windowType}(data::AbstractArray{dataType, 1}, windowFunc::AbstractArray{windowType, 1}, windowStep::Integer, low::Integer, high::Integer, transformFunction::Function)
  nRows = length(data)
  curCol = view(data, 1:nRows, 1)
  centroids = stft_processed(data, windowFunc, windowStep, low, high, transformFunction, (dat, i)->Int(round(centroid(dat))), Int)

  return centroids
end

function centroidfrequencies{dataType}(data::AbstractArray{dataType, 2})
    centroidLength = size(data, 2)
    centroids = Array{dataType, 1}(centroidLength)

    for timeSlice in 1:centroidLength
	centroids[timeSlice] = centroid(view(data, :, timeSlice))
    end

    return centroids
end

export centroidfrequencies