text:
#!/usr/bin/env julia
using HDF5, DataStructures
using Pope: LJH
using DocOpt

doc = """
Pope Checker
Used to find folders where pope analysis didn't work quite right and re-run pope on them. Provide a <dir> to look for ljh containing directories in, and if you want the option to reanalyze the files, provide <pkfilename>. You will have a chance to review the candidates, and review the proposed commands, before anything actually happens. At most this will make a file checker.sh, which you can then run via bash checker.sh.
Usage:
	checker.jl [--missing_hdf5_is_candidate] <dir> [<pkfilename>]

Options:
	--missing_hdf5_is_candidate       Folders that lack a pope hdf5 file algtogther are not considered reanalysis candidates unless you provide this flag.
"""

immutable FolderStatus
	n_ljh_file::Int
	pope_hdf5_found::Bool
	n_channels_in_hdf5::Int
	n_channels_in_hdf5_matching_ljh_found::Int
	n_channels_hdf5_and_ljh_match_length::Int
	n_channels_filt_value_nonzero::Int
end
function Base.show(io::IO,fs::FolderStatus)
	if fs.pope_hdf5_found
		print(io, "$(fs.n_ljh_file) ljh files found, $(fs.n_channels_in_hdf5) channels in hdf5, $(fs.n_channels_in_hdf5_matching_ljh_found) channels match between ljh and hdf5, $(fs.n_channels_hdf5_and_ljh_match_length) channels have matching lengths. $(fs.n_channels_filt_value_nonzero) have nonzero first filt_value.")
	else
		print(io, "$(fs.n_ljh_file) ljh files found, no pope hdf5 found.")
	end
end

function get_pope_hdf5_fname(dir)
	a=filter(x->contains(x,"pope"),readdir(dir))
	if length(a) == 1
		Nullable(joinpath(dir,a[1]))
	else
		Nullable{String}()
	end
end

function comparechannels(popefilename, ljhdict,verbose)
	nchannels = 0
	nchannelsmatchlength = 0
	nchannelsfiltvaluenonzero = 0
	p = h5open(popefilename,"r")
	n_hdf5_channel = length(filter(s->startswith(s,"chan"),names(p)))
	for ch in 1:2:480
		chanstr = "chan$ch"
		ch in keys(ljhdict) || continue
		exists(p,chanstr) || continue
		nchannels+=1
		chan = p[chanstr]
		ljh = Pope.LJH.LJHFile(ljhdict[ch])
		if length(chan["filt_value"]) == length(ljh)
			nchannelsmatchlength+=1
		elseif verbose
			println("$chanstr pope_hdf5 has $(length(chan["filt_value"])) length, $(ljh.filename) has $(length(ljh))")
		end
		chan["filt_value"][1] !=0 && (nchannelsfiltvaluenonzero+=1)
	close(ljh)
	end
	close(p)
	return FolderStatus(length(ljhdict), true, n_hdf5_channel, nchannels, nchannelsmatchlength, nchannelsfiltvaluenonzero)
end

"comparechannels(d;verbose=false)
Look in directory `d` for a pope hdf5 file, and return a FolderStatus object containing info used to determine if the folder is a re-analysis candidate."
function comparechannels(d;verbose=false)
	pope_hdf5_nullable = get_pope_hdf5_fname(d)
	ljhdict = LJH.allchannels(d)
	if isnull(pope_hdf5_nullable)
		return FolderStatus(length(ljhdict), false, 0, 0, 0, 0)
	else
 		return comparechannels(get(pope_hdf5_nullable), ljhdict,verbose)
	end
end
function is_renalysis_candidate(f::FolderStatus, missing_hdf5_is_candidate::Bool)
	f.n_ljh_file==0 && return false
	!f.pope_hdf5_found && return missing_hdf5_is_candidate
	return f.n_channels_hdf5_and_ljh_match_length != f.n_channels_in_hdf5_matching_ljh_found
end

#test for is_reanalysis_candidate
# missing_hdf5_is_candidate
@assert !is_renalysis_candidate(FolderStatus(1,false,1,1,1,1), false)
@assert is_renalysis_candidate(FolderStatus(1,false,1,1,1,1), true)
# must have at least 1 ljh file
@assert !is_renalysis_candidate(FolderStatus(0,false,1,1,1,1), true)
@assert !is_renalysis_candidate(FolderStatus(0,true,1,1,1,1), true)
# length match (4th and 5th arguments match)
@assert !is_renalysis_candidate(FolderStatus(1,true,1,1,1,1), true)
@assert is_renalysis_candidate(FolderStatus(1,true,1,1,0,1), true)

function get_reanalysis_candidates(dirs, missing_hdf5_is_candidate::Bool, verbose::Bool)
	fsdict = OrderedDict(collect(d=>comparechannels(d,verbose=verbose) for d in dirs))
	filter((k,v)-> is_renalysis_candidate(v,missing_hdf5_is_candidate),fsdict)
end

function report_candidates(fsdict)
	println("Review these $(length(fsdict)) re-analysis candidates (must have ljh files in them, also depends on --missing_hdf5_is_candidate flag):")
	for (k,v) in fsdict
		println("$k")
		println("\tWhy? $v")
	end
end

function prompt_yes_no(s;exit_on_no=false)
	while true
		println(s)
		print("[y/n]?")
		r = lowercase(chomp(readline()))
		if startswith("yes",r) && length(r)>0
			return true
		elseif startswith("no",r) && length(r)>0
			if exit_on_no
				exit()
			else
				return false
			end
		end
	end
end

function make_reanalysis_command(pkfilename, candidate)
	k,v=candidate
	chan1 = LJH.fnames(k,1)
	outputpath = LJH.hdf5_name_from_ljh(chan1)
	"./popeonce.jl --overwriteoutput $k $pkfilename $outputpath"
end

arguments = docopt(doc, version=v"0.0.1")

missing_hdf5_is_candidate = arguments["--missing_hdf5_is_candidate"]
pkfilename = arguments["<pkfilename>"] # nothing if missing, String if input
startdir = arguments["<dir>"]
println("Searching $startdir")
lscmd = `ls $startdir`
println(lscmd)
println(chomp(readstring(lscmd)))
searchdirslocal = [d for d in readdir(startdir) if isdir(joinpath(startdir,d))]
searchdirs = [joinpath(startdir,d) for d in searchdirslocal]
println("Found $(length(searchdirs)) subdirs")
candidates = get_reanalysis_candidates(searchdirs,missing_hdf5_is_candidate,false)
#for (k,v) in candidates
#	println(k)
#	iscandidate = is_renalysis_candidate(v,missing_hdf5_is_candidate)
#	println("Is candidate? $iscandidate")
#	println(v)
#end
report_candidates(candidates)
println("Preknowledge filename: $pkfilename")
commands = collect(make_reanalysis_command("$pkfilename",c) for c in candidates)
prompt_yes_no("Should this program proceed to renalyze all shown candidates?",exit_on_no=true)
prompt_yes_no("Use preknowledge file: $pkfilename?",exit_on_no=true)
println("Proposed Commands")
for c in commands println(c) end
if pkfilename == nothing
	println("You must pass pkfilename to re-analyze, commands show have nothing as a stand in name.")
	exit()
end
prompt_yes_no("Do the proposed commands look ok?", exit_on_no=true)
println("Writing proposed commands to checker.sh. You can do `bash checker.sh to execute them.")
open("checker.sh","w") do f
for c in commands write(f,string(c)*"\n") end
end

ext: jl


text:
using LogDensityProblems, Test, Distributions, TransformVariables
import LogDensityProblems: capabilities, dimension, logdensity
using LogDensityProblems: logdensity_and_gradient, LogDensityOrder

import ForwardDiff, Tracker, TransformVariables, Random, Zygote
using UnPack: @unpack

####
#### test setup and utilities
####

###
### reproducible randomness
###

Random.seed!(1)

###
### comparisons (for testing)
###

"""
    a ≅ b

Compare log denfields and types, for unit testing.
"""
≅(::Any, ::Any, atol = 0) = false

function ≅(a::Real, b::Real, atol = 0)
    if isnan(a)
        isnan(b)
    elseif isinf(a)
        a == b
    else
        abs(a - b) ≤ atol
    end
end

function ≅(a::Tuple{Real,Any}, b::Tuple{Real,Any}, atol = 0)
    ≅(first(a), first(b), atol) || return false
    !isfinite(first(a)) || isapprox(last(a), last(b); atol = atol, rtol = 0)
end

@testset "comparisons for unit testing" begin
    @test 1 ≅ 1
    @test !(1 ≅ 2)
    @test Inf ≅ Inf
    @test (1, [1, 2]) ≅ (1, [1, 2])
    @test !((1, [1, 2]) ≅ (1, [1, 3]))
    @test !((3, [1, 2]) ≅ (1, [1, 2]))
    @test (-Inf, [1, 2]) ≅ (-Inf, [1, 2])
    @test (-Inf, [1, 2]) ≅ (-Inf, [1, 3])
    @test (-Inf, [1, 2]) ≅ (-Inf, nothing)
end

###
### simple log density for testing
###

struct TestLogDensity{F}
    ℓ::F
end
logdensity(ℓ::TestLogDensity, x) = ℓ.ℓ(x)
dimension(::TestLogDensity) = 3
test_logdensity1(x) = -2*abs2(x[1]) - 3*abs2(x[2]) - 5*abs2(x[3])
test_logdensity(x) = any(x .< 0) ? -Inf : test_logdensity1(x)
test_gradient(x) = x .* [-4, -6, -10]
TestLogDensity() = TestLogDensity(test_logdensity) # default: -Inf for negative input
Base.show(io::IO, ::TestLogDensity) = print(io, "TestLogDensity")

####
#### traits
####

@test capabilities("a fish") ≡ nothing

@testset "LogDensityOrder" begin
    @test LogDensityOrder(1) == LogDensityOrder(1)
    @test_throws ArgumentError LogDensityOrder(-1)
    @test LogDensityOrder(2) ≥ LogDensityOrder(1)
end

####
#### AD backends
####

@testset "AD via ForwardDiff" begin
    ℓ = TestLogDensity()
    ∇ℓ = ADgradient(:ForwardDiff, ℓ)
    @test repr(∇ℓ) == "ForwardDiff AD wrapper for " * repr(ℓ) * ", w/ chunk size 3"
    @test dimension(∇ℓ) == 3
    @test capabilities(∇ℓ) ≡ LogDensityOrder(1)
    for _ in 1:100
        x = randn(3)
        @test @inferred(logdensity(∇ℓ, x)) ≅ test_logdensity(x)
        @test @inferred(logdensity_and_gradient(∇ℓ, x)) ≅
            (test_logdensity(x), test_gradient(x))
    end
end

@testset "chunk heuristics for ForwardDiff" begin
    @test LogDensityProblems.heuristic_chunks(82) == vcat(1:4:81, [82])
end

@testset "benchmark ForwardDiff chunk size" begin
    ℓ = TransformedLogDensity(as(Array, 20), x -> -sum(abs2, x))
    b = LogDensityProblems.benchmark_ForwardDiff_chunks(ℓ)
    @test b isa Vector{Pair{Int,Float64}}
    @test length(b) ≤ 20
end

@testset "AD via Tracker" begin
    ℓ = TestLogDensity()
    ∇ℓ = ADgradient(:Tracker, ℓ)
    @test repr(∇ℓ) == "Tracker AD wrapper for " * repr(ℓ)
    @test dimension(∇ℓ) == 3
    @test capabilities(∇ℓ) ≡ LogDensityOrder(1)
    for _ in 1:100
        x = randn(3)
        @test @inferred(logdensity(∇ℓ, x)) ≅ test_logdensity(x)
        @test @inferred(logdensity_and_gradient(∇ℓ, x)) ≅ (test_logdensity(x), test_gradient(x))
   end
end

@testset "AD via Zygote" begin
    ℓ = TestLogDensity(test_logdensity1)
    ∇ℓ = ADgradient(:Zygote, ℓ)
    @test repr(∇ℓ) == "Zygote AD wrapper for " * repr(ℓ)
    @test dimension(∇ℓ) == 3
    @test capabilities(∇ℓ) ≡ LogDensityOrder(1)
    for _ in 1:100
        x = randn(3)
        @test @inferred(logdensity(∇ℓ, x)) ≅ test_logdensity1(x)
        @test logdensity_and_gradient(∇ℓ, x) ≅ (test_logdensity1(x), test_gradient(x))
    end
end

@testset "ADgradient missing method" begin
    msg = "Don't know how to AD with Foo, consider `import Foo` if there is such a package."
    P = TransformedLogDensity(as(Array, 1), x -> sum(abs2, x))
    @test_logs((:info, msg), @test_throws(MethodError, ADgradient(:Foo, P)))
end

####
#### transformed Bayesian problem
####

@testset "transformed Bayesian problem" begin
    t = as((y = asℝ₊, ))
    d = LogNormal(1.0, 2.0)
    logposterior = ((x, ), ) -> logpdf(d, x)

    # a Bayesian problem
    p = TransformedLogDensity(t, logposterior)
    @test repr(p) == "TransformedLogDensity of dimension 1"
    @test dimension(p) == 1
    @test p.transformation ≡ t
    @test capabilities(p) == LogDensityOrder(0)

    # gradient of a problem
    ∇p = ADgradient(:ForwardDiff, p)
    @test dimension(∇p) == 1
    @test parent(∇p).transformation ≡ t

    for _ in 1:100
        x = random_arg(t)
        θ, lj = transform_and_logjac(t, x)
        px = logdensity(p, x)
        @test logpdf(d, θ.y) + lj ≈ (px::Real)
        px2, ∇px = logdensity_and_gradient(∇p, x)
        @test px2 == px
        @test ∇px ≈ [ForwardDiff.derivative(x -> logpdf(d, exp(x)) + x, x[1])]
    end
end

@testset "-∞ log densities" begin
    t = as(Array, 2)
    validx = x -> all(x .> 0)
    p = TransformedLogDensity(t, x -> validx(x) ?  sum(abs2, x)/2 : -Inf)
    ∇p = ADgradient(:ForwardDiff, p)

    @test dimension(p) == dimension(∇p) == TransformVariables.dimension(t)
    @test p.transformation ≡ parent(∇p).transformation ≡ t

    for _ in 1:100
        x = random_arg(t)
        px = logdensity(∇p, x)
        px_∇px = logdensity_and_gradient(∇p, x)
        @test px isa Real
        @test px_∇px isa Tuple{Real,Any}
        @test first(px_∇px) ≡ px
        if validx(x)
            @test px ≈ sum(abs2, x)/2
            @test last(px_∇px) ≈ x
        else
            @test isinf(px)
        end
    end
end

@testset "stresstest" begin
    @info "stress testing"
    ℓ = TestLogDensity(x -> all(x .< 0) ? error("invalid") : -sum(abs2, x))
    failures = LogDensityProblems.stresstest(logdensity, ℓ; N = 500)
    @test 50 ≤ length(failures) ≤ 100
    @test all(x -> all(x .< 0), failures)
end


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

res_f_1p = params_intervals(
    [3.],
    1,
    9.,
    f_1p,
    logscale_all = false,
    method = :D2D_PLE
)

res_f_2p_1im = [params_intervals(
    [3., 1.],
    i,
    9.,
    f_2p_1im,
    logscale_all = false,
    method = :D2D_PLE
) for i in 1:2]

res_f_2p = [params_intervals(
    [3., 4.],
    i,
    9.,
    f_2p,
    logscale_all = false,
    method = :D2D_PLE
) for i in 1:2]

res_f_3p_1im = [params_intervals(
    [3., 4., 1.1],
    i,
    9.,
    f_3p_1im,
    #max_iter = 1000000,
    logscale_all = false,
    method = :D2D_PLE
) for i in 1:3]

res_f_4p_2im = [params_intervals(
    [3., 4, 1.1, 7.],
    i,
    9.,
    f_4p_2im,
    logscale_all = false,
    method = :D2D_PLE
) for i in 1:4]

res_f_4p_3im = [params_intervals(
    [3., 4, 1.1, 7.],
    i,
    9.,
    f_4p_3im,
    logscale_all = false,
    method = :D2D_PLE
) for i in 1:4]

res_f_5p_3im = [params_intervals(
    [3., 0.1, 4, 1.1, 8.],
    i,
    9.,
    f_5p_3im,
    logscale_all = false,
    method = :D2D_PLE
) for i in 1:5]

res_f_3p_im = [params_intervals(
    [3., 0.1, 8.],
    i,
    9.,
    f_3p_im,
    scan_bound = [0.,10^3],
    logscale_all = false,
    method = :D2D_PLE
) for i in 1:3]
# tests
@testset "f_1p" begin
    @test all(@. isapprox(res_f_1p.interval, [1.0, 5.0], atol=res_f_1p.input.ptol))
    @test all(@. isapprox(res_f_1p.loss_final, [res_f_1p.input.loss_crit, res_f_1p.input.loss_crit], atol=res_f_1p.input.losstol))
    @test res_f_1p.ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
end

@testset "f_2p_1im" begin
    @test all(@. isapprox(res_f_2p_1im[1].interval, [1.0, 5.0], atol=res_f_2p_1im[1].input.ptol))
    @test all(@. isapprox(res_f_2p_1im[1].loss_final, [res_f_2p_1im[1].input.loss_crit, res_f_2p_1im[1].input.loss_crit],
    atol=res_f_2p_1im[1].input.losstol))
    @test res_f_2p_1im[1].ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
    @test res_f_2p_1im[2].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
end

@testset "f_2p" begin
    @test all(@. isapprox(res_f_2p[1].interval, [1.0, 5.0], atol=res_f_2p[1].input.ptol))
    @test all(@. isapprox(res_f_2p[1].loss_final, [res_f_2p[1].input.loss_crit, res_f_2p[1].input.loss_crit],
    atol=res_f_2p[1].input.losstol))
    @test res_f_2p[1].ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
    @test all(@. isapprox(res_f_2p[2].interval, [2.0, 6.0], atol=res_f_2p[2].input.ptol))
    @test all(@. isapprox(res_f_2p[2].loss_final, [res_f_2p[2].input.loss_crit, res_f_2p[2].input.loss_crit],
    atol=res_f_2p[2].input.losstol))
    @test res_f_2p[2].ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
end

@testset "f_3p_1im" begin
    @test all(@. isapprox(res_f_3p_1im[1].interval, [1.0, 5.0], atol=res_f_3p_1im[1].input.ptol))
    @test all(@. isapprox(res_f_3p_1im[1].loss_final, [res_f_3p_1im[1].input.loss_crit, res_f_3p_1im[1].input.loss_crit],
    atol=res_f_3p_1im[1].input.losstol))
    @test res_f_3p_1im[1].ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
    @test res_f_3p_1im[2].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
    @test res_f_3p_1im[3].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
end

@testset "f_4p_2im" begin
    @test all(@. isapprox(res_f_4p_2im[1].interval, [1.0, 5.0], atol=res_f_4p_2im[1].input.ptol))
    @test all(@. isapprox(res_f_4p_2im[1].loss_final, [res_f_4p_2im[1].input.loss_crit, res_f_4p_2im[1].input.loss_crit],
    atol=res_f_4p_2im[1].input.losstol))
    @test res_f_4p_2im[1].ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
    @test all(@. isapprox(res_f_4p_2im[2].interval, [2.0, 6.0], atol=res_f_4p_2im[2].input.ptol))
    @test all(@. isapprox(res_f_4p_2im[2].loss_final, [res_f_4p_2im[2].input.loss_crit, res_f_4p_2im[2].input.loss_crit],
    atol=res_f_4p_2im[2].input.losstol))
    @test res_f_4p_2im[2].ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
    @test res_f_4p_2im[3].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
    @test res_f_4p_2im[4].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
end

@testset "f_4p_3im" begin
    @test all(@. isapprox(res_f_4p_3im[1].interval, [1.0, 5.0], atol=res_f_4p_3im[1].input.ptol))
    @test all(@. isapprox(res_f_4p_3im[1].loss_final, [res_f_4p_3im[1].input.loss_crit, res_f_4p_3im[1].input.loss_crit],
    atol=res_f_4p_3im[1].input.losstol))
    @test res_f_4p_3im[1].ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
    @test res_f_4p_3im[2].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
    @test res_f_4p_3im[3].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
    @test res_f_4p_3im[4].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
end

@testset "f_5p_3im" begin
    @test all(@. isapprox(log10.(res_f_5p_3im[1].interval), log10.([1.0, 5.0]), atol=res_f_5p_3im[1].input.ptol))
    @test all(@. isapprox(res_f_5p_3im[1].loss_final, [res_f_5p_3im[1].input.loss_crit, res_f_5p_3im[1].input.loss_crit],
    atol=res_f_5p_3im[1].input.losstol))
    @test res_f_5p_3im[1].ret_codes == [:FTOL_REACHED, :FTOL_REACHED]
    @test isapprox(log10(res_f_5p_3im[2].interval[2]), log10(1.1), atol=res_f_5p_3im[2].input.ptol)
    @test res_f_5p_3im[2].ret_codes == [:BOUNDS_REACHED, :FTOL_REACHED]
    @test res_f_5p_3im[3].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
    @test res_f_5p_3im[4].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
    @test res_f_5p_3im[5].ret_codes == [:BOUNDS_REACHED, :BOUNDS_REACHED]
end


text:
using StatsPlots

"""
    coinflip_discrete_bayes_update(p_θ, data)

Given a discrete (possibly unnormalized) prior and some binary data
returns the likelihood and posterior.
"""
function coinflip_discrete_bayes_update(p_θ, data; plot_result=true)
    n_vals = length(p_θ) # number of candidates for θ
    θ = collect(0 : 1/(n_vals -1) : 1) # equally spaced bins
    p_θ = p_θ / sum(p_θ)     # Makes sure that beliefs sum to 1.
    heads = sum(data)
    tails = length(data) - heads
    likelihood = θ.^heads .* (1 .- θ).^tails
    p_data = sum(likelihood .* p_θ)
    posterior = (likelihood .* p_θ) ./ p_data
    return (likelihood = likelihood, posterior = posterior)
end

"""
    plot_coinflip_discrete_bayes_update(p_θ, data)

Reproduces figure 5.1 with custom prior and data.
To reproduce the figure exactly use:
```julia-repl
julia> θ = range(0.1,1, step=0.1)
julia> p_θ = min.(θ, 1 .- θ)
julia> plot_coinflip_discrete_bayes_update(p_θ, [1])
```
"""
function plot_coinflip_discrete_bayes_update(p_θ, data)
    n_vals = length(p_θ)
    θ = collect(0 : 1/(n_vals -1) : 1)
    p_θ = p_θ / sum(p_θ)
    likelihood, posterior = coinflip_discrete_bayes_update(p_θ, data)
    l = @layout (3,1)
    # Plot the prior
    plot1 = bar(θ, p_θ, title="Prior", leg=false, bar_width = 0.01)

    # Plot likelihood
    plot2 = bar(θ, likelihood, title="Likelihood", leg=false, bar_width = 0.01)

    # Plot the posterior
    plot3 = bar(θ, posterior, title="Posterior", leg=false, bar_width = 0.01)
    plt = plot(plot1, plot2, plot3, layout=l, xlims=(-0.02,1.02), xticks=round.(θ,digits=2))
    display(plt)
end


text:
using Base.Iterators: partition
# Flux stuff
using Flux, CUDA
import Flux.Losses: mse
using Flux.Optimise: Momentum, update!
# Image stuff
using Images, ImageIO, ImageMagick
# Plots
using Plots
# Datasets
using MLDatasets
using Statistics
using Parameters
# Other stuff
using LinearAlgebra
using Random
using ProgressBars
import Distributions: Uniform

using Suppressor

@with_kw mutable struct HPARAMS
	batch_size::Int = 64
	epochs::Int = 200
	# feedback frequency
	verbose_freq::Int = 500
	# learning rate
	lr::Float64 = 0.001
	# tolerance on ratios
	tolerance::Float64 = 0.001
	threshold::Float64 = 0.3
	max_loop::Int = 20
end

# load full training set
#train_x, train_y = MNIST.traindata()
train_x, train_y = FashionMNIST.traindata()
# load full test set
#test_x,  test_y  = MNIST.testdata()
test_x,  test_y  = FashionMNIST.testdata()
function model1(hparams)
	return Chain(
		x -> reshape(x,28*28,:),
		Dense(28*28,512),
		x -> leakyrelu.(x,0.3f0),
		Dropout(0.2),
		Dense(512, 32),
		x -> leakyrelu.(x, 0.3f0),
		Dropout(0.2),
		Dense(32,1),
		x -> leakyrelu.(x, 0.2f0)
		)
end

function model2(hparams)
	return Chain(
		x -> reshape(x,28*28,:),
		Dense(28*28,512),
		x -> leakyrelu.(x,0.3f0),
		Dropout(0.2),
		Dense(512, 32),
		x -> leakyrelu.(x, 0.3f0),
		Dropout(0.2),
		Dense(32,1),
		x -> leakyrelu.(x, 0.2f0)

		)
end

#=
function model(hparams)
	return Chain(
		Conv((2,2), 1 => 4; stride=1, pad=2),
		x -> leakyrelu.(x, 0.2f0),
		Dropout(0.2),
		Conv((2,2), 4=>16; stride=1, pad=2),
		x -> leakyrelu.(x, 0.2f0),
		Dropout(0.2),
		Conv((3,3), 16=>8; stride=1, pad=1),
		x -> leakyrelu.(x, 0.2f0),
		Dropout(0.2),
		Conv((3,3), 8=>4; stride=1, pad=1),
		x -> leakyrelu.(x,0.2f0),
		Dropout(0.2),
		Conv((2,2), 4=>1; stride=1, pad=1),
		x -> leakyrelu.(x, 0.2f0),
		Dropout(0.2),
		x -> reshape(x,100,:),
		Dense(100,hparams.batch_size),
		x -> leakyrelu.(x,0.2f0),
		x -> sum(x,dims=2),
		x -> reshape(x,hparams.batch_size),
		Dense(hparams.batch_size,hparams.batch_size),
		x -> sigmoid.(x)
	)
end
=#

function discrete_observation(from_nn,target,hparams)
    max_target = maximum(target)
    output = round.(from_nn .* max_target) .- max_target*0.5

    ratio_positive = length([i for i in output if i>0]) / length(output)
    ratio_negative = length([i for i in output if i<0]) / length(output)
	target_ratio_p = length([i for i in target if i>0]) / length(output)
	target_ratio_n = length([i for i in target if i<0]) / length(target)

    ratio = ratio_positive - ratio_negative
	target_ratio = target_ratio_p - target_ratio_n

	if ratio < hparams.tolerance
		return 1.0
	else
		return (target_ratio-ratio)/(target_ratio + ratio)
	end
end

function criterion(output,target)
    return Flux.mse(output,target)
end

function train_model(x,y,model, optim, switch, hparams)
    ps = Flux.params(model)
    a = model(x)
    y = reshape(y,size(a))
    loss = 0.0

    if switch
        factor = discrete_observation(a,y,hparams)
		counter = 0
		while factor > hparams.threshold
			gs = gradient(ps) do
				loss = criterion(model(x),y)
			end
			update!(optim,ps,gs)
			a = model(x)
			factor = discrete_observation(a,y,hparams)
			counter += 1
			if counter > hparams.max_loop
				return loss
			end
		end
    else
        factor = 1.0
    end

    gs = gradient(ps) do
        loss = criterion(model(x),y)
    end
    update!(optim, ps, gs)

    return loss
end

function less_than(a,b)
	if a < b
		return 1.0
	else
		return 0.0
	end
end

accuracy(x,y,m) = mean(onecold(cpu(m(x)), 1:10) .== onecold(cpu(y), 1:10))

function test(model_interest, x, y, hparams)
    x = reshape(@.(2f0 .* x .- 1f0), 28, 28, :)
    y = reshape(y ./ 9.0, 10000)
    data_in = [x[:, :, r] for r in partition(1:10000, hparams.batch_size)];
    targets = [y[r] for r in partition(1:10000, hparams.batch_size)];
    tolerance = 0.5
	counter = 0
    max_acc = 0.0
    for x in data_in
        counter += 1
		val = targets[counter]
		pred = abs.(model_interest(x)) .* 9.0
		val = abs.(reshape(val,size(pred))) .* 9.0
		acc = Statistics.mean([less_than(i,tolerance) for i in abs.(pred .- val)])
		if max_acc < acc
			max_acc = acc
		end
		if counter > 154
			return 100.0*max_acc
		end
	end

	return 100.0*max_acc
end

function train(; kws...)
    hparams = HPARAMS(; kws...)
    # normalize the data to -1 to 1
	image_tensor = reshape(@.(2f0 .* train_x .- 1f0), 28, 28,:);
	answer_tensor = reshape(train_y/9.0, 60000)
	# parition the above tensor to batches
	data_x = [image_tensor[:, :, r] for r in partition(1:60000, hparams.batch_size)];
	data_y = [answer_tensor[r] for r in partition(1:60000, hparams.batch_size)];
	# models
    without_observation = model1(hparams) |> gpu
    with_observation = model2(hparams) |> gpu
    # optimisers
    no_o_optim = Momentum(hparams.lr,0.9)
    o_optim = Momentum(hparams.lr,0.9)
	train_step = 0

	if ispath("./logs")
		rm("./logs",recursive=true)
		mkpath("./logs")
	else
		mkpath("./logs")
	end

	a = "./logs/logs.csv"
	first_line = "Step, Accuracy with an observer, Accuracy without an observer\n"
	io = open(a, "a")
   for epoch in 1:hparams.epochs
        counter = 0
        for x in ProgressBar(data_x)
            counter = counter + 1 # this will count the entries of data_y
            x = cu(x)
            y = cu(data_y[counter])
            #if counter == length(data_y)-1
            #    break
            #end
            # now pass the data through the networks
        without_loss = train_model(x,y,without_observation,no_o_optim,false,hparams)
        with_loss = train_model(x,y,with_observation, o_optim,true,hparams)

			# iterate the train_step
			train_step += 1

            # test
            if train_step % hparams.verbose_freq == 0
                without = test(without_observation, cu(test_x),cu(test_y), hparams)
                with = test(with_observation, cu(test_x),cu(test_y), hparams)
				println()
				#println("At $(epoch) epoch , $(counter) step.")
                #println("Without observations, the loss is : ", without_loss)
                #println("With observations, the loss is : ", with_loss)
                #println("Without observations, we have an accuracy of $(without)%")
                #println("With observations, we have an accuracy of $(with)%")
				printer = "$(train_step), $(with)%, $(without)%\n"
				write(io,printer)
            end # end of test
        end # end of inner for-loop
    end # end of epoch for-loop
	close(io)
end

@suppress_err begin
    train()
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b11fdebe-5507-11eb-1001-537a20e51d0d
begin
	let
		env = mktempdir()
		import Pkg
		Pkg.activate(env)
		Pkg.Registry.update()
		Pkg.add("FileIO")
		Pkg.add("JLD2")
		Pkg.add("StatsBase")
		Pkg.add("StatsPlots")
		Pkg.add("SimpleANOVA")
		Pkg.add(path="../..")
	end
	using SMLMAssociationAnalysis_NCB
	using FileIO
	using StatsBase
	using Statistics
	using StatsPlots
	using SimpleANOVA
end

# ╔═╡ 56dcb910-5507-11eb-0b0a-a1c2aef89207
md"""
# Analysis of Experiment 1: p53-Mdm2
"""

# ╔═╡ 90d2acb0-5507-11eb-39c4-df93c258dd98
md"""Load packages"""

# ╔═╡ 86502ff0-550d-11eb-1c13-553251c061ee
md"""Load results file and compute normalized values"""

# ╔═╡ 5e3db410-550d-11eb-0331-49a447464b07
begin
	### Load saved data

	samplenames = ["A", "B", "C", "D"]

	nreplicates = 3
	nsamples = 4
	ncells = 10

	outputdir = "../../output"
	datapath = joinpath(outputdir, "results.jld2")

	experimentresults = load(datapath)["experimentresults"]

	medianmeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
	montecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
	positivecontrolmontecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)
	negativecontrolmontecarlomeasurements = Array{Float64,4}(undef, ncells, nreplicates, 4, 2)

	for k ∈ 1:2
		for i ∈ 1:nsamples
			samplemedianresults = Array{Float64,2}(undef, ncells, nreplicates)
			samplelessthan10results = Array{Float64,2}(undef, ncells, nreplicates)
			samplepositivecontrollessthan10results = Array{Float64,2}(undef, ncells, nreplicates)
			samplenegativecontrollessthan10results = Array{Float64,2}(undef, ncells, nreplicates)

			for j ∈ 1:nreplicates
				replicateresults = experimentresults[k][j][i]
				lessthanlimitreplicate = [(x.distances .< 200) .& (x.percentileranks .< 0.1) for x ∈ replicateresults]
				positivecontrollessthanlimitreplicate = [(x.positivecontrol_distances .< 200) .& (x.positivecontrol_percentileranks .< 0.1) for x ∈ replicateresults]
				negativecontrollessthanlimitreplicate = median.([count.([(x.negativecontrol_distances[m] .< 200) .& (x.negativecontrol_percentileranks[m] .< 0.1) for m ∈ eachindex(x.negativecontrol_distances)]) for x ∈ replicateresults])
				lessthan10 = count.(lessthanlimitreplicate) ./ length.(lessthanlimitreplicate)
				positivecontrollessthan10 = count.(positivecontrollessthanlimitreplicate) ./ length.(lessthanlimitreplicate)
				negativecontrollessthan10 = negativecontrollessthanlimitreplicate ./ length.(lessthanlimitreplicate)
				mediandistances = map(x -> x.mediandistance, replicateresults)
				samplemedianresults[:, j] = mediandistances
				samplelessthan10results[:, j] = lessthan10
				samplepositivecontrollessthan10results[:, j] = positivecontrollessthan10
				samplenegativecontrollessthan10results[:, j] = negativecontrollessthan10
			end
			medianmeasurements[:, :, i, k] = samplemedianresults
			montecarlomeasurements[:, :, i, k] = samplelessthan10results
			positivecontrolmontecarlomeasurements[:, :, i, k] = samplepositivecontrollessthan10results
			negativecontrolmontecarlomeasurements[:, :, i, k] = samplenegativecontrollessthan10results
		end
	end

	medianmeasurements = cat(
		cat(medianmeasurements[:, :, 1:2, 1], medianmeasurements[:, :, 3:4, 1], dims = 4),
		cat(medianmeasurements[:, :, 1:2, 2], medianmeasurements[:, :, 3:4, 2], dims = 4),
		dims = 5,
	)
	montecarlomeasurements = cat(
		cat(montecarlomeasurements[:, :, 1:2, 1], montecarlomeasurements[:, :, 3:4, 1], dims = 4),
		cat(montecarlomeasurements[:, :, 1:2, 2], montecarlomeasurements[:, :, 3:4, 2], dims = 4),
		dims = 5,
	)
	positivecontrolmontecarlomeasurements = cat(
		cat(positivecontrolmontecarlomeasurements[:, :, 1:2, 1], positivecontrolmontecarlomeasurements[:, :, 3:4, 1], dims = 4),
		cat(positivecontrolmontecarlomeasurements[:, :, 1:2, 2], positivecontrolmontecarlomeasurements[:, :, 3:4, 2], dims = 4),
		dims = 5,
	)
	negativecontrolmontecarlomeasurements = cat(
		cat(negativecontrolmontecarlomeasurements[:, :, 1:2, 1], negativecontrolmontecarlomeasurements[:, :, 3:4, 1], dims = 4),
		cat(negativecontrolmontecarlomeasurements[:, :, 1:2, 2], negativecontrolmontecarlomeasurements[:, :, 3:4, 2], dims = 4),
		dims = 5,
	)

	normalizedmontecarlomeasurements = (montecarlomeasurements .- negativecontrolmontecarlomeasurements) ./ (positivecontrolmontecarlomeasurements .- negativecontrolmontecarlomeasurements)

	nothing
end

# ╔═╡ 294e6a02-550e-11eb-3450-d9d5989ec08c
md"""
Relevant data are stored in these variables:
- `medianmeasurements` - median distance between the pairs
- `montecarlomeasurements` - fraction of pairs considered associated (p < 0.1, dist < 200 nm)
- `normalizedmontecarlomeasurements` - fraction of pairs considered associated, normalized bassed on simulations of 0% and 100% binding

The data is stored in matrices with these dimensions:
- 1 - cell
- 2 - experimental replicate
- 3 - condition 1 (1 = -Dox, 2 = +Dox)
- 4 - condition 2 (Exp 1: 1 = -Nut, 2 = +Nut; Exp 2: 1 = MEG3, 2 = GAPDH)
- 5 - experiment (1 = p53-Mdm2, 2 = p53-MEG3)
"""

# ╔═╡ 9a94f370-5511-11eb-1f60-7fe83ed11def
md"""
### Planning
Central Limit Theorem - Only 10 measurements per level, but 30 per condition. Might be close enough for a normal distribution.

Group sizes equal

Doxycycline and Nutlin treatments should be independent.
"""

# ╔═╡ 572ba5f2-5513-11eb-1937-3f03cec8252e
md"""
### Medians
#### Check unusual cases
"""

# ╔═╡ 864b181e-5513-11eb-0304-8983d37a8170
let
	p1 = boxplot(medianmeasurements[:, :, 1, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
	p2 = boxplot(medianmeasurements[:, :, 2, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
	p3 = boxplot(medianmeasurements[:, :, 1, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))
	p4 = boxplot(medianmeasurements[:, :, 2, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Median distance (nm)"))

	plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "MDM2-p53")
end

# ╔═╡ 27e24820-5514-11eb-1415-8711a194892d
md"""
Potential outliers in box plots:
"""

# ╔═╡ f0393540-5514-11eb-2f76-3514cc3bc85c
[medianmeasurements[[3,4],2,1,1,1],
 medianmeasurements[8,1,2,1,1],
 medianmeasurements[6,2,2,1,1],
 medianmeasurements[10,3,1,2,1],
 medianmeasurements[5,3,2,2,1]]

# ╔═╡ f02df530-5515-11eb-0d9e-8dbb63ed2023
md"""Evaluating z-scores within each replicate:"""

# ╔═╡ 1bf98170-5516-11eb-0c24-77ac43672fba
mediansflat = [medianmeasurements[:, :, 1, 1, 1] medianmeasurements[:, :, 2, 1, 1] medianmeasurements[:, :, 1, 2, 1] medianmeasurements[:, :, 2, 2, 1]]

# ╔═╡ 25253910-5516-11eb-0a7a-1d73ed51d4a1
medianzscores = [zscore(mediansflat[:,i]) for i ∈ axes(mediansflat, 2)]

# ╔═╡ d8c96c5e-5517-11eb-0861-056d5e425f04
md"""
- In [1,1] Replicate 1: All within 1.96
- In [1,1] Replicate 2: All within 1.96
- In [1,1] Replicate 3: All within 1.96
- In [1,2] Replicate 1: 1 value, Cell 8 (10%) outside 1.96, but just barely. None outside 2.58.
- In [1,2] Replicate 2: 1 value, Cell 6 (10%) outside 1.96, by a lot. None outside 2.58.
- In [1,2] Replicate 3: All within 1.96
- In [2,1] Replicate 1: All within 1.96
- In [2,1] Replicate 2: All within 1.96
- In [2,1] Replicate 3: 1 value, Cell 10 (10%) outside 1.96 significantly. None outside 2.58.
- In [2,2] Replicate 1: All within 1.96
- In [2,2] Replicate 2: All within 1.96
- In [2,2] Replicate 3: 1 value, Cell 5 (10%) outside 1.96 significantly. None outside 2.58.
"""

# ╔═╡ 87eb8f70-5518-11eb-1639-f91b4bcb5c9b
md"""
None of these potential outliers are due to data entry error or measurement error. Sampling a different population is possible but unknown. These potential outliers are relatively modest. In addition, since an ANOVA will be used, must be careful not to delete points from a balanced design. Therefore, no outliers will be removed. However, I will Windsorize the more extreme ones for [1,2] Replicate 2, [2,1] Replicate 3, and [2,2] Replicate 3.
"""

# ╔═╡ 9b7f65a2-5952-11eb-0c21-733df3703bf9
begin
	medianmeasurementsw = copy(medianmeasurements)
	mediansflatw = copy(mediansflat)
	medianmeasurementsw[6,2,2,1,1] = medianmeasurements[7,2,2,1,1]
	medianmeasurementsw[10,3,1,2,1] = medianmeasurements[4,3,1,2,1]
	medianmeasurementsw[5,3,2,2,1] = medianmeasurements[9,3,2,2,1]
	mediansflatw[6,5] = mediansflat[7,5]
	mediansflatw[10,9] = mediansflat[4,9]
	mediansflatw[5,12] = mediansflat[9,12]
end;

# ╔═╡ a87aac80-5518-11eb-38d4-5fe90140b820
md"""
#### Check assumption of homoscedasticity
"""

# ╔═╡ f6e23cd0-5518-11eb-226a-c5751eb65b34
md"""
Based on the ZResid/ZPred plot, the data appears fairly heteroscedastic.

Conducted Levene's test, a 1-way ANOVA of the absolute deviation between each value and the group mean, by replicate.
"""

# ╔═╡ 0ffdd990-5519-11eb-376b-33be0cc966d7
levene(mediansflatw)

# ╔═╡ 28a8da30-5519-11eb-2fb0-efe693c6e44c
md"""
By Levene's test, the data is very likely heteroscedastic. Based on median, still significant at p = 0.011.

Not doing Hartley's Fmax because it isn't as simple.

Group sizes are equal, but samples are on the smaller side. CIs and significance tests may be affected.
"""

# ╔═╡ 5f8de452-5519-11eb-0a2c-051d3e701dd5
md"""
#### Check assumption of normality
"""

# ╔═╡ c4dbef60-55bd-11eb-1822-49c3473000e8
md"""
There are a couple notable deviations but for the *most* part they appear close to normal. The worst ones had the potential outliers identified previously.
"""

# ╔═╡ cf59a26e-55bd-11eb-1c4b-316a603b9007
medianskewness = [skewness(mediansflatw[:,i]) for i ∈ axes(mediansflatw, 2)]

# ╔═╡ 135dbe6e-55be-11eb-3a89-dd3d02c3ad4a
mediankurtosis = [kurtosis(mediansflatw[:,i]) for i ∈ axes(mediansflatw, 2)]

# ╔═╡ 6778a880-55be-11eb-1016-c59335beef0a
md"""
From skewness and kurtosis statistics, all were within 1.96 (5% level).
"""

# ╔═╡ 9bf9aaf0-55be-11eb-17a2-75b3eec550f8
md"""
The primary source of non-normality appear to be outliers, relatively minimal. As a result, I don't intend to do anything special.
"""

# ╔═╡ 265006c0-55c1-11eb-2ba8-3d628d096704
md"""
#### Conduct ANOVA tests
"""

# ╔═╡ 5f1ce3b0-55c1-11eb-3d03-b53516cfc976
medianresult = anova(
    permutedims(medianmeasurementsw[:, :, :, :, 1], (1,3,4,2)),
    [fixed, fixed, subject],
    factornames = ["Doxycycline", "Nutlin-3a", "Replicate"],
)

# ╔═╡ 756bf2ee-55c1-11eb-2363-d51f3f052b4f
md"""
Effect sizes can be classified as Medium, Insignificant, Small, and Large, respectively.
"""

# ╔═╡ b6c33e70-55c1-11eb-3892-d321937a86ed
md"""
#### Check presence of interactions
"""

# ╔═╡ e3220a00-55c1-11eb-080c-451cd462a4d6
plot(medianresult)

# ╔═╡ cbfc9a20-55d0-11eb-3cd3-23a3413a033c
md"""
The lines are quite different, suggesting a potential strong interaction. However, the F statistic for the interaction was not significant. Perhaps the number of replicates isn't enough due to variance between replicates. May be worth checking simple effects anyway.
"""

# ╔═╡ 2ed8ff90-568e-11eb-10a3-0b97b2425d2f
md"""
#### Interpret main effects
Nutlin-3a is just above the significance threshold with a medium effect size while Doxycyline produced no changes whatsoever. Nutlin-3a was thus associated with a drop in median distance, likely due to an increase in expression level.
"""

# ╔═╡ 76abe9e0-6af2-11eb-0dc5-4594d89c209b
md"""
### Figure
"""

# ╔═╡ 81b932c0-6af2-11eb-0a5c-09c0b3f10519
let
	import StatsPlots.mm
	p53_mdm2_median = [medianmeasurementsw[:,:,1,1,1] |> vec; medianmeasurementsw[:,:,2,1,1] |> vec; medianmeasurementsw[:,:,1,2,1] |> vec; medianmeasurementsw[:,:,2,2,1] |> vec]
	groups = repeat([2, 1, 4, 3], inner = 30)
	boxplot(groups, p53_mdm2_median, outliers=false,
		    label=["- Dox", "+ Dox"],
			legend=:none,
			seriescolor=:white,
			size=(256, 512),
			line=(2, 0.75),
		    xgrid=:none,
			xaxis=("Treatment", (1:4, ["-Dox\n-Nut", "+Dox\n-Nut", "-Dox\n+Nut", "+Dox\n+Nut"])),
			yaxis=("Median distance (nm)", [0, 600]))
	p53_mdm2_median_means = dropdims(mean(medianmeasurementsw[:,:,:,:,1], dims = 1), dims = 1) |> vec
	dotplot!(repeat([2, 1, 4, 3], inner=3), p53_mdm2_median_means |> vec, mode = :none, label="", marker=(4, 0.75, :rect, repeat([:orange, :darkblue, :darkred]), stroke(0)))
	dotplot!(groups, p53_mdm2_median, mode = :density, label="", marker=(2, 0.5, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
end

# ╔═╡ b69b87b2-6af5-11eb-30c2-bfa8159a8e0d
md"""Create paper-quality version and save it:"""

# ╔═╡ deafcdc0-6af4-11eb-1673-0101d277b4b3
let
	import StatsPlots.mm
	p53_mdm2_median = [medianmeasurementsw[:,:,1,1,1] |> vec; medianmeasurementsw[:,:,2,1,1] |> vec; medianmeasurementsw[:,:,1,2,1] |> vec; medianmeasurementsw[:,:,2,2,1] |> vec]
	groups = repeat([2, 1, 4, 3], inner = 30)
	boxplot(groups, p53_mdm2_median, outliers=false,
			legend=:none,
		    left_margin=25mm,
			seriescolor=:white,
		    guidefontsize=52,
            tickfontsize=48,
			size=(1024, 2048),
			line=(8, 0.75),
		    xgrid=:none,
			xaxis=("Treatment", (1:4, ["-Dox\n-Nut", "+Dox\n-Nut", "-Dox\n+Nut", "+Dox\n+Nut"])),
			yaxis=("Median distance (nm)", [0, 600]))
	p53_mdm2_median_means = dropdims(mean(medianmeasurementsw[:,:,:,:,1], dims=1), dims=1) |> vec
	dotplot!(repeat([2, 1, 4, 3], inner=3), p53_mdm2_median_means |> vec, mode = :none, label="", marker=(12, 0.75, :rect, repeat([:orange, :darkblue, :darkred]), stroke(0)))
	dotplot!(groups, p53_mdm2_median, mode = :density, label="", marker=(8, 0.5, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
	savefig(joinpath(outputdir, "p53-mdm2-medians.png"))
end

# ╔═╡ 25d0ce1e-5691-11eb-2cc8-49980b6b2161
md"""
### Monte Carlo fraction associated
#### Check unusual cases
"""

# ╔═╡ 8c0ebd10-56a9-11eb-3e8c-95b53bf12061
let
	p1 = boxplot(montecarlomeasurements[:, :, 1, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p2 = boxplot(montecarlomeasurements[:, :, 2, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p3 = boxplot(montecarlomeasurements[:, :, 1, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p4 = boxplot(montecarlomeasurements[:, :, 2, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))

	plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "Mdm2-p53")
end

# ╔═╡ faa812c0-56aa-11eb-3b6b-99fe5f200b67
md"""
Potential outliers in box plots:
"""

# ╔═╡ 012a7c50-56ab-11eb-3d50-c9820925f9ac
[montecarlomeasurements[6,1,2,1,1],
 montecarlomeasurements[1,3,1,2,1],
 montecarlomeasurements[5,3,2,2,1]]

# ╔═╡ 1637a040-56ac-11eb-2433-49c37383e005
md"""Evaluating z-scores within each replicate:"""

# ╔═╡ 87ec3520-56ac-11eb-2a1a-2981155d4168
montecarloflat = [montecarlomeasurements[:, :, 1, 1, 1] montecarlomeasurements[:, :, 2, 1, 1] montecarlomeasurements[:, :, 1, 2, 1] montecarlomeasurements[:, :, 2, 2, 1]]

# ╔═╡ 9ce1f5f0-56ac-11eb-22a5-d9b5d6a18978
montecarlozscores = [zscore(montecarloflat[:,i]) for i ∈ axes(montecarloflat, 2)]

# ╔═╡ ba7d2760-56ac-11eb-11ca-d543a628c68b
md"""
- In [1,1] Replicate 1: All within 1.96
- In [1,1] Replicate 2: All within 1.96
- In [1,1] Replicate 3: All within 1.96
- In [1,2] Replicate 1: 1 value, Cell 6 (10%) outside 1.96 significantly. None outside 2.58.
- In [1,2] Replicate 2: All within 1.96
- In [1,2] Replicate 3: All within 1.96
- In [2,1] Replicate 1: All within 1.96
- In [2,1] Replicate 2: All within 1.96
- In [2,1] Replicate 3: All within 1.96
- In [2,2] Replicate 1: 1 value, Cell 7 (10%) outside 1.96. None outside 2.58.
- In [2,2] Replicate 2: 1 value, Cell 9 (10%) outside 1.96. None outside 2.58.
- In [2,2] Replicate 3: All within 1.96
"""

# ╔═╡ c18cc050-56ad-11eb-261c-4177aa31ca7c
md"""
None of these potential outliers are due to data entry error or measurement error. Sampling a different population is possible but unknown. These potential outliers are relatively modest. In addition, since an ANOVA will be used, must be careful not to delete points from a balanced design. Therefore, no outliers will be removed. However, I shall Windsorize the significant one in [1,2] Replicate 1.
"""

# ╔═╡ 935aa530-6d73-11eb-2012-4d5e8b594f42
begin
	montecarlomeasurementsw = copy(montecarlomeasurements)
	montecarloflatw = copy(montecarloflat)
	montecarlomeasurementsw[6,1,2,1,1] = montecarlomeasurementsw[2,1,2,1,1]
	montecarloflatw[6,4] = montecarloflat[2,4]
end;

# ╔═╡ 2bba9b02-56ae-11eb-2e30-439dc20273f8
md"""
#### Check assumption of homoscedasticity
"""

# ╔═╡ 0f8a00a0-56af-11eb-2dc6-abd2a45298a4
md"""
Based on the ZResid/ZPred plot, the data appears fairly heteroscedastic.

Conducted Levene's test, a 1-way ANOVA of the absolute deviation between each value and the group mean, by replicate.
"""

# ╔═╡ b1ae0750-56af-11eb-0714-35b78e197b5e
levene(montecarloflatw)

# ╔═╡ f9676320-56af-11eb-1481-ef9502edded8
md"""
By Levene's test, the data is likely heteroscedastic.

Group sizes are equal, but samples are on the smaller side. CIs and significance tests may be affected.
"""

# ╔═╡ 4e36f730-56b0-11eb-21e8-0d2e7d0a2562
md"""
#### Check assumption of normality
"""

# ╔═╡ 4423b3d0-56b2-11eb-22be-fd86deabef7f
md"""
There are a couple notable deviations but for the *most* part they appear close to normal. [2,2] replicate 1 is the most deviant.
"""

# ╔═╡ db6922c0-56b2-11eb-2601-5525d3f6bf70
montecarloskewness = [skewness(montecarloflatw[:,i]) for i ∈ axes(montecarloflatw, 2)]

# ╔═╡ 047df370-56b3-11eb-3036-2393031bd967
montecarlokurtosis = [kurtosis(montecarloflatw[:,i]) for i ∈ axes(montecarloflatw, 2)]

# ╔═╡ 3316f790-56b3-11eb-3d9c-e9ebf797c923
md"""
From skewness and kurtosis statistics, none was greater in value than 1.96 (5% level).
"""

# ╔═╡ 6f1da2c0-56b3-11eb-2506-5f9a944735b6
md"""
#### Conduct ANOVA tests
"""

# ╔═╡ 755d8510-56b3-11eb-3a0c-53b643ba2053
montecarloresult = anova(
    permutedims(montecarlomeasurementsw[:, :, :, :, 1], (1,3,4,2)),
    [fixed, fixed, subject],
    factornames = ["Doxycycline", "Nutlin-3a", "Replicate"],
)

# ╔═╡ bf915350-574e-11eb-2f66-236c4ec9e600
md"""
Effect sizes can be classified as Large, Small, Insignificant, and Large, respectively.
"""

# ╔═╡ d79abcc0-574e-11eb-3777-4f8b85036436
md"""
#### Check presence of interactions
"""

# ╔═╡ f5f23950-574e-11eb-02b5-2fb9b6272ce6
plot(montecarloresult)

# ╔═╡ 27747350-5752-11eb-12b9-dfffc84b15c8
md"""
Lines are nearly parallel and the interaction F statistic was solidly non-significant. Checking the simple effects is valid.
"""

# ╔═╡ 248d1560-5753-11eb-2156-89cfb6742205
md"""
#### Interpret main effects
Nutlin-3a is below the significance threshold with a large effect size while Doxycyline produced a small effect with low probability. Nutlin-3a was thus associated with a decrease in fraction associated, while MEG3 expression likely was not associated with a change.
"""

# ╔═╡ 552d3550-6af5-11eb-2fcb-fb755145a298
md"""
### Figure
"""

# ╔═╡ 5ab31670-6af5-11eb-2e08-83b5228005d1
let
	import StatsPlots.mm
	p53_mdm2_montecarlo = [montecarlomeasurementsw[:,:,1,1,1] |> vec; montecarlomeasurementsw[:,:,2,1,1] |> vec; montecarlomeasurementsw[:,:,1,2,1] |> vec; montecarlomeasurementsw[:,:,2,2,1] |> vec]
	groups = repeat([2, 1, 4, 3], inner = 30)
	boxplot(groups, p53_mdm2_montecarlo, outliers=false,
		    label=["- Dox", "+ Dox"],
			legend=:none,
			seriescolor=:white,
		    left_margin=2mm,
			size=(256, 512),
			line=(2, 0.75),
		    xgrid=:none,
			xaxis=("Treatment", (1:4, ["-Dox\n-Nut", "+Dox\n-Nut", "-Dox\n+Nut", "+Dox\n+Nut"])),
			yaxis=("Median distance (nm)", [0, 0.25]))
	p53_mdm2_montecarlo_means = dropdims(mean(montecarlomeasurementsw[:,:,:,:,1], dims = 1), dims = 1) |> vec
	dotplot!(repeat([2, 1, 4, 3], inner=3), p53_mdm2_montecarlo_means |> vec, mode = :none, label="", marker=(4, 0.75, :rect, repeat([:orange, :darkblue, :darkred]), stroke(0)))
	dotplot!(groups, p53_mdm2_montecarlo, mode = :density, label="", marker=(2, 0.5, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
end

# ╔═╡ bc4851c0-6af5-11eb-0db5-cb71d67f93e3
md"""Create paper-quality version and save it:"""

# ╔═╡ af04a0e0-6af5-11eb-05ed-99c27fd1a0ef
let
	import StatsPlots.mm
	p53_mdm2_montecarlo = [montecarlomeasurementsw[:,:,1,1,1] |> vec; montecarlomeasurementsw[:,:,2,1,1] |> vec; montecarlomeasurementsw[:,:,1,2,1] |> vec; montecarlomeasurementsw[:,:,2,2,1] |> vec]
	groups = repeat([2, 1, 4, 3], inner = 30)
	boxplot(groups, p53_mdm2_montecarlo, outliers=false,
		    label=["- Dox", "+ Dox"],
			legend=:none,
			seriescolor=:white,
		    left_margin=25mm,
		    guidefontsize=52,
            tickfontsize=48,
			size=(1024, 2048),
			line=(8, 0.75),
		    xgrid=:none,
			xaxis=("Treatment", (1:4, ["-Dox\n-Nut", "+Dox\n-Nut", "-Dox\n+Nut", "+Dox\n+Nut"])),
			yaxis=("Fraction associated", [0, 0.25]))
	p53_mdm2_montecarlo_means = dropdims(mean(montecarlomeasurementsw[:,:,:,:,1], dims=1), dims=1) |> vec
	dotplot!(repeat([2, 1, 4, 3], inner=3), p53_mdm2_montecarlo_means |> vec, mode = :none, label="", marker=(12, 0.75, :rect, repeat([:orange, :darkblue, :darkred]), stroke(0)))
	dotplot!(groups, p53_mdm2_montecarlo, mode = :density, label="", marker=(8, 0.5, repeat([:orange, :darkblue, :darkred], inner=10), stroke(0)))
	savefig(joinpath(outputdir, "p53-mdm2-montecarlo.png"))
end

# ╔═╡ 78883b10-5756-11eb-32cf-b95fdeebda79
md"""
### Normalized Monte Carlo fraction associated
This approach was abandoned because we determined it was over-normalizing the data, but still included here for completeness.

#### Check unusual cases
"""

# ╔═╡ df7dce20-5756-11eb-3796-037e8e3c49f9
let
	p1 = boxplot(normalizedmontecarlomeasurements[:, :, 1, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p2 = boxplot(normalizedmontecarlomeasurements[:, :, 2, 1, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p3 = boxplot(normalizedmontecarlomeasurements[:, :, 1, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))
	p4 = boxplot(normalizedmontecarlomeasurements[:, :, 2, 2, 1], xaxis = ("Replicates", [1, 2, 3]), yaxis = ("Fraction bound"))

	plot(p1, p2, p3, p4, layout = grid(2, 2), legend = :none, plot_title = "Mdm2-p53")
end

# ╔═╡ faf5aa60-5756-11eb-3b5b-f1ad5e004af2
md"""
Potential outliers in box plots:
"""

# ╔═╡ ff845a40-5756-11eb-2ce8-d5ba8458e2ae
[normalizedmontecarlomeasurements[6,1,2,1,1],
 normalizedmontecarlomeasurements[[2,9],3,1,2,1]]

# ╔═╡ d2282670-5757-11eb-0bd9-b1aaf0b37d70
md"""Evaluating z-scores within each replicate:"""

# ╔═╡ d99fec80-5757-11eb-39cd-4b0febbf61e3
normmontecarloflat = [normalizedmontecarlomeasurements[:, :, 1, 1, 1] normalizedmontecarlomeasurements[:, :, 2, 1, 1] normalizedmontecarlomeasurements[:, :, 1, 2, 1] normalizedmontecarlomeasurements[:, :, 2, 2, 1]]

# ╔═╡ edee22b2-5757-11eb-317a-d3e4d31486cd
normmontecarlozscores = [zscore(normmontecarloflat[:,i]) for i ∈ axes(normmontecarloflat, 2)]

# ╔═╡ a05efb90-5758-11eb-10ca-fb9fb8e05821
md"""
- In [1,1] Replicate 1: All within 1.96
- In [1,1] Replicate 2: All within 1.96
- In [1,1] Replicate 3: 1 value, Cell 2 (10%) outside 1.96 barely. None outside 2.58
- In [1,2] Replicate 1: 1 value, Cell 6 (10%) outside 1.96 significantly. None outside 2.58.
- In [1,2] Replicate 2: All within 1.96
- In [1,2] Replicate 3: All within 1.96
- In [2,1] Replicate 1: All within 1.96
- In [2,1] Replicate 2: All within 1.96
- In [2,1] Replicate 3: 1 value, Cell 8 (10%) outside 1.96. None outside 2.58.
- In [2,2] Replicate 1: All within 1.96
- In [2,2] Replicate 2: All within 1.96
- In [2,2] Replicate 3: All within 1.96
"""

# ╔═╡ 7323e450-5759-11eb-1c62-7ba711bf418b
md"""
None of these potential outliers are due to data entry error or measurement error. Sampling a different population is possible but unknown. These potential outliers are relatively modest. In addition, since an ANOVA will be used, must be careful not to delete points from a balanced design. Therefore, no outliers will be removed. However, I will Windsorize the extreme outlier in [1,2] Replicate 1.
"""

# ╔═╡ 9ecc38b0-6d74-11eb-2d1a-ebfa12f11a64
begin
	normalizedmontecarlomeasurementsw = copy(normalizedmontecarlomeasurements)
	normmontecarloflatw = copy(normmontecarloflat)
	normalizedmontecarlomeasurementsw[6,1,2,1,1] = normalizedmontecarlomeasurementsw[4,1,2,1,1]
	normmontecarloflatw[6,4] = normmontecarloflat[4,4]
end;

# ╔═╡ df4eee90-5759-11eb-02ef-53a0c6f31f07
md"""
#### Check assumption of homoscedasticity
"""

# ╔═╡ 8c2a8790-575b-11eb-34be-99f584826d43
md"""
Based on the ZResid/ZPred plot, the data appears heteroscedastic, although not much relationship between residual and predicted.

Conducted Levene's test, a 1-way ANOVA of the absolute deviation between each value and the group mean, by replicate.
"""

# ╔═╡ 9aade460-575b-11eb-29da-4ff6881ceb60
levene(normmontecarloflatw)

# ╔═╡ b8f00430-575b-11eb-0ada-7f0001326feb
md"""
By Levene's test, the data is likely heteroscedastic.

Group sizes are equal, but samples are on the smaller side. CIs and significance tests may be affected.
"""

# ╔═╡ 1326cdce-575c-11eb-0f72-cf6eebe27bdf
md"""
#### Check assumption of normality
"""

# ╔═╡ 3c6921b0-575d-11eb-160f-6f06db5965f3
md"""
There are a couple notable deviations but for the *most* part they appear close to normal.
"""

# ╔═╡ 59314bb0-575d-11eb-2b54-e1bf609e6d36
normmontecarloskewness = [skewness(normmontecarloflat[:,i]) for i ∈ axes(normmontecarloflat, 2)]

# ╔═╡ 682c7090-575d-11eb-3b12-d50d450d3fce
normmontecarlokurtosis = [kurtosis(normmontecarloflat[:,i]) for i ∈ axes(normmontecarloflat, 2)]

# ╔═╡ 80f22522-575d-11eb-211b-59fe55838280
md"""
From skewness and kurtosis statistics, none was greater in value than 1.96 (5% level).
"""

# ╔═╡ 91c676d0-575d-11eb-2fd5-97ca14bfdec0
md"""
#### Conduct ANOVA tests
"""

# ╔═╡ b7624d5e-575d-11eb-1221-974537c3b48e
normmontecarloresult = anova(
    permutedims(normalizedmontecarlomeasurementsw[:, :, :, :, 1], (1,3,4,2)),
    [fixed, fixed, subject],
    factornames = ["Doxycycline", "Nutlin-3a", "Replicate"],
)

# ╔═╡ 8f0dad80-575f-11eb-2312-adae92d50e83
md"""
### Check presence of interactions
"""

# ╔═╡ a5489880-575f-11eb-1afa-c3bfaf5dccc8
plot(normmontecarloresult)

# ╔═╡ 4b211122-575f-11eb-0aa9-59fc4e9a2df7
md"""
The lines cross, but the values are so small that it's likely nothing. F-statistic is not significant. Interpreting simple effects valid.
"""

# ╔═╡ 408e5d30-5764-11eb-124f-4f2630d1d908
md"""
#### Interpret main effects
Nutlin-3a did not have a significant effect on fraction bound (from 1.7% to 1.7%). This is unexpected for an inhibitor of Mdm2-p53 binding, but it may reflect that p53 and Mdm2 levels both increase until a steady state is reached, maintaining a similar level of association.

Doxycycline did not have a significant effect, though with a small effect size, though may be due to lack of power in the experiment. There was still a consistent increase with doxycycline exposure (from 1.4% to 1.9%). Thus, MEG3 overexpression does not cause significantly less binding between p53 and Mdm2.
"""

# ╔═╡ d83ab120-6be5-11eb-0b44-176f26137172
md"""
## Localization/molecule maps
"""

# ╔═╡ e5020610-6be5-11eb-25ab-adb9479d6b16
let
	# Exp 2 (p53-Mdm2)
	# A8
	insetx, insety = [15100, 16500], [30000, 31400]
	r = experimentresults[1][3][1][8]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "2 - A8 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "2 - A8 dSTORM 1400nm.png"))
	p1 = localizationsplot(r, insetbox = [insetx, insety], include_scalebar = true)
	p2 = insetplot(r, insetx, insety, include_scalebar = true)
	plot(p1, p2, layout = grid(1,2), size=(1024, 512), fmt = :png)
end

# ╔═╡ 840e9e20-6be7-11eb-2b80-37ffe2f4495d
let
	# Exp 2 (p53-Mdm2)
	# B5
	insetx, insety = [16050, 17450], [21100, 22500]
	r = experimentresults[1][3][2][5]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "2 - B5 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "2 - B5 dSTORM 1400nm.png"))
	p1 = localizationsplot(r, insetbox = [insetx, insety], include_scalebar = true)
	p2 = insetplot(r, insetx, insety, include_scalebar = true)
	plot(p1, p2, layout = grid(1,2), size=(1024, 512), fmt = :png)
end

# ╔═╡ 7affa0d0-6be8-11eb-3395-555814cab9be
let
	# Exp 2 (p53-Mdm2)
	# C9
	insetx, insety = [16400, 17800], [22697, 24096]
	r = experimentresults[1][2][3][9]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "2 - C9 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "2 - C9 dSTORM 1400nm.png"))
	p1 = localizationsplot(r, insetbox = [insetx, insety], include_scalebar = true)
	p2 = insetplot(r, insetx, insety, include_scalebar = true)
	plot(p1, p2, layout = grid(1,2), size=(1024, 512), fmt = :png)
end

# ╔═╡ 13ecfdf0-6bea-11eb-3af4-8d627cccdd2b
let
	# Exp 2 (p53-Mdm2)
	# D4
	insetx, insety = [24100, 25500], [14000, 15400]
	r = experimentresults[1][2][4][4]
	localizationsplot(r, insetbox = [insetx, insety], forprint = true)
	savefig(joinpath(outputdir, "2 - D4 dSTORM.png"))
	insetplot(r, insetx, insety, include_scalebar = false, forprint = true)
	savefig(joinpath(outputdir, "2 - D4 dSTORM 1400nm.png"))
	p1 = localizationsplot(r, insetbox = [insetx, insety], include_scalebar = true)
	p2 = insetplot(r, insetx, insety, include_scalebar = true)
	plot(p1, p2, layout = grid(1,2), size=(1024,512), fmt = :png)
end

# ╔═╡ 5db1f140-6ce7-11eb-29f1-ffb93c61c906
let
	# Plot all
	r = localizationsplot.(experimentresults[1][j][i][k] for k ∈ 1:10 for i ∈ 1:4 for j ∈ 1:3)
	p = plot(r..., size=(2048, 2048), layout=grid(10,12), fmt=:png)
	savefig(joinpath(outputdir, "p53-mdm2-all-localizations.png"))
	p
end

# ╔═╡ 9c92b290-56ae-11eb-2595-85845fe95f0e
md"""
## Functions
"""

# ╔═╡ 2623d480-56af-11eb-2ab8-0d0afd8c45dc
function zresid_zpred(data)
	z = zscore(data)
	zpred = repeat(mean(z, dims = 1), 10)
	zresid = z .- zpred
	zresid, zpred
end

# ╔═╡ a655a120-56ae-11eb-30be-459d7c3067dc
function zresid_zpred_plot(data)
	zresid, zpred = zresid_zpred(data)
	scatter(
		zresid,
		zpred,
		xaxis = ("Standardized Residual (ZResid)"),
		yaxis = ("Standardized Predicted Value (ZPred)"),
		legend = :none,
	)
end

# ╔═╡ c8f83d60-5518-11eb-2643-a30df2dd9fb9
zresid_zpred_plot(mediansflatw)

# ╔═╡ 3a3fa580-56ae-11eb-1f19-f3c06ef2d9d7
zresid_zpred_plot(montecarloflatw)

# ╔═╡ 59522180-575a-11eb-39c4-f9eaa27f7901
zresid_zpred_plot(normmontecarloflatw)

# ╔═╡ 7922d1d0-56b0-11eb-261c-5f87ca587394
function qqnormplot(data)
	zresid, _ = zresid_zpred(data)
	p = [qqnorm(zresid[:, i]) for i ∈ 1:12]
	plot(p..., layout = grid(4, 3), legend = :none, plot_title = "MDM2-p53 qqnorm")
end

# ╔═╡ 6f727ed0-5519-11eb-1462-01c83aa975b6
qqnormplot(mediansflatw)

# ╔═╡ 5ef3f050-56b0-11eb-31f9-419a4161b67e
qqnormplot(montecarloflat)

# ╔═╡ 1823a8d0-575c-11eb-0fa1-893a914f69f7
qqnormplot(normmontecarloflat)

# ╔═╡ Cell order:
# ╟─56dcb910-5507-11eb-0b0a-a1c2aef89207
# ╟─90d2acb0-5507-11eb-39c4-df93c258dd98
# ╠═b11fdebe-5507-11eb-1001-537a20e51d0d
# ╟─86502ff0-550d-11eb-1c13-553251c061ee
# ╠═5e3db410-550d-11eb-0331-49a447464b07
# ╟─294e6a02-550e-11eb-3450-d9d5989ec08c
# ╟─9a94f370-5511-11eb-1f60-7fe83ed11def
# ╟─572ba5f2-5513-11eb-1937-3f03cec8252e
# ╟─864b181e-5513-11eb-0304-8983d37a8170
# ╟─27e24820-5514-11eb-1415-8711a194892d
# ╠═f0393540-5514-11eb-2f76-3514cc3bc85c
# ╟─f02df530-5515-11eb-0d9e-8dbb63ed2023
# ╟─1bf98170-5516-11eb-0c24-77ac43672fba
# ╟─25253910-5516-11eb-0a7a-1d73ed51d4a1
# ╟─d8c96c5e-5517-11eb-0861-056d5e425f04
# ╟─87eb8f70-5518-11eb-1639-f91b4bcb5c9b
# ╠═9b7f65a2-5952-11eb-0c21-733df3703bf9
# ╟─a87aac80-5518-11eb-38d4-5fe90140b820
# ╟─c8f83d60-5518-11eb-2643-a30df2dd9fb9
# ╟─f6e23cd0-5518-11eb-226a-c5751eb65b34
# ╟─0ffdd990-5519-11eb-376b-33be0cc966d7
# ╟─28a8da30-5519-11eb-2fb0-efe693c6e44c
# ╟─5f8de452-5519-11eb-0a2c-051d3e701dd5
# ╟─6f727ed0-5519-11eb-1462-01c83aa975b6
# ╟─c4dbef60-55bd-11eb-1822-49c3473000e8
# ╟─cf59a26e-55bd-11eb-1c4b-316a603b9007
# ╟─135dbe6e-55be-11eb-3a89-dd3d02c3ad4a
# ╟─6778a880-55be-11eb-1016-c59335beef0a
# ╟─9bf9aaf0-55be-11eb-17a2-75b3eec550f8
# ╟─265006c0-55c1-11eb-2ba8-3d628d096704
# ╠═5f1ce3b0-55c1-11eb-3d03-b53516cfc976
# ╟─756bf2ee-55c1-11eb-2363-d51f3f052b4f
# ╟─b6c33e70-55c1-11eb-3892-d321937a86ed
# ╟─e3220a00-55c1-11eb-080c-451cd462a4d6
# ╟─cbfc9a20-55d0-11eb-3cd3-23a3413a033c
# ╟─2ed8ff90-568e-11eb-10a3-0b97b2425d2f
# ╟─76abe9e0-6af2-11eb-0dc5-4594d89c209b
# ╟─81b932c0-6af2-11eb-0a5c-09c0b3f10519
# ╟─b69b87b2-6af5-11eb-30c2-bfa8159a8e0d
# ╠═deafcdc0-6af4-11eb-1673-0101d277b4b3
# ╟─25d0ce1e-5691-11eb-2cc8-49980b6b2161
# ╟─8c0ebd10-56a9-11eb-3e8c-95b53bf12061
# ╟─faa812c0-56aa-11eb-3b6b-99fe5f200b67
# ╠═012a7c50-56ab-11eb-3d50-c9820925f9ac
# ╟─1637a040-56ac-11eb-2433-49c37383e005
# ╟─87ec3520-56ac-11eb-2a1a-2981155d4168
# ╟─9ce1f5f0-56ac-11eb-22a5-d9b5d6a18978
# ╟─ba7d2760-56ac-11eb-11ca-d543a628c68b
# ╟─c18cc050-56ad-11eb-261c-4177aa31ca7c
# ╠═935aa530-6d73-11eb-2012-4d5e8b594f42
# ╟─2bba9b02-56ae-11eb-2e30-439dc20273f8
# ╟─3a3fa580-56ae-11eb-1f19-f3c06ef2d9d7
# ╟─0f8a00a0-56af-11eb-2dc6-abd2a45298a4
# ╟─b1ae0750-56af-11eb-0714-35b78e197b5e
# ╟─f9676320-56af-11eb-1481-ef9502edded8
# ╟─4e36f730-56b0-11eb-21e8-0d2e7d0a2562
# ╟─5ef3f050-56b0-11eb-31f9-419a4161b67e
# ╟─4423b3d0-56b2-11eb-22be-fd86deabef7f
# ╟─db6922c0-56b2-11eb-2601-5525d3f6bf70
# ╟─047df370-56b3-11eb-3036-2393031bd967
# ╟─3316f790-56b3-11eb-3d9c-e9ebf797c923
# ╟─6f1da2c0-56b3-11eb-2506-5f9a944735b6
# ╠═755d8510-56b3-11eb-3a0c-53b643ba2053
# ╟─bf915350-574e-11eb-2f66-236c4ec9e600
# ╟─d79abcc0-574e-11eb-3777-4f8b85036436
# ╟─f5f23950-574e-11eb-02b5-2fb9b6272ce6
# ╟─27747350-5752-11eb-12b9-dfffc84b15c8
# ╟─248d1560-5753-11eb-2156-89cfb6742205
# ╟─552d3550-6af5-11eb-2fcb-fb755145a298
# ╟─5ab31670-6af5-11eb-2e08-83b5228005d1
# ╟─bc4851c0-6af5-11eb-0db5-cb71d67f93e3
# ╠═af04a0e0-6af5-11eb-05ed-99c27fd1a0ef
# ╟─78883b10-5756-11eb-32cf-b95fdeebda79
# ╟─df7dce20-5756-11eb-3796-037e8e3c49f9
# ╟─faf5aa60-5756-11eb-3b5b-f1ad5e004af2
# ╠═ff845a40-5756-11eb-2ce8-d5ba8458e2ae
# ╟─d2282670-5757-11eb-0bd9-b1aaf0b37d70
# ╟─d99fec80-5757-11eb-39cd-4b0febbf61e3
# ╟─edee22b2-5757-11eb-317a-d3e4d31486cd
# ╟─a05efb90-5758-11eb-10ca-fb9fb8e05821
# ╟─7323e450-5759-11eb-1c62-7ba711bf418b
# ╠═9ecc38b0-6d74-11eb-2d1a-ebfa12f11a64
# ╟─df4eee90-5759-11eb-02ef-53a0c6f31f07
# ╟─59522180-575a-11eb-39c4-f9eaa27f7901
# ╟─8c2a8790-575b-11eb-34be-99f584826d43
# ╟─9aade460-575b-11eb-29da-4ff6881ceb60
# ╟─b8f00430-575b-11eb-0ada-7f0001326feb
# ╟─1326cdce-575c-11eb-0f72-cf6eebe27bdf
# ╟─1823a8d0-575c-11eb-0fa1-893a914f69f7
# ╟─3c6921b0-575d-11eb-160f-6f06db5965f3
# ╟─59314bb0-575d-11eb-2b54-e1bf609e6d36
# ╟─682c7090-575d-11eb-3b12-d50d450d3fce
# ╟─80f22522-575d-11eb-211b-59fe55838280
# ╟─91c676d0-575d-11eb-2fd5-97ca14bfdec0
# ╟─b7624d5e-575d-11eb-1221-974537c3b48e
# ╟─8f0dad80-575f-11eb-2312-adae92d50e83
# ╟─a5489880-575f-11eb-1afa-c3bfaf5dccc8
# ╟─4b211122-575f-11eb-0aa9-59fc4e9a2df7
# ╟─408e5d30-5764-11eb-124f-4f2630d1d908
# ╟─d83ab120-6be5-11eb-0b44-176f26137172
# ╠═e5020610-6be5-11eb-25ab-adb9479d6b16
# ╠═840e9e20-6be7-11eb-2b80-37ffe2f4495d
# ╠═7affa0d0-6be8-11eb-3395-555814cab9be
# ╠═13ecfdf0-6bea-11eb-3af4-8d627cccdd2b
# ╠═5db1f140-6ce7-11eb-29f1-ffb93c61c906
# ╟─9c92b290-56ae-11eb-2595-85845fe95f0e
# ╟─2623d480-56af-11eb-2ab8-0d0afd8c45dc
# ╟─a655a120-56ae-11eb-30be-459d7c3067dc
# ╟─7922d1d0-56b0-11eb-261c-5f87ca587394


text:
using Color
using FixedPointNumbers
using HDF5, JLD

require("interp.jl")

rawdata{T<:FixedPoint}(img::Array{Gray{T}, 2}) = convert(Array{Float64, 2}, img)
rawdata(img::Image) = rawdata(data(img)')

# Create dataset consiting of nonzero pixels of face images
# Returns:
#   dat : Matrix(# of nonzeros x # of filenames) - data matrix
#   nonzero_idxs : Vector((i, j)) - nonzero pixel coordinates
function facedata(datadir=CK_DATA_DIR_FACE_ALIGNED, imsize=IMSIZE)
    filenames = readdir(datadir)
    refmat = rawdata(imread(datadir * "/" * filenames[1]))
    refmat = imresize(refmat, imsize)
    # refmat = refmat ./ 256  -- imresize already does it
    nonzero_idxs = nonzero_indexes(refmat)
    dat = zeros(length(nonzero_idxs), length(filenames))
    for (i, fname) in enumerate(filenames)
        if i % 500 == 0 println(i, " ", fname) end
        mat = rawdata(imread(datadir * "/" * fname))
        mat = imresize(mat, imsize)
        # mat = mat ./ 256
        dat[:, i] = collect_nonzeros(mat, nonzero_idxs)
    end
    return dat, nonzero_idxs
end

namebase(filename) = join(split(filename, ['_', '.'])[1:3], "_")
labname(nbase) = nbase * "_emotion.txt"
facename(nbase) = nbase * ".png"

function read_resize_image(filename, imsize)
    mat = rawdata(imread(filename))
    return imresize(mat, imsize)
end

function read_aligned(datadir=CK_DATA_DIR, imsize=IMSIZE; labeled_only=false)
    facedir = datadir * FACES_ALIGNED_SUBDIR
    labdir = datadir * LABELES_SUBDIR
    namebases = (labeled_only ?
                 map(namebase, readdir(labdir)) :
                 map(namebase, readdir(facedir)))
    refimgpath = facedir * "/" * readdir(facedir)[1]
    refmat = read_resize_image(refimgpath, imsize)
    nzs = nonzero_indexes(refmat)
    img_dat = zeros(length(nzs), length(namebases))
    lab_dat = zeros(Int8, length(namebases))
    for (i, nbase) in enumerate(namebases)
        if i % 500 == 0 println("$(i) images processed") end
        mat = read_resize_image(facedir * "/" * facename(nbase), imsize)
        img_dat[:, i] = collect_nonzeros(mat, nzs)
        if isfile(labdir * "/" * labname(nbase))
            lab_dat[i] = int(float(open(readall, labdir * "/" * labname(nbase))))
        else
            lab_dat[i] = -1
        end
    end
    return img_dat, lab_dat, nzs
end

max_stars_repo_name: dfdx/masque

max_issues_repo_name: dfdx/masque

max_forks_repo_name: dfdx/masque


text:
using Test
using BeastUtils
using BeastUtils.TreeUtils, BeastUtils.Simulation, BeastUtils.MatrixUtils,
      BeastUtils.SlowLikelihoods, BeastUtils.XMLConstructor, BeastUtils.Logs,
      BeastUtils.RunBeast
using PhyloNetworks, LinearAlgebra
import Random
Random.seed!(666)

const tol = 1e-10

N = 10
K = 2
P = 5

taxa = ["taxon_" * "$i" for i = 1:N]

tree = rtree(N, labels = taxa)
newick = writeTopology(tree)

## Testing Basic MBD Model
Σ = randpd(P)

tdm = TreeDiffusionModel(tree, Σ)
tsm = TraitSimulationModel(taxa, tdm)

data = simulate(tsm)

ll = loglikelihood(tsm, data)

# TODO: move above to separate file

## Testing MBD Models
Σ = randpd(P)
Γ = randpd(P)

tdm = TreeDiffusionModel(tree, Σ)
rvm = ResidualVarianceModel(Γ)
tsm = TraitSimulationModel(taxa, tdm, rvm)

data = simulate(tsm)

bx = make_residual_xml(data, taxa, newick)
set_mbd_precision(bx, inv(Σ))
set_residual_precision(bx, inv(Γ))

filename = "test_residual"
xml_path = "$filename.xml"
save_xml(xml_path, bx)
run_beast(xml_path)
log_path = "$filename.log"
col, likelihoods = get_log_match(log_path, "likelihood", burnin=0.0)

standardize_data!(data)
ll_slow = loglikelihood(tsm, data, rescale_tree = true)
ll_beast = likelihoods[1]

@test ll_slow ≈ ll_beast atol=tol

## Testing Factor Models
Σ = Diagonal(ones(K))
Λ = Diagonal(rand(P))
L = randn(K, P)

tdm = TreeDiffusionModel(tree, Σ)
lfm = LatentFactorModel(L, Λ)
tsm = TraitSimulationModel(taxa, tdm, lfm)

data = simulate(tsm)

bx = make_pfa_xml(data, taxa, newick, K, useHMC = false)
set_residual_precision(bx, inv(Λ))
set_loadings(bx, L)

filename = "test_factors"
xml_path = "$filename.xml"
save_xml(xml_path, bx)
run_beast(xml_path)
log_path = "$filename.log"
col, likelihoods = get_log_match(log_path, "likelihood", burnin=0.0)

standardize_data!(data)
ll_slow = loglikelihood(tsm, data, rescale_tree = true)
ll_beast = likelihoods[1]

@test ll_slow ≈ ll_beast atol=tol


text:
using Test
using Algorithm101
const ALGO = Algorithm101

@testset "Matrix Chain Multiplication" begin
    @testset "Utils" begin
        chain = [
            ALGO.Dim("A", 2, 6),
            ALGO.Dim("B", 6, 3),
            ALGO.Dim("C", 3, 5),
            ALGO.Dim("D", 5, 10),
            ALGO.Dim("E", 10, 2),
        ]

        @test ALGO.check_dim(chain)
    end

    @testset "MCM" begin
        chain = [
            ALGO.Dim("A", 30, 35),
            ALGO.Dim("B", 35, 15),
            ALGO.Dim("C", 15, 5),
            ALGO.Dim("D", 5, 10),
            ALGO.Dim("E", 10, 20),
            ALGO.Dim("F", 20, 25),
        ]
        @test ALGO.multiplice_matrix_chain(chain) == ("((A(BC))((DE)F))", 15125)
    end
end


text:
using Test
using ConvDiffMIPDECO
using jInv.Mesh

domain = [0 1. 0 1 0 1]
n  = [7 8 9]
M2D  = getRegularMesh(domain[1:4],n[1:2])

gx = X -> X[:,1]
gy = X -> X[:,2]
gz = X -> X[:,3]

xa,ya = getCellCenteredAxes(M2D)
(gl,gb,gr,gt) = getBoundaryCondition(gx, M2D)
@test all(gl.==0.0)
@test all(gb.==xa)
@test all(gr.==1.0)
@test all(gt.==xa)
@test length(gl)==n[2]
@test length(gb)==n[1]
@test length(gr)==n[2]
@test length(gt)==n[1]

(gl,gb,gr,gt) = getBoundaryCondition(gy, M2D)
@test all(gl.==ya)
@test all(gb.==0.0)
@test all(gr.==ya)
@test all(gt.==1.0)

M3D   = getRegularMesh(domain,n)
xa,ya,za = getCellCenteredAxes(M3D)
(gx1,gy1,gz1,gx2,gy2,gz2) = getBoundaryCondition(gx, M3D)
@test all(gx1.==0.0)
@test all(gx2.==1.0)
@test all(gy1.==repeat(xa,1,n[3]))
@test all(gy2.==repeat(xa,1,n[3]))
@test all(gz1.==repeat(xa,1,n[2]))
@test all(gz2.==repeat(xa,1,n[2]))

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
## Types

struct ProblemSpec1D
  ranges::Array{Tuple{String,Float64,Float64}}
  f::Function
end

## Polynomial definitions

# The 2 functions from Farzad
function poly1(x, y, z, w)
  -32 * w^3 * x * z^3
    + 44 * w^5 * x * z^3
    - 16 * w^7 * x * z^3
    + 4 * w^9 * x * z^3
    + 16 * w^3 * y * z^3
    - 100 * w^5 * y * z^3
    + 128 * w^7 * y * z^3
    - 44 * w^9 * y * z^3
    - 56 * w^5 * x * y * z^3
    + 48 * w^7 * x * y * z^3
    - 16 * w^9 * x * y * z^3
    - 56 * w^5 * y^2 * z^3
    + 48 * w^7 * y^2 * z^3
    - 16 * w^9 * y^2 * z^3
    + 32 * w^3 * x * z^5
    + 16 * w^7 * x * z^5
    - 64 * w^9 * x * z^5
    + 16 * w^11 * x * z^5
    + 16 * w^3 * x^2 * z^5
    - 64 * w^5 * x^2 * z^5
    + 96 * w^7 * x^2 * z^5
    - 64 * w^9 * x^2 * z^5
    + 16 * w^11 * x^2 * z^5
    - 16 * w^3 * y * z^5
    + 144 * w^5 * y * z^5
    - 128 * w^7 * y * z^5
    - 16 * w^9 * y * z^5
    + 16 * w^11 * y * z^5
    + 16 * w^3 * x * y * z^5
    + 48 * w^5 * x * y * z^5
    - 16 * w^7 * x * y * z^5
    + 80 * w^9 * x * y * z^5
    - 32 * w^11 * x * y * z^5
    + 112 * w^5 * y^2 * z^5
    - 112 * w^7 * y^2 * z^5
    + 144 * w^9 * y^2 * z^5
    - 48 * w^11 * y^2 * z^5
    - 44 * w^5 * x * z^7
    - 16 * w^7 * x * z^7
    + 80 * w^11 * x * z^7
    - 20 * w^13 * x * z^7
    + 8 * w^5 * x^2 * z^7
    - 32 * w^7 * x^2 * z^7
    + 48 * w^9 * x^2 * z^7
    - 32 * w^11 * x^2 * z^7
    + 8 * w^13 * x^2 * z^7
    - 44 * w^5 * y * z^7
    - 64 * w^7 * y * z^7
    + 144 * w^9 * y * z^7
    - 64 * w^11 * y * z^7
    + 28 * w^13 * y * z^7
    - 48 * w^5 * x * y * z^7
    + 48 * w^7 * x * y * z^7
    - 192 * w^9 * x * y * z^7
    + 48 * w^11 * x * y * z^7
    - 56 * w^5 * y^2 * z^7
    + 80 * w^7 * y^2 * z^7
    - 240 * w^9 * y^2 * z^7
    + 80 * w^11 * y^2 * z^7
    - 8 * w^13 * y^2 * z^7
    + 16 * w^7 * x * z^9
    + 64 * w^9 * x * z^9
    - 80 * w^11 * x * z^9
    + 64 * w^7 * y * z^9
    - 80 * w^9 * y * z^9
    + 64 * w^11 * y * z^9
    - 48 * w^13 * y * z^9
    - 16 * w^7 * x * y * z^9
    + 112 * w^9 * x * y * z^9
    - 16 * w^11 * x * y * z^9
    + 16 * w^13 * x * y * z^9
    - 16 * w^7 * y^2 * z^9
    + 112 * w^9 * y^2 * z^9
    - 16 * w^11 * y^2 * z^9
    + 16 * w^13 * y^2 * z^9
    - 4 * w^9 * x * z^11
    - 16 * w^11 * x * z^11
    + 20 * w^13 * x * z^11
    - 4 * w^9 * y * z^11
    - 16 * w^11 * y * z^11
    + 20 * w^13 * y * z^11
    - 16 * w^11 * x * y * z^11
    - 8 * w^13 * x * y * z^11
    - 16 * w^11 * y^2 * z^11
    - 8 * w^13 * y^2 * z^11
end

function poly2(x, y, z, w)
  20 * w * x * z^3
    - 48 * w^3 * x * z^3
    + 28 * w^5 * x * z^3
    + 8 * w * x^2 * z^3
    - 16 * w^3 * x^2 * z^3
    + 8 * w^5 * x^2 * z^3
    + 20 * w * y * z^3
    - 20 * w^5 * y * z^3
    + 8 * w * x * y * z^3
    - 16 * w^3 * x * y * z^3
    - 8 * w^5 * y^2 * z^3
    - 16 * w * x * z^5
    + 64 * w^3 * x * z^5
    - 64 * w^5 * x * z^5
    + 16 * w^7 * x * z^5
    + 16 * w * x^2 * z^5
    + 16 * w^3 * x^2 * z^5
    - 80 * w^5 * x^2 * z^5
    + 48 * w^7 * x^2 * z^5
    - 16 * w * y * z^5
    - 80 * w^3 * y * z^5
    + 80 * w^5 * y * z^5
    + 16 * w^7 * y * z^5
    + 16 * w * x * y * z^5
    + 16 * w^3 * x * y * z^5
    - 48 * w^5 * x * y * z^5
    + 32 * w^7 * x * y * z^5
    + 32 * w^5 * y^2 * z^5
    - 16 * w^7 * y^2 * z^5
    - 4 * w * x * z^7
    - 80 * w^3 * x * z^7
    + 144 * w^5 * x * z^7
    - 16 * w^7 * x * z^7
    - 44 * w^9 * x * z^7
    - 112 * w^3 * x^2 * z^7
    + 240 * w^5 * x^2 * z^7
    - 144 * w^7 * x^2 * z^7
    + 16 * w^9 * x^2 * z^7
    - 4 * w * y * z^7
    + 64 * w^3 * y * z^7
    - 64 * w^7 * y * z^7
    + 4 * w^9 * y * z^7
    - 112 * w^3 * x * y * z^7
    + 192 * w^5 * x * y * z^7
    - 80 * w^7 * x * y * z^7
    + 16 * w^9 * x * y * z^7
    - 48 * w^5 * y^2 * z^7
    + 64 * w^7 * y^2 * z^7
    + 64 * w^3 * x * z^9
    - 64 * w^5 * x * z^9
    - 128 * w^7 * x * z^9
    + 128 * w^9 * x * z^9
    + 16 * w^3 * x^2 * z^9
    - 80 * w^5 * x^2 * z^9
    + 112 * w^7 * x^2 * z^9
    - 48 * w^9 * x^2 * z^9
    + 16 * w^3 * y * z^9
    - 16 * w^5 * y * z^9
    + 16 * w^7 * y * z^9
    - 16 * w^9 * y * z^9
    + 16 * w^3 * x * y * z^9
    - 48 * w^5 * x * y * z^9
    + 16 * w^7 * x * y * z^9
    - 48 * w^9 * x * y * z^9
    + 32 * w^5 * y^2 * z^9
    - 96 * w^7 * y^2 * z^9
    - 44 * w^5 * x * z^11
    + 144 * w^7 * x * z^11
    - 100 * w^9 * x * z^11
    + 56 * w^5 * x^2 * z^11
    - 112 * w^7 * x^2 * z^11
    + 56 * w^9 * x^2 * z^11
    - 44 * w^5 * y * z^11
    + 44 * w^9 * y * z^11
    + 48 * w^5 * x * y * z^11
    - 48 * w^7 * x * y * z^11
    + 56 * w^9 * x * y * z^11
    - 8 * w^5 * y^2 * z^11
    + 64 * w^7 * y^2 * z^11
    - 16 * w^7 * x * z^13
    + 16 * w^9 * x * z^13
    + 32 * w^7 * y * z^13
    - 32 * w^9 * y * z^13
    - 16 * w^7 * x * y * z^13
    - 16 * w^7 * y^2 * z^13
end

function poly3(x)
  sum(x .* x)
end

poly1def = ProblemSpec1D(
  [("x", -1., 1.), ("y", -1., 1.), ("z", -1., 1.), ("w", -1., 1.)],
  x -> poly1(x[1], x[2], x[3], x[4]))

poly4def = ProblemSpec1D(
  [("x", -1., 1.), ("y", -1., 1.)],
  x -> poly1(x[1], x[2], 2000., 1000.))

poly2def = ProblemSpec1D(
  [("x", -1., 1.), ("y", -1., 1.), ("z", -1., 1.), ("w", -1., 1.)],
  x -> poly2(x[1], x[2], x[3], x[4]))

poly3def = ProblemSpec1D(
  [("x", -1., 1.), ("y", -1., 1.), ("z", -1., 1.), ("w", -1., 1.)],
  poly3)

## Sliceplorer definition stuff
spec = poly3def
using Sobol

function func1d(f, pt, d)
  function(x)
    xx = copy(pt)
    xx[d] = x
    f(xx)
  end
end

function func_samples(f, rng, pt, d)
  f1d = func1d(f, pt, d)
  samples = range(rng[2], stop=rng[3], length=51) # 50 samples for now
  [(x, f1d(x)) for x = samples]
end

function sample_dim(f, rng, fps, d)
  [func_samples(f, rng, fp, d) for fp = fps]
end

function sliceplorer_samps(spec, n=500, ranges=spec.ranges)
  # Create focus points for slicing
  dmins = [x[2] for x in ranges]
  dmaxs = [x[3] for x in ranges]
  seq = SobolSeq(dmins, dmaxs)
  focuspoints = hcat([next!(seq) for i = 1:n])

  [sample_dim(spec.f, ranges[d], focuspoints, d) for d=eachindex(ranges)]
end

# Now we visualize the result
using Compose
using ColorTypes
using Plots

function plot_dim(dimname, curve_samps, y_rng)
  xs = [v[1] for v = curve_samps[1]]
  ys = [[vv[2] for vv = v] for v = curve_samps]
  plot(xs, hcat(ys...),
       #linealpha=0.6, linecolor=:black, linewidth=0.5,
       linealpha=1.0, linecolor=:black, linewidth=0.5,
       legend=false,
       xlabel=dimname)
       #ylims=y_rng, xlabel=dimname)
end

function sliceplorer_plot(spec, n=500, ranges=spec.ranges)
  samps = sliceplorer_samps(spec, n, ranges)

  # figure out the min/max y values for the scale
  y_vals = [x[2] for x = hcat(hcat(values(samps)...)...)]
  y_rng = (minimum(y_vals), maximum(y_vals))

  plots = [plot_dim(ranges[d][1], samps[d], y_rng) for d = keys(samps)]
  plot(plots..., layout=(length(plots), 1))
end

## The GUI stuff

using Interact

function plot_ui(spec)
  # The main plot
  plt = Observable(sliceplorer_plot(spec, 500))

  # Sliders for each dimension
  sliders = [dim_slider(r[1], r[2], r[3]) for r = spec.ranges]
  # turn sliders into an observable dict
  oo = map((x...) -> collect(x), sliders[1], sliders[2:length(sliders)]...)
  sv = map((vs, ns) -> [(x[1],x[2][1],x[2][2]) for x = zip(ns,vs)], oo, [x[1] for x=spec.ranges])
  #sv = map((v, n) -> Dict(zip(n,v)), oo, [x[1] for x=spec.ranges])

  # Slider for slice count
  ns = slider(1:1000, 500)
  nslider = Widget{:nslider}(OrderedDict(:label => "Slices", :widget => ns), output=observe(ns))
  @layout! nslider hbox(:label, :widget, observe(_))

  # Button to trigger the updates
  update = button("Update")

  # Update the plot when the button is clicked
  map!(t -> sliceplorer_plot(spec, nslider[], sv[]), plt, update)

  # The final layout of all the widgets
  vbox(plt, Interact.hline(), hbox(sliders...), nslider, update)
end

function dim_slider(dn, xmn, xmx)
  #d = OrderedDict(:label => dn, :slider => rangepicker(range(xmn, stop=xmx, length=101), readout=false))
  d = OrderedDict(:label => dn, :lb => widget(-1.), :ub => widget(1.))
  o = Observable((-1., 1.))
  map!((mn,mx) -> (mn, mx), o, d[:lb], d[:ub])
  w = Widget{:dimslider}(d, output=observe(o))
  @layout! w hbox(:label, :lb, :ub)
end

max_stars_repo_name: gabysbrain/multid-slicing

max_issues_repo_name: gabysbrain/multid-slicing

max_forks_repo_name: gabysbrain/multid-slicing


text:
"""
    planar push block
        particle with contacts at each corner
"""
struct PlanarPush{T} <: Model{T}
	# dimensions
	nq::Int # generalized coordinates
	nu::Int # controls
	nw::Int # parameters
	nc::Int # contact points

    mass_block::T
	mass_pusher::T

    inertia::T
    μ_surface::Vector{T}
	μ_pusher::T
    gravity::T

    contact_corner_offset::Vector{Vector{T}}
end

# Kinematics
r_dim = 0.1

function sd_2d_box(p, pose)
	x, y, θ = pose
	Δ = rotation_matrix(-θ) * (p - pose[1:2])
	s = 10
	sum(Δ.^s)^(1/s) - r_dim
end

# contact corner
cc1 = [r_dim, r_dim]
cc2 = [-r_dim, r_dim]
cc3 = [r_dim, -r_dim]
cc4 = [-r_dim, -r_dim]

contact_corner_offset = [cc1, cc2, cc3, cc4]

# Parameters
μ_surface = 0.5  # coefficient of friction
μ_pusher = 0.5
gravity = 9.81
mass_block = 1.0   # mass
mass_pusher = 10.0
inertia = 1.0 / 12.0 * mass_block * ((2.0 * r_dim)^2 + (2.0 * r_dim)^2)

# rnd = 0.01
# dim = [r_dim, r_dim]
# dim_rnd = [r_dim - rnd, r_dim - rnd]

# Methods
M_func(model::PlanarPush, q) = Diagonal([model.mass_block, model.mass_block,
	model.inertia, model.mass_pusher, model.mass_pusher])

function C_func(model::PlanarPush, q, q̇)
	[0.0, 0.0, 0.0, 0.0, 0.0]
end

function rotation_matrix(x)
	[cos(x) -sin(x); sin(x) cos(x)]
end

function ϕ_func(model::PlanarPush, q)
    p_block = q[1:3]
	p_pusher = q[4:5]

	sdf = sd_2d_box(p_pusher, p_block)

    [sdf]
end

function B_func(model::PlanarPush, q)
	[0.0 0.0;
	 0.0 0.0;
	 0.0 0.0;
	 1.0 0.0;
	 0.0 1.0]
end

function N_func(model::PlanarPush, q)
	ϕ = ϕ_func(model, q)
	vec(Symbolics.jacobian(ϕ, q))
end

function p_func(model, x)
    pos = x[1:2]
    θ = x[3]
    R = rotation_matrix(θ)

    [(pos + R * model.contact_corner_offset[1])[1:2];
     (pos + R * model.contact_corner_offset[2])[1:2];
     (pos + R * model.contact_corner_offset[3])[1:2];
     (pos + R * model.contact_corner_offset[4])[1:2]]
end

function P_func(model::PlanarPush, q)
	pf = p_func(model, q)
	P_block = Symbolics.jacobian(pf, q)

	p_block = q[1:3]
	p_pusher = q[4:5]

	ϕ = ϕ_func(model, q)
	N = vec(Symbolics.jacobian(ϕ, q))

	N_pusher = N[4:5]

	n_dir = N_pusher ./ sqrt(N_pusher[1]^2.0 + N_pusher[2]^2.0)
	t_dir = [-n_dir[2]; n_dir[1]]

	r = p_pusher - p_block[1:2]
	m = r[1] * t_dir[2] - r[2] * t_dir[1]

	P = [t_dir[1]; t_dir[2]; m; -t_dir[1]; -t_dir[2]]

	return [P_block; transpose(P)]
end

function residual(model, z, θ, κ)
    nq = model.nq
    nu = model.nu
    nc = model.nc
    nc_impact = 1

    nb = 3 * 4 + 2 * 1

    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    h = θ[2nq + nu .+ (1:1)]

    q2 = z[1:nq]
    γ1 = z[nq .+ (1:nc_impact)]
	s1 = z[nq + nc_impact .+ (1:nc_impact)]

	ψ1 = z[nq + 2 * nc_impact .+ (1:5)]
	b1 = z[nq + 2 * nc_impact + 5 .+ (1:9)]
	sψ1 = z[nq + 2 * nc_impact + 5 + 9 .+ (1:5)]
	sb1 = z[nq + 2 * nc_impact + 5 + 9 + 5 .+ (1:9)]

	ϕ = ϕ_func(model, q2)
	N = vec(Symbolics.jacobian(ϕ, q2))

	# λ1 = [b1; γ1]
	P = P_func(model, q2)
    vT = P * (q2 - q1) ./ h[1]

	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(a -> M_func(model, a), (a, b) -> C_func(model, a, b), qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(a -> M_func(model, a), (a, b) -> C_func(model, a, b), qm2, vm2)

    d = (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2#
            + B_func(model, qm2) * u1
            + N * γ1[1]
            + transpose(P) * b1)

    [
	 d;

	 s1 .- ϕ;

	 ψ1[1] .- model.μ_surface[1] * model.mass_block * model.gravity * h[1] * 0.25;

	 ψ1[2] .- model.μ_surface[2] * model.mass_block * model.gravity * h[1] * 0.25;

	 ψ1[3] .- model.μ_surface[3] * model.mass_block * model.gravity * h[1] * 0.25;

	 ψ1[4] .- model.μ_surface[4] * model.mass_block * model.gravity * h[1] * 0.25;

	 ψ1[5] .- model.μ_pusher * γ1[1];

	 vT - sb1;

	 γ1 .* s1 .- κ[1];
	 cone_product([ψ1[1]; b1[1:2]], [sψ1[1]; sb1[1:2]]) - [κ[1]; 0.0; 0.0];
	 cone_product([ψ1[2]; b1[2 .+ (1:2)]], [sψ1[2]; sb1[2 .+ (1:2)]]) - [κ[1]; 0.0; 0.0];
	 cone_product([ψ1[3]; b1[4 .+ (1:2)]], [sψ1[3]; sb1[4 .+ (1:2)]]) - [κ[1]; 0.0; 0.0];
	 cone_product([ψ1[4]; b1[6 .+ (1:2)]], [sψ1[4]; sb1[6 .+ (1:2)]]) - [κ[1]; 0.0; 0.0];
	 cone_product([ψ1[5]; b1[8 .+ (1:1)]], [sψ1[5]; sb1[8 .+ (1:1)]]) - [κ[1]; 0.0];
    ]
end

# Dimensions
nq = 5 # configuration dimension
nu = 2 # control dimension
nc = 5 # number of contact points
nc_impact = 1
nf = 3 # number of faces for friction cone pyramid
nb = (nc - nc_impact) * nf + (nf - 1) * nc_impact

planarpush = PlanarPush(nq, nu, 0, nc,
			mass_block, mass_pusher,
			inertia, [μ_surface for i = 1:nc], μ_pusher, gravity,
			contact_corner_offset)


text:
# abstract type Model end

struct Dimensions
    q::Int         # configuration
    u::Int         # control
	w::Int         # disturbance
	c::Int         # contact points
    quat::Int      # number of quaternions in the configuration
end

# https://github.com/HarvardAgileRoboticsLab/drake/blob/75b260c9eb250d08ffbbf3daa80758e4fe558d7f/drake/matlab/solvers/trajectoryOptimization/VariationalTrajectoryOptimization.m
function lagrangian_derivatives(model::Model, q, v)
	D1L = -1.0 * C_fast(model, q, v)
    D2L = M_fast(model, q) * v
	return D1L, D2L
end

function dynamics(model::Model, h, q0, q1, u1, w1, Λ1, q2)
	# evalutate at midpoint
	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	# return 0.0
	return (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2
		+ transpose(B_fast(model, qm2)) * u1
		+ transpose(A_fast(model, qm2)) * w1
		# + transpose(J_fast(model, q2)) * λ1
		+ Λ1
		- h[1] * model.joint_friction .* vm2)

end

function dynamics(model::Model, env::Environment, h, q0, q1, u1, w1, λ1, q2)
	Λ1 = transpose(J_func(model, env, q2)) * λ1 #@@@@ maybe need to use J_fast
	dynamics(model, h, q0, q1, u1, w1, Λ1, q2)
end

mutable struct BaseMethods
	L::Any
	M::Any
	B::Any
	A::Any
	# J::Any
	C::Any
	k::Any
end

function BaseMethods()
	function f()
		error("Not Implemented: use instantiate_base!")
		return nothing
	end
	return BaseMethods(fill(f, 6)...)
end

mutable struct DynamicsMethods
	d::Any
	dy::Any
	dθ::Any
	dq0::Any
	dq1::Any
	du1::Any
	dw1::Any
	dγ1::Any
	db1::Any
	dΛ1::Any
	∂q2::Any
end

function DynamicsMethods()
	function f()
		error("Not Implemented: use instantiate_dynamics!")
		return nothing
	end
	return DynamicsMethods(fill(f, 11)...)
end

"""
	get_model(name::String, surf::String)
	Helper function that provides a model where fast functions have been instantiated.
"""
function get_model(name::String;
	model_name::String = name, approx = false)

	#TODO: assert model exists
	path = joinpath(@__DIR__, name)
	model = eval(Symbol(model_name * (surf != "flat" ? "_" * surf : "")))
	instantiate_base!(model, joinpath(path, dynamics, "base.jld2"))
	instantiate_dynamics!(model, joinpath(path, dynamics, "dynamics.jld2"),
		derivs = approx)

	return model
end

function get_gait(name::String, gait::String)
	#TODO: assert model exists
	path = joinpath(@__DIR__, name)
	gait_path = joinpath(path, "gaits/" * gait * ".jld2")

	res = JLD2.jldopen(gait_path) # z̄ x̄ ū h̄ q u γ b

	return res["q"], res["u"], res["γ"], res["b"], mean(res["h̄"])
end

function model_name(model::Model)
    name = Symbol(string(typeof(model).name)[10:end-1])
    return name
end


text:
epochs = ProgressBar(iter:n_iter_max);
loss_epoch = zeros(Float32, n_exp);
grad_norm = zeros(Float32, n_exp_train);
for epoch in epochs
    global p
    if conf["epoch_size"] != -1
        epoch_samples = sample(1:n_exp_train, conf["epoch_size"], replace=false)
    else
        epoch_samples = sample(1:n_exp_train, n_exp_train, replace=false)
    end
    for i_exp in epoch_samples
        batch = rand(batch_size:ntotal)
        grad = Zygote.gradient(x -> loss_neuralode(x, i_exp, batch), p)[1]
        grad_norm[i_exp] = norm(grad, 2)
        if grad_norm[i_exp] > grad_max
            grad = grad ./ grad_norm[i_exp] .* grad_max
        end
        update!(opt, p, grad)
    end
    for i_exp in 1:n_exp
        loss_epoch[i_exp] = loss_neuralode(p, i_exp)
    end
    loss_train = mean(loss_epoch[1:n_exp_train]);
    loss_val = mean(loss_epoch[n_exp_train+1 : n_exp_train+n_exp_val]);
    loss_test = mean(loss_epoch[n_exp_train+n_exp_val+1:end]);
    g_norm = mean(grad_norm)

    set_description(epochs, string(@sprintf("Loss train %.4e tol %d/%d lr %.1e",
                            loss_train, no_change, n_iter_tol, opt[1].eta)))
    cb(p, loss_train, loss_val, loss_test, g_norm);

    if n_iter_tol > 0
        if checkconvergence()
            break
        end
    end
end

if !parsed_args["disable-display"]
    for i_exp in n_exp_train:n_exp
        @suppress_err cbi(p, i_exp)
    end
end


text:
function summary(io::IO, x::fdawarp)
    @printf("-------------------------------------\n")
    @printf("Number of Functions: %d\n", size(x.f,2))
    @printf("Number of time points: %d\n", size(x.f,1))
    @printf("-------------------------------------\n")
    @printf("Alignment with lambda: %f\n", x.lambda)
    @printf("Method: %s\n", x.method)
    @printf("Optimization Method: %s\n", x.omethod)
    @printf("Original Variance: %f\n", x.orig_var)
    @printf("Amplitude Variance: %f\n", x.amp_var)
    @printf("Phase Variance: %f\n", x.phase_var)
    @printf("Number of Iterations: %d\n", length(x.cost))
end


text:
"""
    SVM

A super type for both classifiers and regressors that are implemented as Support Vector Machines.
"""
abstract type SVM end

"""
    LSSVC <: SVM
    LSSVC(; kernel=:rbf, γ=1.0, σ=1.0, degree=0)

The type to hold a Least Squares Support Vector Classifier.

# Fields
- `kernel::Symbol`: The kind of kernel to use for the non-linear mapping of the data. Can be one of the following: `:rbf`, `:linear`, or `:poly`.
- `γ::Float64`: The gamma hyperparameter that is intrinsic of the Least Squares version of the Support Vector Machines.
- `σ::Float64`: The hyperparameter for the RBF kernel.
- `degree::Int`: The degree of the polynomial kernel. Only used if `kernel` is `:poly`.

"""
mutable struct LSSVC <: SVM
    kernel::Symbol
    γ::Float64
    σ::Float64
    degree::Int
end

LSSVC(; kernel=:rbf, γ=1.0, σ=1.0, degree::Int=0) = LSSVC(kernel, γ, σ, degree)

"""
    LSSVR <: SVM
    LSSVR(; kernel=:rbf, γ=1.0, σ=1.0, degree=0)

The type to hold a Least Squares Support Vector Regressor.

# Fields
- `kernel::Symbol`: The kind of kernel to use for the non-linear mapping of the data. Can be one of the following: `:rbf`, `:linear`, or `:poly`.
- `γ::Float64`: The gamma hyperparameter that is intrinsic of the Least Squares version of the Support Vector Machines.
- `σ::Float64`: The hyperparameter for the RBF kernel.
- `degree::Int`: The degree of the polynomial kernel. Only used if `kernel` is `:poly`.

"""
mutable struct LSSVR <: SVM
    kernel::Symbol
    γ::Float64
    σ::Float64
    degree::Int
end

LSSVR(; kernel=:rbf, γ=1.0, σ=1.0, degree::Int=0) = LSSVR(kernel, γ, σ, degree)

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
_brightness_contrast(color, matrix::AbstractMatrix{T}, brightness, contrast) where {T} =
    @. color(matrix / 255 * T(contrast) + T(brightness) / 255)

"""
    brightness_contrast(image::AstroImage; brightness_range = 0:255, contrast_range = 1:1000, header_number = 1)

Visualize the fits image by changing the brightness and contrast of image.

Users can also provide their own range as keyword arguments.
"""
function brightness_contrast(img::AstroImage{T,C,N}; brightness_range = 0:255,
                             contrast_range = 1:1000, header_number = 1) where {T,C,N}
    @manipulate for brightness  in brightness_range, contrast in contrast_range
        _brightness_contrast(C, img.data[header_number], brightness, contrast)
    end
end

# This is used in Jupyter notebooks
Base.show(io::IO, mime::MIME"text/html", img::AstroImage; kwargs...) =
    show(io, mime, brightness_contrast(img), kwargs...)


text:
#=
Contains functions/type prototypes for the mesh.
=#

#= Mesh =#
#Subtypes should contain all variables that define a mesh.
#Should make an effort to follow the naming conventions established
#in existing methods in order to give compatibility between methods
abstract type AbstractMesh end

mutable struct Mesh{T} <: AbstractMesh
    N::Int

    δ_min::Vector{T}
    digits::Vector{Union{Int, Nothing}}

    l::Int
    δ::Vector{T}
    δ⁰::Vector{T}
    Δ::Vector{T}
    ρ::Vector{T}
    ρ_min::T
    a::Vector{Int}
    b::Vector{Int}
    a⁰::Vector{Int}
    b⁰::Vector{Int}

    is_anisotropic::Bool
    aₜ::T
    only_granular::Bool

    # Override constructor for different default meshes for
    # different poll techniques.
    Mesh(N::Int64) = Mesh{Float64}(N)
    function Mesh{T}(N::Int64) where T
        mesh = new()

        mesh.N = N

        mesh.aₜ = 0.1

        mesh.l = 0
        mesh.a⁰ = Vector{Int}(undef, N)
        mesh.b⁰ = Vector{Int}(undef, N)
        mesh.a = Vector{Int}(undef, N)
        mesh.b = Vector{Int}(undef, N)
        mesh.δ = Vector{T}(undef, N)
        mesh.Δ = Vector{T}(undef, N)
        return mesh
    end
end

"""
    MeshSetup!(p::DSProblem)

Sets up the Mesh with the parameters defined for problem.
"""
MeshSetup!(p::DSProblem) = MeshSetup!(p, p.config.mesh)
function MeshSetup!(p::DSProblem, m::Mesh)
    m.is_anisotropic = p.config.poll isa OrthoMADS

    m.δ_min = p.granularity
    m.digits = map(get_decimal_places, p.granularity)
    m.only_granular = all(p.granularity .> 0)

    init_a_and_b!(p, m)

    SetMeshParameters!(m)
    m.δ⁰ = copy(m.δ)
end

"""
    MeshUpdate!(mesh::Mesh, ::AbstractPoll, result::IterationOutcome, dir::Union{Vector,Nothing})

Implements update rule from Audet, Le Digabel & Tribes 2019 adapted for progressive
barrier constraints with Audet & Dennis 2009 expression 2.4.

`dir` is the direction of success of the iteration, equal to `nothing`, if there is none.
"""
function MeshUpdate!(m::Mesh, ::AbstractPoll, result::IterationOutcome, dir::Union{Vector,Nothing})
    if result == Unsuccessful
        for i=1:m.N
            decrease_a_and_b!(m, i)
        end
        m.l -= 1
        SetMeshParameters!(m)
    elseif result == Dominating
        for i=1:m.N
            increase_a_and_b!(m, i, dir)
        end
        m.l += 1
        SetMeshParameters!(m)
    end
end

"""
    SetMeshParameters!(m::Mesh)

Sets the mesh parameters: mesh size vector δ, poll size vector Δ, and ratio vector ρ.
"""
function SetMeshParameters!(m::Mesh)
    SetMeshSizeVector!(m)
    SetPollSizeVector!(m)
    SetRatioVector!(m)
end

"""
    SetMeshSizeVector!(m::Mesh)

Sets the mesh size vector δ.
"""
function SetMeshSizeVector!(m::Mesh)
    for i = 1:m.N
        m.δ[i] = 10.0^(m.b[i] - abs(m.b[i] - m.b⁰[i]))

        if m.δ_min[i] > 0
            m.δ[i] = m.δ_min[i] * max.(1, m.δ[i])
        end
    end
end

"""
    SetPollSizeVector!(m::Mesh)

Sets the poll size vector Δ.
"""
function SetPollSizeVector!(m::Mesh)
    for i = 1:m.N
        m.Δ[i] = m.a[i] * (10.0 ^ m.b[i])

        if m.δ_min[i] > 0
            m.Δ[i] = m.δ_min[i] * m.Δ[i]
        end
    end
end

"""
    SetRatioVector!(m::Mesh)

Sets the ratio vector ρ.
"""
function SetRatioVector!(m::Mesh{T}) where T
    m.ρ = m.Δ ./ m.δ

    if any(m.δ_min .== 0)
        m.ρ_min = minimum(m.ρ[m.δ_min .== 0])
    else
        m.ρ_min = typemin(T)
    end
end

"""
    init_a_and_b!(p::DSProblem, m::Mesh)

Initializes the parameters `a` and `b`, which define the mesh and poll sizes.
"""
function init_a_and_b!(p::DSProblem, m::Mesh)
    for i = 1:p.N
        α = get_poll_size_estimate(p.user_initial_point[i], p.lower_bounds[i], p.upper_bounds[i])

        if m.δ_min[i] > 0
            α_scaled = α/m.δ_min[i]
        else
            α_scaled = α
        end

        b⁰ = floor(Int, log10(α_scaled))

        # b cannot has to be nonnegative for granular indexes
        if b⁰ < 0 && m.δ_min[i] > 0
            b⁰ = 0
        end

        # the NOMAD way
        # if b⁰ < 0
        #     b⁰ = 0
        # end

        a⁰ = α_scaled / (10.0^b⁰)

        if a⁰ < 1.5
            a⁰ = 1
        elseif a⁰ < 3.5
            a⁰ = 2
        else
            a⁰ = 5
        end

        m.a⁰[i] = a⁰
        m.b⁰[i] = b⁰

        m.a[i] = a⁰
        m.b[i] = b⁰
    end
end

"""
    get_poll_size_estimate(x⁰::T, lower_bound::Union{T,Nothing}, upper_bound::Union{T,Nothing})::T where T

Calculates the initial poll size, as given in Audet, Le Digabel & Tribes 2019 expression 3.3,
using the initial point and variable bounds, if specified.
"""
function get_poll_size_estimate(x⁰::T, lower_bound::Union{T,Nothing}, upper_bound::Union{T,Nothing})::T where T
    if !(lower_bound === nothing) && !(upper_bound === nothing) && (lower_bound < upper_bound)
        (upper_bound - lower_bound)/10
    else
        if !(lower_bound === nothing) && (upper_bound === nothing)
            finite_bound = lower_bound
        elseif !(upper_bound === nothing) && (lower_bound === nothing)
            finite_bound = upper_bound
        else
            finite_bound = nothing
        end

        if !(finite_bound === nothing) && !(finite_bound === x⁰)
            abs(x⁰ - finite_bound)/10
        elseif !(x⁰ == 0)
            abs(x⁰)/10
        else
            1
        end
    end
end

function increase_a_and_b!(m::Mesh, i::Int, dir::Union{Vector,Nothing})
    if dir === nothing || !m.is_anisotropic || (abs(dir[i])/m.ρ[i] > m.aₜ) || (m.δ_min[i] == 0  && m.δ[i] < m.δ⁰[i] && m.ρ[i] > m.ρ_min^2)
        if m.a[i] === 1
            m.a[i] = 2
        elseif m.a[i] === 2
            m.a[i] = 5
        else
            m.a[i] = 1
            m.b[i] += 1
        end
    end
end

function decrease_a_and_b!(m::Mesh, i::Int)
    if m.a[i] === 1
        if m.b[i] > 0 || m.δ_min[i] == 0
            m.a[i] = 5
            m.b[i] -= 1
        end
    elseif m.a[i] === 2
        m.a[i] = 1
    else
        m.a[i] = 2
    end
end

function get_decimal_places(granularity)::Union{Int, Nothing}
    granularity == 0 && return nothing

    str_split = split(string(granularity), ".")

    if length(str_split) == 1 || (length(str_split[2]) == 1 && str_split[2][1] == '0')
        return 0
    else
        return length(str_split[2])
    end
end


text:
@testset "CoreModel generic interface" begin
    model = load_model(CoreModel, model_paths["e_coli_core.mat"])

    @test reaction_stoichiometry(model, "EX_ac_e") == Dict("ac_e" => -1)
    @test reaction_stoichiometry(model, 44) == Dict("ac_e" => -1)
end


text:
using Revise
using LibPQ
using DataKnots: DataKnot
using FunSQL: FunSQL, SQLTable, From, Select, Where, Join, Group, Agg, Fun, Get, to_sql, normalize
using PostgresCatalog
using Tables

FunSQL.SQLTable(tbl::PostgresCatalog.PGTable) =
    SQLTable(Symbol(tbl.schema.name), Symbol(tbl.name), Symbol[Symbol(col.name) for col in tbl])

conn = LibPQ.Connection("postgresql:///synpuf-hcfu?host=/var/run/postgresql")

function run(q)
    sql = to_sql(normalize(q))
    println(sql, ";")
    println()
    result = execute(conn, sql)
    println(convert(DataKnot, result))
    println('-'^80)
end

cat = PostgresCatalog.introspect(conn)
for tbl in cat["public"]
    @eval $(Symbol(tbl.name)) = SQLTable($tbl)
end

q = From(person)
run(q)

q = person |>
    Select(Get.person_id, Get.year_of_birth)
run(q)

q = person |>
    Group(Get.gender_concept_id)
run(q)

q = person |>
    Group(Get.gender_concept_id) |>
    Select(
        Get.gender_concept_id,
        Agg.Count())
run(q)

q = person |>
    Group() |>
    Select(Agg.Count())
run(q)

q = person |>
    Group() |>
    Select(
        Agg.Count(),
        Agg.Min(Get.year_of_birth),
        Agg.Max(Get.year_of_birth))
run(q)

q = (:p_group => person |> Group()) |>
    Select(
        Agg.Count(over=Get.p_group),
        Agg.Min(Get.year_of_birth, over=Get.p_group),
        Agg.Max(Get.year_of_birth, over=Get.p_group))
run(q)

q = concept |>
    Join(concept_ancestor,
         Fun."="(Get.concept_id, Get.descendant_concept_id)) |>
    Select(Get.concept_id,
           Get.descendant_concept_id,
           Get.ancestor_concept_id)
run(q)

q = concept |>
    Join(concept_ancestor,
         Fun."="(Get.concept_id, Get.descendant_concept_id)) |>
    Join(:ancestor => concept |> Where(Fun."="(Get.concept_id, "312327")),
         Fun."="(Get.ancestor_concept_id, Get.ancestor.concept_id)) |>
    Select(Get.concept_id,
           Get.descendant_concept_id,
           Get.ancestor_concept_id,
           Get.ancestor.concept_id)
run(q)

q = person |>
    Join(:condition => condition_occurrence |> Group(Get.person_id),
         Fun."="(Get.person_id, Get.condition.person_id)) |>
    Select(
        Get.person_id,
        Agg.Count(),
        Agg.Min(Get.condition_start_date),
        Agg.Max(Get.condition_start_date))
run(q)

q = person |>
    Join(:condition => condition_occurrence |> Group(Get.person_id),
         Fun."="(Get.person_id, Get.condition.person_id)) |>
    Where(Fun.">"(Agg.Count(), 2)) |>
    Select(Get.person_id, Agg.Count())
run(q)

q = person |>
    Join(:condition => condition_occurrence |> Group(Get.person_id),
         Fun."="(Get.person_id, Get.condition.person_id)) |>
    Group("# conditions" => Agg.Count()) |>
    Select(
        Get."# conditions",
        "# patients" => Agg.Count(),
        "min year of birth" => Agg.Min(Get.year_of_birth))
run(q)

q = person |>
    Join(:condition => condition_occurrence |> Group(Get.person_id),
         Fun."="(Get.person_id, Get.condition.person_id)) |>
    Join(:visit => visit_occurrence |> Group(Get.person_id),
         Fun."="(Get.person_id, Get.visit.person_id)) |>
    Select(
        Get.person_id,
        Agg.Count(over=Get.condition),
        Agg.Min(Get.condition_start_date, over=Get.condition),
        Agg.Max(Get.condition_start_date, over=Get.condition),
        Agg.Count(over=Get.visit),
        Agg.Min(Get.visit_start_date, over=Get.visit),
        Agg.Max(Get.visit_start_date, over=Get.visit))
run(q)


text:
using CGE

edges, weights, vweights, comm, clusters, embed, verbose, land, forced, method, directed, split, seed, samples = parseargs()
distances = zeros(length(vweights))
init_edges = Array{Int,2}(undef,0,0)
init_vweights = Vector{Float64}()
init_eweights = Vector{Float64}()
init_embed = Array{Float64,2}(undef,0,0)
v_to_l = Int[]
if land != -1
    init_edges = copy(edges)
    init_vweights = copy(vweights)
    init_eweights = copy(weights)
    init_embed = copy(embed)
    distances, embed, comm, edges, weights, vweights, v_to_l = landmarks(edges,
            weights, vweights, clusters, comm, embed, verbose, land, forced, method, directed)
end
if directed
    results = wGCL_directed(edges, weights, comm, embed, distances, vweights, init_vweights,
                v_to_l, init_edges, init_eweights, init_embed, split, seed, samples, verbose)
else
    results = wGCL(edges, weights, comm, embed, distances, vweights, init_vweights,
                v_to_l, init_edges, init_eweights, init_embed, split, seed, samples, verbose)
end
println(results)


text:
"""
"""
function fix_type
end

fix_type(x::Any, T) = convert(T, fix_type(x))

fix_type(x::Any) = x

function fix_type(x::AbstractArray)::Array
    result = reshape(
        [[fix_type(element) for element in x]...],
        size(x),
        )
    return result
end

function fix_type(
        x::AbstractDict;
        default_key_type::Type=Any,
        default_value_type::Type=Any,
        )::Dict
    if length(x) == 0
        result = Dict{default_key_type, default_value_type}()
    else
        keys_eltype=eltype( fix_type( collect( keys(x) ) ) )
        values_eltype=eltype( fix_type( collect( values(x) ) ) )
        result = Dict{keys_eltype, values_eltype}()
        for k in keys(x)
            result[k] = x[k]
        end
    end
    return result
end