text:
module arg_handler

const doc = """This is the Julia vfld program
I can help you with:
Convert VFLD to SQLite (cmd: vfld_to_sqlite)
Rename files from date+leadtime to date only (cmd: date_naming)

date_naming: Before use, please ensure if you want to distunguish between analysis and forecasts. If you want that, then date_naming is not for you.

Usage:
  vfld.jl [vfld_to_sqlite] [date_naming] [options]
  vfld.jl -h | --help

Options:
  -h --help                     Show this screen.
  --starttime=<YYYY-MM-DD-HH>   start time.
  --endtime=<YYYY-MM-DD-HH>     end time.
  --file-prefix=<str>           Prefix of vfld files (eg. "vfldER5") [default: "vfld"]
  --file-postfix=<str>          Postfix of vfld files (eg. "01") [default: ""]
  --indir=<str>                 Input directory [default: ~/].

vfld_to_sqlite:
  --sqlite-file=<str>           Output to files [default: ~/out.db]

date_naming:
  --max-leadtime=<int>          Only use up x hours leadtime (eg. 2020-01-01-00+03 will be 2020-01-01-03, only if --max-leadtime is above 2) [default: 2]
  --leadtime-characters=<int>   How many of the last digits is the leadtime info? [default: 2]
  --outdir=<str>                Which directory to write renamed files to [default: .]

Development shortcuts:
    julia --project=~/git/vfld/ vfld_util.jl vfld_to_sqlite --starttime=2013-10-01-00 --endtime=2013-10-02-00 --file-prefix=vfld --sqlite-file=/home/kah/git/vfld/scr/test.db --indir=/home/kah/tmp/vfld/ERA5/
    julia --project=~/git/vfld/ vfld_util.jl date_naming --starttime=2013-10-01-00 --endtime=2020-11-01-00 --file-prefix=vfldDKREA --leadtime-characters=2 --max-leadtime=2 --indir=/home/kah/tmp/vfld/DKREA/201310/ --outdir=/home/kah/tmp/vfld/DKREA_RENAMED/

Notes:
    Assumes the naming of vfld files follows a structure like so:
    vfldER5201310312300 (must end with 10 or 12 date characters)

    It is recommended not to cover more than one month with "vfld_to_sqlite" command.
"""

const avail_commands = ["vfld_to_sqlite", "date_naming"]

function main_args(args)
    help_msg(args)
    cmds = make_dict(args)
    cmd_message = command_msg(cmds)
    return cmd_message
end

function help_msg(args)
    """Checks if just the doc string should be printed"""
    if "-h" in args || "--help" in args
        print(doc)
        exit(0)
    end
end

function make_dict(args)
    """Makes a dictionary of input arguments"""
    cmds = Dict()

    for k in args
        if startswith(k,"--")       # Identifier for an option
            opt = split(k,"=")
            cmds[opt[1]] = opt[2]
        elseif startswith(k,"-")    # Identifier for a flag
            cmds[k] = "flag"
        else                        # Identifier for a command
            cmds[k] = "command"
        end
    end
    return cmds
end

function command_msg(cmds)
    """Checks which command we are supposed to run"""

    for cmd in avail_commands

        # command: vfld_to_sqlite
        if cmd == "vfld_to_sqlite" && haskey(cmds,cmd) # Determine which command to run
            cmd_message = check_options(cmds,cmd)
            return cmd_message # We only allow one command at a time for now
        elseif cmd == "date_naming" && haskey(cmds,cmd)
            cmd_message = check_options(cmds,cmd)
            return cmd_message
        end

    end
end

function check_options(cmds, cmd)
     """Further check for necessary options for this specific command"""
    if cmd == "vfld_to_sqlite"

        starttime = key_check("--starttime", cmds)
        endtime = key_check("--endtime", cmds)

        if starttime === missing || endtime === missing
            raise_missing_option("--starttime and --endtime")
        end

        file_prefix = key_check("--file-prefix", cmds)
        file_prefix===missing ? file_prefix = "vfld" : nothing

        file_postfix = key_check("--file-postfix", cmds)
        file_postfix===missing ? file_postfix = "" : nothing

        indir = key_check("--indir", cmds)
        indir===missing ? indir = "~/" : nothing

        sqlitefile = key_check("--sqlite-file", cmds)
        sqlitefile===missing ? sqlitefile = "~/out.db" : nothing

        cmd_message = String.((cmd, starttime, endtime, file_prefix, file_postfix, indir, sqlitefile))

        return cmd_message
    elseif  cmd == "date_naming"

        starttime = key_check("--starttime", cmds)
        starttime===missing ? starttime = "1970-01-01-00" : nothing

        endtime = key_check("--endtime", cmds)
        endtime===missing ? endtime = "2100-01-01-00" : nothing

        file_prefix = key_check("--file-prefix", cmds)
        file_prefix===missing ? file_prefix = "vfld" : nothing

        file_postfix = key_check("--file-postfix", cmds)
        file_postfix===missing ? file_postfix = "" : nothing

        indir = key_check("--indir", cmds)
        indir===missing ? indir = "~/" : nothing

        maxleadtime = key_check("--max-leadtime", cmds)
        maxleadtime===missing ? maxleadtime = 2 : nothing

        leadtimechars = key_check("--leadtime-characters", cmds)
        leadtimechars===missing ? leadtimechars = 2 : nothing

        outputdir = key_check("--outdir", cmds)
        outputdir===missing ? outputdir = "." : nothing

        cmd_message = String.((cmd, starttime, endtime, file_prefix, file_postfix, indir, maxleadtime, leadtimechars, outputdir))

        return cmd_message
    end

end

function key_check(key,dict)
    """Checks if a key exist in a dict and if so, return the value"""
    if haskey(dict,key)
        value = dict[key]
        return value
    else
        return missing
    end
end

function raise_missing_option(option)
    println(option*" necessary for this command")
    exit(0)
end

end


text:
module G
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    [1700    0    0;
0 1700    0;
0    0 1700]
end
end

ext: jl


text:
"""
    fitexp(A::NMRData{T,2}, selector)

# Examples

```julia
y = loadnmr("exampledata/pseudo2D_XSTE/test.ft1")
y = settval(y, LinRange(0.05, 0.95, 10))
rate, plt = fitexp(y, Between(7.5,9))
display(plt)
```
"""
fitexp(A::D, selector) where {D<:NMRData{T,2} where T} =
        fitexp(SimpleTraits.trait(HasPseudoDimension{D}), A, selector)

function fitexp(::Type{HasPseudoDimension{D}}, A::D, selector) where {D<:NMRData{T, 2} where T}
    # @show A[selector,:]
    # @show data(A[selector,:])
    integrals = vec(data(sum(A[selector,:],dims=X)))
    integrals /= maximum(integrals)

    t = tval(A)
    tscale = maximum(t) / 2

    @. model(t, p) = p[1] * exp(-p[2] * t)
    p0 = [1.0, 1.0] # take a rough guess of initial rate based on observed time points
    fit = curve_fit(model, t / tscale, integrals, p0)
    pfit = coef(fit)
    err = stderror(fit)
    k = (pfit[2] ± err[2]) / tscale

    x = LinRange(0, maximum(t)*1.1, 50)
    yfit = model(x, [pfit[1], pfit[2] / tscale])

    p1 = plot(x, yfit)
    scatter!(p1, t, integrals, legend=nothing)
    xlabel!(p1, label(A,Ti))
    ylabel!(p1, "Integrated signal")
    ylims!(p1, 0, maximum(yfit)*1.1)
    title!(p1, "")

    p2 = plot(A[:,1],linecolor=:grey)
    p2 = plot!(p2, A[selector,1], fill=(0,:orange), linecolor=:red)
    title!(p2, "")
    #xlims!(p2, (0,10))

    plt = plot(p1, p2, layout=(2,1), size=(600,600))

    return k, plt
end

"""
    fitdiffusion(A::NMRData{T,2}, selector; δ=0.004, Δ=0.1, σ=0.9, Gmax=0.55)

# Examples

```julia
y = loadnmr("exampledata/pseudo2D_XSTE/test.ft1")
y = settval(y, LinRange(0.05, 0.95, 10))
rH, D = fitdiffusion(y, Between(7.5,9))
```
"""
fitdiffusion(spec::X, selector; δ=0.004, Δ=0.1, σ=0.9, Gmax=0.55, solvent=:h2o, T=298, showplot=true) where {X<:NMRData{T,2} where T} =
        fitdiffusion(SimpleTraits.trait(HasPseudoDimension{X}), spec, selector, δ, Δ, σ, Gmax, solvent, T, showplot)

function fitdiffusion(::Type{HasPseudoDimension{X}}, spec::X, selector, δ, Δ, σ, Gmax, solvent, T, showplot) where {X<:NMRData{T, 2} where T}
    g = tval(spec) * Gmax
    γ = 2.675e8;
    q = (γ*δ*σ*g).^2 * (Δ - δ/3)
    lab = label(spec)

    y = settval(spec, q) # create copy of data containing q values
    label!(y,Ti,"q / s m-2")

    D, plt = fitexp(y, selector)

    if solvent==:h2o
        A = 802.25336
        a = 3.4741e-3
        b = -1.7413e-5
        c = 2.7719e-8
        gamma = 1.53026
        T0 = 225.334
    elseif solvent==:d2o
        A = 885.60402
        a = 2.799e-3
        b = -1.6342e-5
        c = 2.9067e-8
        gamma = 1.55255
        T0 = 231.832
    else
        @error "fitdiffusion: solvent not recognised (should be :h2o or :d2o)"
    end

    DT = T - T0
    k = 1.38e-23
    η = A * (DT + a*DT^2 + b*DT^3 + c*DT^4)^(-gamma)
    rH = k*T / (6π*η*0.001 * D) * 1e10 # in Å

    plt.subplots[1].attr[:title] = "$(lab): rH = $(rH) Å"

    showplot && display(plt)

    return rH, D, plt
end


text:
# ------------------------------------------------------------------------------
# --- DESCRIPTION OF MODEL -----------------------------------------------------
# ------------------------------------------------------------------------------
# Discrete space, discrete time model for examining stochastic evolutionary
# processes during range expansion.
#
# Original R Code by Ben Phillips (2015)
# Translated to Julia by Brad Ochocki (2016)
# Adapted for current project by Brad Ochocki (2017)

# ------------------------------------------------------------------------------
# --- PACKAGE MANAGEMENT -------------------------------------------------------
# ------------------------------------------------------------------------------
installPackages = false
updatePackages  = false
loadPackages    = true

# Install required packages ----------------------------------------------------
if installPackages
	Pkg.add("Distributions")
	Pkg.add("DataFrames")
end

# Update packages --------------------------------------------------------------
if updatePackages
	Pkg.update()
end

# Set packages to use ----------------------------------------------------------
if loadPackages
	using Distributions
	using DataFrames
end

# Include utility functions ----------------------------------------------------
include("utilities.jl")

# ------------------------------------------------------------------------------
# --- FUNCTION DEFINITIONS -----------------------------------------------------
# ------------------------------------------------------------------------------

### HELPER FUNCTIONS ###########################################################

# function for naming output files
function outname(p, r)
	return string("correlated_traits_", p, "_" , r, ".csv")
end

### GROWTH FUNCTIONS ###########################################################

# Beverton-Holt per-capits population growth -----------------------------------
function growth(N::Int64, r::Float64, b::Float64)
	# Given [N, r, b] returns expected [per-capita growth rate]
  # N : number of individuals
	# r : low-density per-capita growth rate phenotype
	# b : Beverton-Holt parameter

	# use a logistic link function to map growth phenotype to a growth rate that
	# is bounded between 0 and K.
  if N == 0
    return 0
  else
    return r / (1 + b .* N)
  end
end

### GENETIC VARIANCE FUNCTIONS #################################################

# Initialize population --------------------------------------------------------
function init_inds(n, M, V, ρ, H)
	# Given [n, M, V, C, H] returns [X, D, PD, K, PK, r, Pr], which is
	# referred to in other functions as 'popmatrix'.
  # n   : number of starting individuals
	# M   : vector of mean phenotypes [μD μr]
	# V   : vector of total phenotypic variance [VD Vr]
	# ρ   : vector of correlations [ρA ρE]
	# H   : vector of heritabilities [h²D h²r]

	# X   : spatial location of each individual
	# D   : dispersal genotype
	# PD  : dispersal phenotype (accounts for environmental variance)
	# r   : low-density growth rate (r) genotype
	# Pr  : r phenotype (accounts for environmental variance)
	# ΣA  : additive genetic covariance matrix
	# ΣE  : environmental covariance matrix

	# set all starting locations to zero
	X = zeros(Int64, n)

	ΣA = VA(V, ρ[1], H)
	ΣE = VE(V, ρ[2], H)
  G = rand(MvNormal(M, ΣA), n)'
	P = rand(MvNormal([0.0, 0.0], ΣE), n)' + G

	return DataFrame(X=X, D=G[:, 1], PD=P[:, 1], r=G[:, 2], Pr=P[:, 2])
end

# Calculate non-additive (environmental) phenotypic covariance matrix ----------
function VE(V, ρ, H)
	# Given [V, ρ, H], returns ΣE
	# V  : vector of total phenotypic variance [VD Vr]
  # ρ  : covariance between D and r
	# H  : vector of heritabilities [h²D h²r]

	# ΣE : environmental variance matrix

	return sqrt.((1.0 - H) .* V) * sqrt.((1.0 - H) .* V)' .* [1.0 ρ; ρ 1.0]
end

# Calculate additive genetic covariance matrix ---------------------------------
function VA(V, ρ, H)
	# Given [V, ρ, H], returns ΣA
	# V  : vector of total phenotypic variance [VD Vr]
  # ρ  : covariance between D and r
	# H  : vector of heritabilities [h²D h²r]

	# ΣA : environmental variance matrix

  return sqrt.(H .* V) * sqrt.(H .* V)' .* [1.0 ρ; ρ 1.0]
end

### FUNCTIONS FOR PATCH-RELATED CALCULATIONS ###################################

# Determine which spatial locations to sample from -----------------------------
function accordian_bins(A::Vector{Int64})

	ls = minimum(A)
	rs = maximum(A)

	# center
	C = collect(-50:1:50)

	# left side
	if ls < -1150
		L  = collect(ls:1:(ls + 99))
		LC = floor(Int64, linspace(L[100], -50, 1002)[2:1001])
	else
		L  = collect(-1150:1:-1051)
		LC = collect(-1050:1:-51)
	end

	# right side
	if rs > 1150
		R  = collect((rs - 99):1:rs)
		CR = floor(Int64, linspace(50, R[1], 1002)[2:1001])
	else
		R  = collect(1051:1:1150)
		CR = collect(51:1:1050)
	end

  pa    = Vector{Int64}(2301)
	pa[:] = vcat(L, LC, C, CR, R)

	return pa
end

### METRICS FUNCTIONS ##########################################################

# Calculate breeding value metrics wrt individuals -----------------------------
function metrics(popmatrix, M, phenotypes = false, patches = "occupied")
	# Given   [popmatrix, phenotypes, occupied]
	# returns [Nx MeanD Meanr V_D V_r C_Dr pID]
	# popmatrix : see init_inds
	# phenotypes: should metrics calculate phenotypes instead of genotypes?
	# patches   : one of the following:
	#             "occupied" : all patches occupied by at least one individual.
	#             "accordian": left edge, center, right edge, & some evenly-
	#                          distributed intermediate patches.
	#             "edges"    : left-most, center-most, & right-most patches only.

	# Nx        : population density in each patch
	# MeanD			: mean dispersal genotype/phenotype in each individual's patch
	# Meanr			: mean low-density growth rate genotype/phenotype
	# V_D       : variance in D
	# V_r       : variance in r
	# C_Dr      : covariance between D and r

  # If the population is extinct, return NULL
  ni = nrow(popmatrix)
	if ni < 1
   	return []
  end

	# get individual-level data
	x  = convert(Vector{Int64}, popmatrix[:X])

  # center traits at zero
	if phenotypes
    D = vector_add(convert( Vector{Float64}, popmatrix[:PD]),
		               -convert(Vector{Float64}, popmatrix[ :D]))
		r = vector_add(convert( Vector{Float64}, popmatrix[:Pr]),
		               -convert(Vector{Float64}, popmatrix[ :r]))
	else
		D = vector_scalar_add(convert(Vector{Float64}, popmatrix[:D]), -M[1])
		r = vector_scalar_add(convert(Vector{Float64}, popmatrix[:r]), -M[2])
	end

	# which patches should be analyzed?
	if patches == "occupied"
    pa = Vector{Int64}()
    sizehint!(pa, ni)
		pa = unique(x)
	elseif patches == "accordian"
		pa = accordian_bins(x)
	elseif patches == "edges"
		pa  = [minimum([minimum(x), -1.0]), 0, maximum([maximum(x), 1.0])]
	end

  # total number of patches
  np = length(pa)

	# make a patch ID dictionary to denote array indices for each patch location
  pID = Dict{Int64, Int64}()
  sizehint!(pID, np)
  for i = 1:np
    pID[pa[i]] = i
  end

  # map each individual's patch to an appropriate array index using pID
	ix = Vector{Int64}(ni)
  map_dict!(ix, x, pID)

  # arrays to store patch-specific metrics
  Nx     = zeros(Int64,   np)
  sumD   = zeros(Float64, np)
  sumr   = zeros(Float64, np)
  spd_D  = zeros(Float64, np)
  spd_r  = zeros(Float64, np)
  spd_Dr = zeros(Float64, np)

  # get patch-specific totals
  running_sum!(sumD, D, ix)
  running_sum!(sumr, r, ix)
  running_count!(Nx, ix)

  # calculate patch-specific means
  Mean_D = sumD ./ Nx
  Mean_r = sumr ./ Nx

  # patch-specific sums of the products of deviations from means
  # (used for calculating variances and covariances)
  sum_prod_deviates!(spd_D,  Mean_D, Mean_D, D, D, ix)
  sum_prod_deviates!(spd_r,  Mean_r, Mean_r, r, r, ix)
  sum_prod_deviates!(spd_Dr, Mean_D, Mean_r, D, r, ix)

  # calculate the uncorrected (co)variances.
  # using uncorrected formula since population mean is (exactly) known.
  V_D  = spd_D  ./ Nx
  V_r  = spd_r  ./ Nx
  C_Dr = spd_Dr ./ Nx

	# center traits on (expected) population mean
	if phenotypes
		genotypes = metrics(popmatrix, M, false, patches)
		vector_add!(Mean_D, convert(Vector{Float64}, genotypes[:Mean_D]))
		vector_add!(Mean_r, convert(Vector{Float64}, genotypes[:Mean_r]))
	else
		vector_scalar_add!(Mean_D, M[1])
		vector_scalar_add!(Mean_r, M[2])
	end

	return DataFrame(X=pa, Nx=Nx, Mean_D=Mean_D, Mean_r=Mean_r, V_D=V_D, V_r=V_r,
	                 C_Dr=C_Dr, pID=collect(pID))
end

### LIFE HISTORY FUNCTIONS #####################################################

# Random mating among individuals that share a patch ---------------------------
function find_mates!(m::Vector{Int64}, pai::Array{Int64, 2},
	                   sx::Array{Int64, 2}, np::Int64)
	# given [m, pai, sx, np] returns [m]
	# m   : a vector of mate matches for each individual
	# pai : an array to track starting and ending indices for individuals that
	#       share a patch. each row corresponds to a row in pa.
	# sx  : and array where the first column is a sorted vector of x and the
	#       second column provides the entry's original index in x.
	# np  : the number of unique occupied patches

	@inbounds for i = 1:np

		# potential mates in patch i
		pm  = sx[pai[i, 1] : pai[i, 2], 2]
		lpm = length(pm)

		if lpm == 1
			# lone individuals do not find mates
			m[sx[pai[i, 1], 2]] = -99
		else
			# individuals sample mates with replacement
			for j = 1:lpm
				m[sx[pai[i, 1] + (j - 1), 2]] = rand(pm[vcat(1:(j - 1), (j + 1) : lpm)])
			end
		end
	end
	nothing
end

# Sexual mating ----------------------------------------------------------------
function mate(popmatrix)
	# Given [popmatrix], returns [m].
	# popmatrix : see init_inds
	# m         : a vector of mate matches (sampled with replacement)

	# If the population is extinct, return NULL
	ni = nrow(popmatrix)
  if ni < 1
  	return []
  end

  # get individual locations
  x  = convert(Vector{Int64}, popmatrix[:X])

	# get a sorted vector of unique occupied patches
	pa = Vector{Int64}()
	sizehint!(pa, ni)
	pa = sort(unique(x))

  # total number of patches
  np = length(pa)

	# make a sorted vector of x and track the original indices
	sx = [x collect(1:ni)]
	sortrowsi!(sx)

	# an array to track starting and ending indices for individuals that share a
	# patch. each row corresponds to a row in pa.
	pai = Array{Int}(np, 2)
	group_index!(pai, sx, pa, np, ni)

  # mate sampling
	m  = Array{Int}(ni)
	find_mates!(m, pai, sx, np)

  return m
end

# Reproduction and dispersal ---------------------------------------------------
function reproduce(Nx::Vector{Int64}, r::Vector{Float64}, m::Vector{Int64},
	                 b::Float64)
  # Given [Nx, r, m, K, a], return [Noff]
	# Nx   : patch density for individual i
	# r    : growth rate for individual i
	# m    : vector of mates
	# b    : Beverton-Holt parameter

	# Noff : the number of offspring produced by individual i

	Noff = Vector{Int64}(length(Nx))
  for i = 1:length(Nx)
		if m[i] == -99
			Noff[i] = 0      # no mate, no offspring
		else
			Noff[i] = rand(Poisson(growth(Nx[i], exp.(r[i]), b)))
		end
  end
	return Noff
end

function disperse(x::Vector{Int64}, λ::Vector{Float64})
	# Given [x, λ], returns [y].
	# x : an integer vector of each individual's current spatial location
	# λ : the Poisson parameter

	# y : an integer vector of new spatial locations for each individual

	lx = length(x)
	y = Vector{Int64}(lx)
	n = Vector{Int64}(lx)

	# since the Poisson only gives positive numbers, flip a coin to determine
	# which direction each individual moves
	d = rand(Bernoulli(), lx)

	for i = 1:lx
		n[i] = rand(Poisson(λ[i]))
		y[i] = x[i] + (d[i] * 2 - 1) * n[i]
	end
	return y
end

function repro_disp(popmatrix, M, V, ρ, H, b, return_dx = false)
	# Given [popmatrix(t), H, V, bw], returns [popmatrix(t+1)]
	# popmatrix : see init_inds
  # V : vector of total phenotypic variance [VD Vr]
  # ρ : vector of additive and non-additive genetic correlations between D and r
	# H : vector of heritabilities [h²D h²r]
	# b : Beverton-Holt parameter

	# If the population is extinct, or if there is only 1 individual, return NULL
	ni = nrow(popmatrix)
	if ni < 2
		return []
  end

	# convert some popmatrix columns to standalone vectors
	x  = convert(Vector{Int64},   popmatrix[:X])
	D  = convert(Vector{Float64}, popmatrix[:D])
	r  = convert(Vector{Float64}, popmatrix[:r])
	Pr = convert(Vector{Float64}, popmatrix[:Pr])

	# Calculate traits through space
  mets = metrics(popmatrix, M)

	# map patch-specific metrics to individuals
  ix  = Vector{Int64}(ni)
	pID = Dict(convert(Vector, mets[:pID]))
	map_dict!(ix, x, pID)

	Nx   = map_valsi(convert(Vector{Int64},  mets[:Nx]),   ix)
	V_D  = map_vals(convert(Vector{Float64}, mets[:V_D]),  ix)
	V_r  = map_vals(convert(Vector{Float64}, mets[:V_r]),  ix)
	C_Dr = map_vals(convert(Vector{Float64}, mets[:C_Dr]), ix)

	# (1) Density-dependent reproduction -----------------------------------------

	m    = mate(popmatrix)         # find a mate (or not) for each individual
	Noff = reproduce(Nx, Pr, m, b) # number of offspring for each parent
	no   = sum(Noff)               # total number of offspring

  # If there are no offspring, return NULL
  if no == 0
		return []
  end

	# (2) Inheritence of genotypes from parents to offspring ---------------------

	# make vectors of values to describe the parents of each offspring
	did   = rep_eachi(collect(1:ni), Noff)  # 'dam' ids
	sid   = rep_eachi(m, Noff)              # 'sire' ids
	dx    = map_valsi(x, did)               # dam location
	dD    = map_vals(D,  did)               # dam dispersal genotype
	sD    = map_vals(D,  sid)               # sire dispersal genotypes
	dr    = map_vals(r,  did)               # dam growth genotype
	sr    = map_vals(r,  sid)               # sire growth genotype

	# get patch-level additive genetic covariance matrix elements for offspring
	oV_D  = rep_each(V_D,  Noff)            # patch-level variance in D
	oV_r  = rep_each(V_r,  Noff)            # patch-level variance in r
	oC_Dr = rep_each(C_Dr, Noff)            # patch-level covariance in D and r

	# multiply elements of additive genetic covariance matrix by 0.5; this is the
	# variance among full-siblings when sampling traits (Roughgarden 1979)
	vector_scalar_mult!(oV_D,  0.5)
	vector_scalar_mult!(oV_r,  0.5)
	vector_scalar_mult!(oC_Dr, 0.5)

  # calculate midparent values for each offspring
	mpvD = vector_add(dD, sD)
	mpvr = vector_add(dr, sr)
	vector_scalar_mult!(mpvD, 0.5)
	vector_scalar_mult!(mpvr, 0.5)

	# initiate some variables
	ΣE   = VE(V, ρ[2], H)        # non-additive (environmental) covariance
	ΣA   = Array{Float64}(2,  2) # additive genetic covariance (placeholder)
	oDr  = Array{Float64}(no, 2) # genotype deviates from mean
	oPDr = Array{Float64}(no, 2) # phenotype deviates from mean

	# loop over all offspring
	for i = 1:no

		ΣA[1, 1] = oV_D[i]
		ΣA[1, 2] = oC_Dr[i]
		ΣA[2, 1] = oC_Dr[i]
		ΣA[2, 2] = oV_r[i]

		# make sure ΣA is positive definite. If not, it's most likely because
		# the matrix is singular, so nudge the corresponding correlation matrix
		# slightly to make the matrix positive definite.
		pos_def_check!(ΣA)

		# sample additive and non-additive deviates
		oDr[i,  1:2] = rand(MvNormal(ΣA))
		oPDr[i, 1:2] = rand(MvNormal(ΣE))
	end

	# sum means and deviates to get trait values
	oD  = vector_add( oDr[:, 1], mpvD)
	or  = vector_add( oDr[:, 2], mpvr)
	oPD = vector_add(oPDr[:, 1], oD)
	oPr = vector_add(oPDr[:, 2], or)

	# (3) Offspring dispersal ----------------------------------------------------
  ox = disperse(dx, exp.(oPD))

	if return_dx
		return DataFrame(X=ox, dX=dx, D=oD, PD=oPD, r=or, Pr=oPr)
	else
		return DataFrame(X=ox, D=oD, PD=oPD, r=or, Pr=oPr)
	end
end

### SIMULATION FUNCTIONS #######################################################

# Save output ------------------------------------------------------------------
function generate_output(popmatrix, p, r, M, H, ρ, i, flag = "full")
	# Given [popmatrix, p, r, H, C, ρ, i, bw], generate summary stats for
	# generation i
	# popmatrix : see init_inds
	# p         : population number for this simulation
	# r         : replicate number for this simulation
	# H         : heritability vector
	# ρ         : correlation vector
	# i         : current generation

	# make blank dataframes for 'zero' and 'time' flags
	zero_DFa = DataFrame(X=-1150:1150,
	                     Nx=0,    Mean_D=0, Mean_r=0, V_D=0, V_r=0, C_Dr=0)
	zero_DFe = DataFrame(X=-1:1,
	                     Nx=0,    Mean_D=0, Mean_r=0, V_D=0, V_r=0, C_Dr=0)
	time_DFe = DataFrame(X=-1:1,
	                     Nx=-999, Mean_D=0, Mean_r=0, V_D=0, V_r=0, C_Dr=0)

	if flag == "full"
    # simulation is running as planned, get metrics
		ni = nrow(popmatrix)                           # population size
	  ag = metrics(popmatrix, M, false, "accordian") # genotype summary metrics
		ap = metrics(popmatrix, M, true,  "accordian") # phenotype summary metrics
		eg = metrics(popmatrix, M, false, "edges")     # edge patches gen. summary
		ep = metrics(popmatrix, M, true,  "edges")     # edge patches phen. summary
	elseif flag == "zero"
		# population has gone extinct, write a bunch of zeros
		ni = 0
	  ag = zero_DFa
		ap = zero_DFa
		eg = zero_DFe
		ep = zero_DFe
	elseif flag == "timeout"
		# simulation has timed out, write some nonsense-values
		ni = -i
	  ag = zero_DFa
		ap = zero_DFa
		eg = time_DFe
		ep = time_DFe
	end

	# parameter info
	parameters = repmat([p, r, M[1], M[2], H[1], H[2], ρ[1], ρ[2], i, ni]', 12, 1)

	# reshaped summary statistics from sum_metrics and LE_sum_metrics
  row_id    = ["patch", "N",
	             "Mean_D",  "Mean_r",  "V_D",  "V_r",  "C_Dr",
							 "Mean_PD", "Mean_Pr", "V_PD", "V_Pr", "C_PDr"]
	edges     = vcat(eg[:X]', eg[:Nx]',
	                 eg[:Mean_D]', eg[:Mean_r]', eg[:V_D]', eg[:V_r]', eg[:C_Dr]',
								   ep[:Mean_D]', ep[:Mean_r]', ep[:V_D]', ep[:V_r]', ep[:C_Dr]')
	accordian = vcat(ag[:X]', ag[:Nx]',
	                 ag[:Mean_D]', ag[:Mean_r]', ag[:V_D]', ag[:V_r]', ag[:C_Dr]',
								   ap[:Mean_D]', ap[:Mean_r]', ap[:V_D]', ap[:V_r]', ap[:C_Dr]')

	return hcat(row_id, parameters, edges, accordian)
end

# Run the simulation -----------------------------------------------------------
function runsim(popmatrix, ngens, M, V, ρ, H, b, p, r, bat_time)
	# Given [popmatrix, ngens, M, V, ρ, H, b, p, r, bat_time]
	# returns changes in population size, extent, and genetic variance over time.
  # popmatrix : see init_inds
	# n  	      : initial number of starting individuals
	# spX       : initial spatial locations
	# ngens     : number of generations to simulate
	# M         : vector of mean phenotypes [μD μK μr]
	# V         : vector of total phenotypic variance [VD VK Vr]
	# ρ         : vector of correlations [ρA ρE]
	# H         : vector of heritabilities [h²D h²r]
	# b         : Beverton-Holt parameter
	# p         : population number for this generation
	# r         : replicate number for this simulation
	# bat_time  : time that the current batch process started

	# start timers
	tic()
	sim_time = time()
	sim_flag = ""
	met_flag = "full"

	# set a unique seed based on H, p, and r, so simulations can be repeated.
	srand(plant_seed(ρ, H, r, p))

	# initialize output
	batchoutput = Array{Any}(ngens * 12, 2315)

	# loop over generations
  for i = 1:ngens

		# an array index for the output
		ind = (1:12) + (12) * (i - 1)

		if met_flag == "full"
			# do reproduction and dispersal
			popmatrix = repro_disp(popmatrix, M, V, ρ, H, b)
		end

		# check if future generations should be simulated
		if popmatrix == []
			# stop simulating if population went extinct
			met_flag = "zero"
		elseif (time() - bat_time) > 16.0 * 3600.0
			# stop simulating if the batch has run longer than 16 hours
			met_flag = "timeout"
			sim_flag = "!BATCH_RUN_TIMEOUT!"
		elseif (time() - sim_time) > 10.0 * 3600.0
			# stop simulating if the simulation has run longer than 10 hours
			met_flag = "timeout"
			sim_flag = "!SIM_RUN_TIMEOUT!"
		else
		end

		# record data
		batchoutput[ind, :] = generate_output(popmatrix, p, r, M, H, ρ, i, met_flag)
	end

	# track simulation info and run times
	print(  "rep ",       (@sprintf "%2.0f" r),
				", pop ",       (@sprintf "%2.0f" p),
				", h²D ",       (@sprintf "%0.2f" H[1]),
				", h²r ",       (@sprintf "%0.2f" H[2]),
				", ρA ",        (@sprintf "%5.2f" ρ[1]),
				", ρE ",        (@sprintf "%5.2f" ρ[2]),
				", time ",      (clocktime(toq())),
				", final pop ", (@sprintf "%-9.0f" batchoutput[end, 11]),
				" | total time ", (clocktime(time() - bat_time)),
				" | ", sim_flag, "\n")

	return batchoutput
end

max_stars_repo_name: bochocki/correlatedtraits

max_issues_repo_name: bochocki/correlatedtraits

max_forks_repo_name: bochocki/correlatedtraits


text:
using HermiteGF
using Test

include("trapz.jl")
include("test_interpolation1d.jl")
include("test_interpolation2d.jl")
include("test_interpolation3d.jl")
include("test_interpolation4d.jl")
include("test_interpolation5d.jl")


text:
  cd(@__DIR__)
  using Pkg
  pkg"activate ."

  function main()
    include(joinpath("src", "TODOList.jl"))
  end

  main()

max_stars_repo_name: ky-s/MyGenieApp

max_issues_repo_name: ky-s/MyGenieApp

max_forks_repo_name: ky-s/MyGenieApp


text:
# Mel Generalized-Log Spectrum Approximation Digital Filter

# MGLSABaseFilter represents a base filter of MGLSAFilter.
mutable struct MGLSABaseFilter <: Filter
    order::Int   # order of mel generalized cepstrum
    α::Float64  # all-pass constant
    delay::Vector{Float64}

    function MGLSABaseFilter(order::Int, α::Float64)
        new(order, α, zeros(order+1))
    end
end

allpass_alpha(f::MGLSABaseFilter) = f.α

function filt!(bf::MGLSABaseFilter, x::Float64, coef::Vector{Float64})
    d = delay(bf)
    y = d[1] * coef[2]
    α = allpass_alpha(bf)

    for i=2:length(coef)-1
        @inbounds d[i] += α*(d[i+1] - d[i-1])
        @inbounds y += d[i] * coef[i+1]
    end

    result = x - y

    # t <- t+1 in time
    for i=length(d):-1:2
        @inbounds d[i] = d[i-1]
    end
    d[1] = α*d[1] + (1.0 - α*α)*result

    result
end

# MGLSAFilter represents a Mel Generalized Log Spectrum Approximation Digital
# Filter.
mutable struct MGLSADF <: MelGeneralizedCepstrumSynthesisFilter
    filters::Vector{MGLSABaseFilter}

    function MGLSADF(order::Int, α::Float64, nstage::Int)
        @assert nstage > 0
        filters = Vector{MGLSABaseFilter}(nstage)
        for i=1:length(filters)
            filters[i] = MGLSABaseFilter(order, α)
        end
        new(filters)
    end
end

allpass_alpha(f::MGLSADF) = allpass_alpha(f.filters[1])
glog_gamma(f::MGLSADF) = -1.0/length(f.filters)
nstage(f::MGLSADF) = length(f.filters)

function filt!(f::MGLSADF, x::Float64, coef::Vector{Float64})
    y = x
    for i=1:length(f.filters)
        y = filt!(f.filters[i], y, coef)
    end
    y
end

function to_filtcoef{T<:MelGeneralizedCepstrum}(f::MGLSADF,
                                                mgc::SpectralParamState{T})
    def = paramdef(mgc)
    @assert allpass_alpha(def) == allpass_alpha(f)
    @assert glog_gamma(def) == glog_gamma(f)
    rawdata(mgc2b(mgc))
end


text:
using Convex, LinearAlgebra
using UnbalancedOptimalTransport: AbstractDivergence, KL, TV, RG, Balanced

"""
    divergence(D::AbstractDivergence, a, b)

Computes the the divergence `D` between `a` and `b`, where `a` can be a
Convex.jl expression, and `b` is a constant.
"""
function divergence end

lg(x) = x <= 0 ? zero(x) : log(x)

# need to type-pirate a method for non-Convex.jl objects
Convex.entropy(a::AbstractArray) = sum(x -> -x * lg(x), a)

function divergence(::KL{ρ}, a, b) where {ρ}
    ρ * (-entropy(a) - dot(a, lg.(b)) + sum(b - a))
end

function divergence(::TV{ρ}, a, b) where {ρ}
    ρ * norm(a - b, 1)
end

# Computes the objective of the unbalanced regularized optimal transport problem
# This is a separate method so it can be reused for testing `optimal_coupling`.
function objective(
    π,
    D::AbstractDivergence,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ;
    C = (x, y) -> norm(x - y),
)
    π_1 = sum(π, dims = 2)
    π_2 = vec(sum(π, dims = 1))

    # instantiate the cost matrix
    C_matrix = [C(xᵢ, yⱼ) for xᵢ in a.set, yⱼ in b.set]

    obj = dot(C_matrix, π) + ϵ * divergence(KL(), vec(π), kron(b.density, a.density))

    if (D isa KL) || (D isa TV)
        obj += divergence(D, π_1, a.density) + divergence(D, π_2, b.density)
    end

    return obj
end

function OT_convex(
    D::AbstractDivergence,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ;
    solver,
    verbose = false,
    C = (x, y) -> norm(x - y),
)
    ϵ >= 0 || throw(DomainError(ϵ, "Need ϵ >= 0"))

    # coupling and its marginals
    π = Variable(length(a), length(b))
    π_1 = sum(π, dims = 2)
    π_2 = vec(sum(π, dims = 1))

    obj = objective(π, D, a, b, ϵ; C = C)

    if D isa Balanced
        constraints = [π >= 0, π_1 == a.density, π_2 == b.density]
    elseif D isa RG
        l, u = params(D)
        constraints = [
            π >= 0,
            l * a.density <= π_1,
            π_1 <= u * a.density,
            l * b.density <= π_2,
            π_2 <= u * b.density,
        ]
    else
        constraints = [π >= 0]
    end

    problem = minimize(obj, constraints)

    solve!(problem, solver; verbose = verbose)

    return (optimal_value = problem.optval, optimal_coupling = evaluate(π))
end

params(::RG{l,u}) where {l,u} = (l, u)


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/LGforcing.jl")

    add_comp!(m, LGforcing)

    set_param!(m, :LGforcing, :c_LGconcentration, readpagedata(m,"test/validationdata/$valdir/c_LGconcentration.csv"))

    # run Model
    run(m)

    forcing=m[:LGforcing,:f_LGforcing]
    forcing_compare=readpagedata(m,"test/validationdata/$valdir/f_LGforcing.csv")

    @test forcing ≈ forcing_compare rtol=1e-3
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
module DestructuringKoans

#=
Links:
    - https://docs.julialang.org/en/v1/manual/functions/#Multiple-Return-Values-1
    - https://docs.julialang.org/en/v1/manual/functions/#Argument-destructuring-1
=#

function add_separately(x, y)
    x_real = real(x)
    x_im = imag(x)
    y_real = real(y)
    y_im = imag(y)
    x_real + y_real, x_im + y_im
end

function add_coordinates((x1, y1), (x2, y2))
    x1 + x2, y1 + y2
end

end


text:
### Abstract MALA state

abstract MALAState{F<:VariateForm} <: LMCSamplerState{F}

### MALA state subtypes

## UnvMALAState holds the internal state ("local variables") of the MALA sampler for univariate parameters

type UnvMALAState <: MALAState{Univariate}
  pstate::ParameterState{Continuous, Univariate} # Parameter state used internally by MALA
  tune::MCTunerState
  ratio::Real
  μ::Real

  function UnvMALAState(pstate::ParameterState{Continuous, Univariate}, tune::MCTunerState, ratio::Real, μ::Real)
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    new(pstate, tune, ratio, μ)
  end
end

UnvMALAState(pstate::ParameterState{Continuous, Univariate}, tune::MCTunerState=BasicMCTune()) =
  UnvMALAState(pstate, tune, NaN, NaN)

## MuvMALAState holds the internal state ("local variables") of the MALA sampler for multivariate parameters

type MuvMALAState <: MALAState{Multivariate}
  pstate::ParameterState{Continuous, Multivariate} # Parameter state used internally by MALA
  tune::MCTunerState
  ratio::Real
  μ::RealVector

  function MuvMALAState(pstate::ParameterState{Continuous, Multivariate}, tune::MCTunerState, ratio::Real, μ::RealVector)
    if !isnan(ratio)
      @assert ratio > 0 "Acceptance ratio should be positive"
    end
    new(pstate, tune, ratio, μ)
  end
end

MuvMALAState(pstate::ParameterState{Continuous, Multivariate}, tune::MCTunerState=BasicMCTune()) =
  MuvMALAState(pstate, tune, NaN, Array(eltype(pstate), pstate.size))

### Metropolis-adjusted Langevin Algorithm (MALA)

immutable MALA <: LMCSampler
  driftstep::Real

  function MALA(driftstep::Real)
    @assert driftstep > 0 "Drift step is not positive"
    new(driftstep)
  end
end

MALA() = MALA(1.)

### Initialize MALA sampler

## Initialize parameter state

function initialize!{F<:VariateForm}(
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::MALA
)
  parameter.uptogradlogtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
  @assert all(isfinite(pstate.gradlogtarget)) "Gradient of log-target not finite: initial values out of support"
end

## Initialize MALA state

sampler_state(
  parameter::Parameter{Continuous, Univariate},
  sampler::MALA,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Univariate},
  vstate::VariableStateVector
) =
  UnvMALAState(generate_empty(pstate), tuner_state(parameter, sampler, tuner))

sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::MALA,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector
) =
  MuvMALAState(generate_empty(pstate), tuner_state(parameter, sampler, tuner))

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Univariate},
  x::Real,
  parameter::Parameter{Continuous, Univariate},
  sampler::MALA
)
  pstate.value = x
  parameter.uptogradlogtarget!(pstate)
end

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::MALA
)
  pstate.value = copy(x)
  parameter.uptogradlogtarget!(pstate)
end

Base.show(io::IO, sampler::MALA) = print(io, "MALA sampler: drift step = $(sampler.driftstep)")

Base.writemime(io::IO, ::MIME"text/plain", sampler::MALA) = show(io, sampler)


text:
function coldist(Z,i,j,N)
    s = 0
    @inbounds @simd for k in 1:N
        s += Z[k,i] != Z[k,j]
    end
    s
end

function compute_theta(Z::Matrix{Ti}) where Ti <: Integer
    N,M = size(Z)
    meanfracid = @distributed (+) for i = 1:M-1
        s = 0.0
        for j = i+1:M
            s += 1.0 - coldist(Z,i,j,N) / N
        end
        s
    end
    meanfracid /= 0.5 * M * (M-1)
    theta = min(0.5, 0.38 * 0.32 / meanfracid)
    return theta
end

"""
    compute_weights(Z::Matrix{Ti}, theta::Real) where Ti <: Integer
Compute the normalized counts of the number of sequences at hamming distance ≤ `theta` from any given sequence
in `Z`.
"""
function compute_weights(Z::Matrix{Ti}, theta::Real) where Ti <: Integer
    N, M = size(Z)
    thresh = floor(N * theta)
    W = zeros(M) |> SharedArray
    @sync @distributed for i in 1:M
        @simd for j in 1:M
            dij = coldist(Z,i,j,N)
            val = (dij < thresh)
            W[i] += val
        end
    end

    @sync @distributed for i in 1:M
        W[i] = 1.0 / W[i]
    end
    Meff = sum(W)
    println("M = $M N = $N Meff = $Meff theta = $theta")
    return sdata(W)./Meff , Meff
end

"""
    compute_weights(Z::Matrix{Ti}, theta::Symbol) where Ti<:Integer
Compute the normalized counts of the number of sequences at hamming distance ≤ of a precomputed optimal threshold.
See `Fast and accurate multivariate Gaussian modeling of protein families: Predicting residue contacts and
protein-interaction partners` [https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0092721]
and in particular the Supplementary Information at section `2 Reweighting Scheme` for details.
"""
function compute_weights(Z::Matrix{Ti}, theta::Symbol) where Ti<:Integer
    theta !== :auto && error("only theta=:auto is supported")
    thr = compute_theta(Z)
    return compute_weights(Z,thr)
end

max_stars_repo_name: pagnani/CorrDCA

max_issues_repo_name: pagnani/CorrDCA

max_forks_repo_name: pagnani/CorrDCA


text:
"""
   update_opinion!(state, agent_idx, config)

Update an agent's opinion depending on its perceived public opinion

# Arguments
- `state`: a tuple of the current graph and agent_list
- `agent_idx`: agent index
- `config`: Config object as provided by Config()

See also: [`Config`](@ref), [`update__perceiv_publ_opinion!`](@ref)
"""
function update_opinion!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    config::Config
)

    agent_list = state[2]
    this_agent = agent_list[agent_idx]

    if (
        abs(this_agent.opinion - this_agent.perceiv_publ_opinion)
        < config.opinion_threshs.backfire
    )
        this_agent.opinion = (
            config.agent_props.own_opinion_weight * this_agent.opinion
            + (
                (1 - config.agent_props.own_opinion_weight)
                * this_agent.perceiv_publ_opinion
            )
        )
    else
        if (
            this_agent.opinion * this_agent.perceiv_publ_opinion > 0
            && (
                abs(this_agent.opinion) - abs(this_agent.perceiv_publ_opinion)
                < 0
            )
        )
            this_agent.opinion = (
                config.agent_props.own_opinion_weight * this_agent.opinion
            )
        else
            this_agent.opinion = (
                (2 - config.agent_props.own_opinion_weight)
                * this_agent.opinion
            )
        end

        if this_agent.opinion > 1
            this_agent.opinion = 1
        elseif this_agent.opinion < -1
            this_agent.opinion = -1
        end
    end

    return state
end

# suppress output of include()
;


text:
include("../src/includes.jl")

F = rand(2,2)
D = rand(3,3)

qaph = HopfieldQAP(D, F)

# using Main.HopfieldNets
using Test
using LinearAlgebra
# using Pkg
# include(Pkg.dir("Main.HopfieldNets", "demo", "letters.jl"))
include("../demo/letters.jl")

patterns = hcat(X, O)

n = size(patterns, 1)

dnet = DiscreteHopfieldNet(n)
cnet = ContinuousHopfieldNet(n)

threex3 = DiscreteHopfieldNet((3,3))
twox3 = DiscreteHopfieldNet((2,3))

# N = (2,3)

# a = Vector{Int}()

# for i in N
#      push!(a, repeat([i], length(N))...)
# end
# a

cycle!(twox3)

print(threex3.s[1,1])
update!(threex3, (1,1))
print(threex3.s[1,1])

# N = size(threex3.s)
# N_2 = size(threex3.W)
# randInd = [shuffle(1:i) for i in N]

# D = rand(Float64, 2, 2)
# F = rand(Float64, 3, 3)
# A = B = 2

# # D & F are square matrices always
# N = (size(D)[1],size(F)[1]) #dimensions of the Hopfield Network
# tnet = DiscreteHopfieldNet(N)

# iters = Base.Iterators.product([collect(1:i) for i in N]...)
# for i in iters
#     for j in iters
#           tnet.W[i...,j...] = -2()

# N = size(twox3.s)
# N_2 = size(twox3.W)
# randInd = [shuffle(1:i) for i in N]

for net in [dnet, cnet]
    train!(net, patterns)

    e0 = energy(net)
    settle!(net, 1_000, false)
    eFinal = energy(net)
    @assert e0 != eFinal

    Xcorrupt = copy(X)
    for i = 2:7
         Xcorrupt[i] = 1
    end
    Xrestored = associate!(net, Xcorrupt)
    @test norm(Xcorrupt - Xrestored) > 1e-4
    @test norm(X - Xrestored) < 1e-4

    Ocorrupt = copy(O)
    for i = 2:7
         Ocorrupt[i] = -1
    end
    Orestored = associate!(net, Ocorrupt)
    @test norm(Ocorrupt - Orestored) > 1e-4
    @test norm(O - Orestored) < 1e-4
end


text:
"""
Created on Fri 15 May 2020

@author: Yoann Pradat

Functions for simulation different matrices for the factorization V =  WH + E
"""

using Combinatorics
using LinearAlgebra

function scale_col(X; p=1)
    norms_X = mapslices(x -> norm(x, p), X, dims=1)
    X .* repeat(norms_X .^ -1, size(X, 1), 1)
end

function random_fixed_total(n::Integer, total::Integer, rng::AbstractRNG)
    μ = total÷n
    T = zeros(Integer, n)
    T[1] = total

    while sum(T[1:n-1]) >= total
        T[1:n-1] .= rand(rng, 1:2μ,  n-1)
    end

    T[n] = total - sum(T[1:n-1])

    T
end

"""
    simulate_X_clustered_a(K, N, n_clusters, dist_X, rng)

Exaclty one cluster of rows is active in one cluster of columns. Number of clusters of rows is equal to the number of
clusters of columns.
"""
function simulate_X_clustered_a(K::Integer, N::Integer, n_clusters::Integer, dist_X::Distribution, rng::AbstractRNG)
    X = zeros(Float64, K,N)

    lims_K = cumsum(random_fixed_total(n_clusters, K, rng))
    lims_N = cumsum(random_fixed_total(n_clusters, N, rng))
    pushfirst!(lims_K, 0)
    pushfirst!(lims_N, 0)

    for i = 1:n_clusters
        range_K = (lims_K[i]+1):lims_K[i+1]
        range_N = (lims_N[i]+1):lims_N[i+1]
        X[range_K, range_N] .= rand(rng, dist_X, length(range_K), length(range_N))
    end

    X, lims_K, lims_N, Matrix{Int}(I,n_clusters, n_clusters)
end

# """
#     simulate_X_clustered_b(K, N, n_clusters, dist_X, rng)
#
# One or more clusters of rows are active in one cluster of columns. The number of clusters of rows is equal to the number
# of clusters of columns.
# """
# function simulate_X_clustered_b(K::Integer, N::Integer, n_clusters::Integer, dist_X::Distribution, rng::AbstractRNG)
#     X = zeros(Float64, K,N)
#
#     lims_K = cumsum(random_fixed_total(n_clusters, K, rng))
#     lims_N = cumsum(random_fixed_total(n_clusters, N, rng))
#     pushfirst!(lims_K, 0)
#     pushfirst!(lims_N, 0)
#
#     active_clusters = zeros(Int, n_clusters, n_clusters)
#     combinations_K  = collect(combinations(1:n_clusters))
#     randperm_comb_K = randperm(rng, 2^n_clusters-1)
#     for j = 1:n_clusters
#         active_clusters[combinations_K[randperm_comb_K[j]], j] .= 1
#     end
#
#     for j = 1:n_clusters
#         range_N = (lims_N[j]+1):lims_N[j+1]
#         for i in 1:n_clusters
#             range_K = (lims_K[i]+1):lims_K[i+1]
#             if active_clusters[i,j] == 1
#                 X[range_K, range_N] .= rand(rng, dist_X, length(range_K), length(range_N))
#             end
#         end
#     end
#
#     X, lims_K, lims_N, active_clusters
# end

"""
    simulate_X_clustered_b(K, N, n_clusters_K, n_clusters_N, dist_X, rng)

One or more clusters of rows are active in one cluster of columns. Number of clusters of rows is inferior or equal to
the number of clusters of columns.
"""
function simulate_X_clustered_b(K::Integer, N::Integer, n_clusters_K::Integer, n_clusters_N::Integer, dist_X::Distribution, rng::AbstractRNG)
    X = zeros(Float64, K, N)

    n_clusters_K <= n_clusters_N || throw(ArgumentError("Please choose n_clusters_K <= n_clusters_N"))

    lims_K = cumsum(random_fixed_total(n_clusters_K, K, rng))
    lims_N = cumsum(random_fixed_total(n_clusters_N, N, rng))
    pushfirst!(lims_K, 0)
    pushfirst!(lims_N, 0)

    active_clusters = zeros(Int, n_clusters_K, n_clusters_N)
    combinations_K  = collect(combinations(1:n_clusters_K))
    randperm_comb_K = randperm(rng, 2^n_clusters_K-1)
    for j = 1:n_clusters_N
        active_clusters[combinations_K[randperm_comb_K[j]], j] .= 1
    end

    for j = 1:n_clusters_N
        range_N = (lims_N[j]+1):lims_N[j+1]
        for i in 1:n_clusters_K
            range_K = (lims_K[i]+1):lims_K[i+1]
            if active_clusters[i,j] == 1
                X[range_K, range_N] .= rand(rng, dist_X, length(range_K), length(range_N))
            end
        end
    end

    X, lims_K, lims_N, active_clusters
end

"""
    simulate_WH_clustered(; F, N, K, n_clusters_K, n_clusters_N, dist_W, dist_H, H_cluster_mode, rng)

Simulate W and H with clustering on H.
"""
function simulate_WH_clustered(;F::Integer, N::Integer, K::Integer, n_clusters_K::Integer=0, n_clusters_N::Integer, dist_W::Distribution, dist_H::Distribution, H_cluster_mode::Symbol, rng::AbstractRNG)
    K = K

    W = rand(rng, dist_W, F, K)
    W = scale_col(W, p=1)

    if H_cluster_mode == :a
        H, lims_K, lims_N, active_clusters = simulate_X_clustered_a(K, N, n_clusters_N, dist_H, rng)
    elseif H_cluster_mode == :b
        H, lims_K, lims_N, active_clusters = simulate_X_clustered_c(K, N, n_clusters_K, n_clusters_N, dist_H, rng)
    end

    W, H, lims_K, lims_N, active_clusters
end

"""
    simulate_WH(F, N; n_factors, dist_W, dist_H, rng)

Simulate W and H.
"""
function simulate_WH(;F::Integer, N::Integer, n_factors::Integer, dist_W::Distribution, dist_H::Distribution, rng::AbstractRNG)
    K = n_factors

    W = rand(rng, dist_W, F, K)
    W = scale_col(W, p=1)

    H = rand(rng, dist_H, K, N)

    W, H
end


text:
using Documenter
using Molly

makedocs(
    format = :html,
    sitename = "Molly.jl",
    pages = [
        "Home"                      => "index.md",
        "Documentation"             => "docs.md",
        "Differentiable simulation" => "differentiable.md",
        "Examples"                  => "examples.md",
        "API"                       => "api.md",
    ],
    authors = "Joe G Greener"
)

deploydocs(
    repo = "github.com/JuliaMolSim/Molly.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
function calculate_initial_conditions!(sim::Simulation, inputs::SimulationInputs)
    @debug "Start state intialization routine"
    sys = get_system(sim)
    initial_guess = sim.x0_init
    res = PSY.solve_powerflow!(sys)
    if !res
        sim.status = BUILD_FAILED
        @error("PowerFlow failed to solve")
        return false
    end
    var_count = get_variable_count(inputs)

    @debug "Setting up initilization indexing"
    bus_size = get_bus_count(inputs)
    injection_start = get_injection_pointer(inputs)
    branches_start = get_branches_pointer(inputs)

    @debug "Updating Voltage guess"
    for bus in PSY.get_components(PSY.Bus, sys)
        @debug PSY.get_name(bus)
        @debug "V_r" PSY.get_magnitude(bus) * cos(PSY.get_angle(bus))
        @debug "V_i" PSY.get_magnitude(bus) * sin(PSY.get_angle(bus))
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(inputs)[bus_n]
        initial_guess[bus_ix] = PSY.get_magnitude(bus) * cos(PSY.get_angle(bus))
        initial_guess[bus_ix + bus_size] = PSY.get_magnitude(bus) * sin(PSY.get_angle(bus))
    end

    @debug "Updating Source internal voltages"
    sources = PSY.get_components(PSY.Source, sys)
    if !isempty(sources)
        for s in sources
            initialize_device!(s)
        end
    end

    @debug "Updating Dynamic Injection Component Initial Guess"
    for d in get_injectors_data(inputs)
        dynamic_device = PSY.get_dynamic_injector(d)
        @debug PSY.get_name(d) typeof(d) typeof(dynamic_device)
        n_states = PSY.get_n_states(dynamic_device)
        ix_range = range(injection_start, length = n_states)
        injection_start = injection_start + n_states
        x0_device = initialize_device!(d)
        @assert length(x0_device) == n_states
        initial_guess[ix_range] = x0_device
    end

    @debug "Updating Component Initial Guess"
    dyn_branches = PSY.get_components(PSY.DynamicBranch, sys)
    if !isempty(dyn_branches)
        for br in dyn_branches
            @debug PSY.get_name(br) typeof(br)
            n_states = PSY.get_n_states(br)
            ix_range = range(branches_start, length = n_states)
            branches_start = branches_start + n_states
            x0_branch = initialize_device!(br)
            @assert length(x0_branch) == n_states
            initial_guess[ix_range] = x0_branch
        end
    else
        @debug "No Dynamic Branches in the system"
    end

    if any(!isfinite, initial_guess)
        i = findall(!isfinite, initial_guess)
        invalid_initial_guess = String[]
        for (device, states) in get_global_index(sim)
            for state in states
                if state.second ∈ i
                    push!(invalid_initial_guess, "$device - $(p.first)")
                end
            end
        end
        error("Invalid initial guess values $invalid_initial_guess")
        sim.status = BUILD_FAILED
        return false
    end

    dx0 = zeros(var_count) #Define a vector of zeros for the derivative
    inif! = (out, x) -> system!(
        out,    #output of the function
        dx0,    #derivatives equal to zero
        x,      #states
        inputs,    #Parameters
        0.0,    #time equals to zero.
    )
    #Refine initial solution
    @debug "Start NLSolve System Run"
    @debug initial_guess
    try
        sys_solve = NLsolve.nlsolve(
            inif!,
            initial_guess,
            xtol = STRICT_NL_SOLVE_TOLERANCE,
            ftol = STRICT_NL_SOLVE_TOLERANCE,
            method = :trust_region,
        ) #Solve using initial guess x0
    catch e
        bt = catch_backtrace()
        @error "NLsolve failed to solve" exception = e, bt
        sim.status = BUILD_FAILED
        return false
    end
    if !NLsolve.converged(sys_solve)
        @warn(
            "Initialization failed, initial conditions do not meet conditions for an stable equilibrium.\nTrying to solve again reducing numeric tolerance from $(STRICT_NL_SOLVE_TOLERANCE):"
        )
        sys_solve = NLsolve.nlsolve(
            inif!,
            initial_guess,
            xtol = RELAXED_NL_SOLVE_TOLERANCE,
            ftol = RELAXED_NL_SOLVE_TOLERANCE,
            method = :trust_region,
        ) #Solve using initial guess x0
        if NLsolve.converged(sys_solve)
            @info(
                "Initialization succeeded with a relaxed tolerance of $(RELAXED_NL_SOLVE_TOLERANCE). Saving solution"
            )
        else
            @warn(
                "Initialization failed again. Initial conditions do not meet conditions for an stable equilibrium\nSaving best result."
            )
        end
    end
    @debug "Write result to initial guess vector"
    initial_guess .= sys_solve.zero
    return NLsolve.converged(sys_solve)
end

"""
Returns a Dictionary with the resulting initial conditions of the simulation
"""
function get_initial_conditions(sim::Simulation)
    system = get_system(sim)
    bus_size = get_bus_count(sim.simulation_inputs)
    V_R = Dict{Int, Float64}()
    V_I = Dict{Int, Float64}()
    Vm = Dict{Int, Float64}()
    θ = Dict{Int, Float64}()
    for bus in PSY.get_components(PSY.Bus, system)
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(sim.simulation_inputs)[bus_n]
        V_R[bus_n] = sim.x0_init[bus_ix]
        V_I[bus_n] = sim.x0_init[bus_ix + bus_size]
        Vm[bus_n] = sqrt(sim.x0_init[bus_ix]^2 + sim.x0_init[bus_ix + bus_size]^2)
        θ[bus_n] = angle(sim.x0_init[bus_ix] + sim.x0_init[bus_ix + bus_size] * 1im)
    end
    results = Dict{String, Any}("V_R" => V_R, "V_I" => V_I, "Vm" => Vm, "θ" => θ)
    for device in PSY.get_components(PSY.DynamicInjection, system)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        global_index = get_global_index(sim.simulation_inputs)[name]
        x0_device = Dict{Symbol, Float64}()
        for s in states
            x0_device[s] = sim.x0_init[global_index[s]]
        end
        results[name] = x0_device
    end
    dyn_branches = PSY.get_components(PSY.DynamicBranch, system)
    if !isempty(dyn_branches)
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            global_index = get_global_index(sim.simulation_inputs)[name]
            x0_br = Dict{Symbol, Float64}()
            for s in states
                x0_br[s] = sim.x0_init[global_index[s]]
            end
            printed_name = "Line " * name
            results[printed_name] = x0_br
        end
    end
    return results
end


text:
module execute
using LPDM
using POMDPs, Parameters, StaticArrays, D3Trees, Distributions, SparseArrays
using Combinatorics
using StatsFuns
using Printf
using Dates

include("LightDarkLPDM.jl")
using LightDarkLPDM

include("1D/LPDMBounds1d.jl")
include("2D/LPDMBounds2d.jl")

struct ProblemConfig
    n_bins                  ::Int64
    max_actions             ::Int64
    n_new_actions           ::Int64 # number of new actions to select at the same time in SA
    action_range_fraction   ::Float64
    max_exploit_visits      ::Int64
    max_belief_clusters     ::Int64
end

struct LPDMTest
    action_mode         ::Symbol
    obs_mode            ::Symbol
    reward_mode         ::Symbol
    action_std          ::Float64
    problem_config      ::ProblemConfig
    n_sims              ::Int64
end

struct LPDMScenario{S}
    s0::S
end

# reward function options - :quadratic or :fixed
function batch_execute(;dims::Int64=1, debugger::Bool = false)

    # General test parameters
    reward_mode                 = :quadratic
    n_sims                      = debugger ? 1 : 100
    vis::Vector{Int64}          = debugger ? Int64[1] : Int64[]
    debug                       = debugger ? 3 : 0

    # 1D problem configuration
    n_bins1d                    = 10
    max_actions1d               = debugger ? 10 : 50
    n_new_actions1d             = 2
    action_range_fraction1d     = 0.25
    max_exploit_visits1d        = 0
    max_belief_clusters1d       = 10

    # 2D problem configuration
    n_bins2d                    = 10 # per dimension
    max_actions2d               = 100
    n_new_actions2d             = 2
    action_range_fraction2d     = 0.25
    max_exploit_visits2d        = 0
    max_belief_clusters2d       = 10

    pconfig1d = ProblemConfig(n_bins1d, max_actions1d, n_new_actions1d, action_range_fraction1d, max_exploit_visits1d, max_belief_clusters1d)
    pconfig2d = ProblemConfig(n_bins2d, max_actions2d, n_new_actions2d, action_range_fraction2d, max_exploit_visits2d, max_belief_clusters2d)

    # 1D solver configuration (some are left as defaults)
    sconfig1d = LPDMConfig(
            search_depth        = 30,
            seed                = 0xffffffff, # assigned per scenario (will cause an error if not assigned)
            time_per_move       = -1.0, # should normally be 1.0
            n_particles         = 25,
            sim_len             = debugger ? 1 : -1,
            max_trials          = 10000,
            debug               = debug,
            action_mode         = :tbd,     # assigned per test (will cause an error if not assigned)
            obs_mode            = :tbd)     # assigned per test (will cause an error if not assigned)

    sconfig1d_debug = LPDMConfig(
            search_depth        = 2,
            seed                = 0xffffffff, # assigned per scenario (will cause an error if not assigned)
            time_per_move       = 1.0, # should normally be 1.0
            n_particles         = 10,
            sim_len             = 1,
            max_trials          = 3,
            debug               = debug,
            action_mode         = :tbd,     # assigned per test (will cause an error if not assigned)
            obs_mode            = :tbd)     # assigned per test (will cause an error if not assigned)

    # 2D solver configuration
    sconfig2d = LPDMConfig(
            search_depth        = 30,
            seed                = 0xffffffff,  # assigned per scenario (will cause an error if not assigned)
            time_per_move       = 2.0,
            # time_per_move       = -1.0,
            n_particles         = 50,
            sim_len             = debugger ? 1 : -1,
            max_trials          = -1,
            # max_trials          = 100,
            debug               = debug,
            action_mode         = :tbd,     # assigned per test (will cause an error if not assigned)
            obs_mode            = :tbd)     # assigned per test (will cause an error if not assigned)

    if dims == 1
        S = LD1State
        A = LD1Action
        O = LD1Obs
        B = LDBounds1d{S,A,O}
        pconfig = pconfig1d
        sconfig = debugger ? sconfig1d_debug : sconfig1d
    elseif dims == 2
        S = LD2State
        A = LD2Action
        O = LD2Obs
        B = LDBounds2d{S,A,O}
        pconfig = pconfig2d
        sconfig = sconfig2d
    else
        error("Invalid number of dimensions $dims")
    end

    test = Array{LPDMTest}(undef,0)

    # DISCRETE OBSERVATIONS
    # push!(test, LPDMTest(:standard, :discrete, reward_mode, pconfig, n_sims))
    # push!(test, LPDMTest(:extended, :discrete, reward_mode, pconfig, n_sims))
    # push!(test, LPDMTest(:blind_vl, :discrete, reward_mode, pconfig, n_sims))

    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.0, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.1, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.2, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.3, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.4, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.5, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.6, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.7, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.8, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 0.9, pconfig, n_sims))
    push!(test, LPDMTest(:adaptive, :discrete, reward_mode, 1.0, pconfig, n_sims))

    # CONTINUOUS OBSERVATIONS
    # push!(test, LPDMTest(:standard, :continuous, reward_mode, 0.0, pconfig, n_sims))
    # push!(test, LPDMTest(:extended, :continuous, reward_mode, pconfig, n_sims))
    # push!(test, LPDMTest(:blind_vl, :continuous, reward_mode, pconfig, n_sims))
    # push!(test, LPDMTest(:adaptive, :continuous, reward_mode, pconfig, n_sims))

    scen=Array{LPDMScenario{S}}(undef,0)

    if dims == 1
        # push!(scen, LPDMScenario(LD1State(-2*π)))
        push!(scen, LPDMScenario(LD1State(π/2)))
        # push!(scen, LPDMScenario(LD1State(3/2*π)))
        # push!(scen, LPDMScenario(LD1State(2*π)))
    elseif dims == 2
        # push!(scen, LPDMScenario(LD2State(-2*π, π)))
        push!(scen, LPDMScenario(LD2State(π/2, -π/2)))
        # push!(scen, LPDMScenario(LD2State(π, 2*π)))
        # push!(scen, LPDMScenario(LD2State(2*π, -π)))
    end

    # Dummy execution, just to make sure all the code is compiled and loaded,
    # to improve uniformity of subsequent executions.
    # execute(solv_mode = :despot, action_space_type = :small, n_sims = 1, s0 = LD2State(π,π), debug = 0)
    # execute(solv_mode = :lpdm, action_space_type = :adaptive, n_sims = 1, s0 = LD2State(π,π), debug = 0)

    f = open("results_" * Dates.format(now(),"yyyy-mm-dd_HH_MM") * ".txt", "w")

    # Printf.@printf(f,"GENERAL SOLVER PARAMETERS\n")
    # Printf.@printf(f,"\tsteps:\t\t\t%d\n", steps)
    # Printf.@printf(f,"\ttime per move:\t\t\t%f\n", time_per_move)
    # Printf.@printf(f,"\tsearch depth:\t\t\t%d\n", search_depth)
    # Printf.@printf(f,"\tN particles:\t\t\t%d\n", n_particles)
    # Printf.@printf(f,"\tmax trials:\t\t\t%d\n\n", max_trials)
    print(f,sconfig)
    Printf.@printf(f,"\nPROBLEM PARAMETERS\n")
    Printf.@printf(f,"\tdimensions:\t\t\t%d\n",             dims)
    Printf.@printf(f,"\tN bins (per dim):\t\t\t%d\n",       pconfig.n_bins)
    Printf.@printf(f,"\tmax actions:\t\t\t%d\n",            pconfig.max_actions)
    Printf.@printf(f,"\tN new actions:\t\t\t%d\n",          pconfig.n_new_actions)
    Printf.@printf(f,"\taction range fraction:\t\t\t%f\n",  pconfig.action_range_fraction)
    Printf.@printf(f,"\tmax exploit. visits:\t\t\t%d\n",    pconfig.max_exploit_visits)
    Printf.@printf(f,"\tmax belief clusters:\t\t\t%d\n",    pconfig.max_belief_clusters)
    Printf.@printf(f,"\ttests per scenario:\t\t\t%d\n\n",   n_sims)

    for i in 1:length(scen)
        if debug >= 0
            println("")
            println("SCENARIO $i, s0 = $(scen[i].s0)")
            println("------------------------")
        end

        Printf.@printf(f,"SCENARIO %d, s0 = %s\n", i, "$(scen[i].s0)")
        Printf.@printf(f,"==========================================================================\n")
        Printf.@printf(f,"ACT. MODE\t\tOBS. MODE\t\t\tACT. STD\t\t\tSTEPS (STD)\t\t\t\tREWARD (STD)\n")
        Printf.@printf(f,"==========================================================================\n")
        for t in test
            if debug >= 0
                println("dimensions: $dims, actions: $(t.action_mode), observations: $(t.obs_mode), action std: $(t.action_std), rewards: $(t.reward_mode)")
            end
            steps_avg, steps_std, reward_avg, reward_std =
                        run_scenario(scen[i].s0, t, A, O, B, sconfig,
                                     dims      = dims,
                                     n_sims    = n_sims,
                                     vis       = vis)

            Printf.@printf(f,"%s\t\t%s\t\t\t%0.2f\t\t\t\t\t%05.2f (%06.2f)\t\t%06.2f (%06.2f)\n",
                        string(t.action_mode), string(t.obs_mode), t.action_std, steps_avg, steps_std, reward_avg, reward_std)
            flush(f)
            debug >=0 && println("\n*******************************************************************************************************************************")
            debug >=0 && println("STATS: actions = $(t.action_mode), observations = $(t.obs_mode), action std = $(t.action_std), steps = $(steps_avg) ($steps_std), reward = $(reward_avg) ($reward_std)")
            debug >=0 && println("*******************************************************************************************************************************\n")
        end
        Printf.@printf(f,"==========================================================================\n\n")
    end
    close(f)
end

function run_scenario(s0::S, test::LPDMTest, A::Type, O::Type, B::Type, sconfig::LPDMConfig;
                      dims::Int64                 = 1,
                      n_sims::Int64               = 1,
                      vis::Vector{Int64}          = Int64[]
                      ) where {S}

    if dims == 1
        p = LightDark1DLpdm(action_mode             = test.action_mode,
                            obs_mode                = test.obs_mode,
                            reward_mode             = test.reward_mode,
                            action_std              = test.action_std,
                            n_bins                  = test.problem_config.n_bins,
                            max_actions             = test.problem_config.max_actions,
                            n_new_actions           = test.problem_config.n_new_actions,
                            action_range_fraction   = test.problem_config.action_range_fraction,
                            max_exploit_visits      = test.problem_config.max_exploit_visits,
                            max_belief_clusters     = test.problem_config.max_belief_clusters
                            )
    elseif dims == 2
        p = LightDark2DLpdm(action_mode             = test.action_mode,
                            obs_mode                = test.obs_mode,
                            reward_mode             = test.reward_mode,
                            action_std              = test.action_std,
                            n_bins                  = test.problem_config.n_bins,
                            max_actions             = test.problem_config.max_actions,
                            n_new_actions           = test.problem_config.n_new_actions,
                            action_range_fraction   = test.problem_config.action_range_fraction,
                            max_exploit_visits      = test.problem_config.max_exploit_visits,
                            max_belief_clusters     = test.problem_config.max_belief_clusters
                            )
    end

    sim_rewards = Vector{Float64}(undef,n_sims)
    sim_steps   = Vector{Int64}(undef,n_sims)

    for sim in 1:n_sims

        world_rng = RNGVector(1, UInt32(sim))
        LPDM.set!(world_rng, 1)
        step_rewards::Array{Float64}     = Vector{Float64}(undef,0)

        sconfig.seed        = UInt32(2*sim+1)
        sconfig.action_mode = test.action_mode
        sconfig.obs_mode    = test.obs_mode
        solver = LPDM.LPDMSolver{S, A, O, B, RNGVector}(config = sconfig)

    #---------------------------------------------------------------------------------
        # Belief
        bu = LPDMBeliefUpdater(p,
                               n_particles = sconfig.n_particles,
                               seed = UInt32(3*sim+1));  # initialize belief updater
        initial_states = state_distribution(p, s0, sconfig, world_rng)     # create initial  distribution
        current_belief = LPDM.create_belief(bu)                       # allocate an updated belief object
        LPDM.initialize_belief(bu, initial_states, current_belief)    # initialize belief
        updated_belief = LPDM.create_belief(bu)

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

        policy::LPDMPolicy = POMDPs.solve(solver, p)
        s = s0
        step::Int64 = 0
        r::Float64 = 0.0

        if sconfig.debug >= 0
            println("")
            println("*** SIM $sim (dimensions: $dims, action mode: $(test.action_mode), obs mode: $(test.obs_mode), action std: $(test.action_std), s0: $s0)***")
            println("")
        end

        val, run_time, bytes, gctime, memallocs =
        @timed while !isterminal(p, s) && (sconfig.sim_len == -1 || step < sconfig.sim_len)
            step += 1

            if sconfig.debug >= 1
                println("")
                println("=============== Step $step ================")
                show(current_belief)
            end

            a = POMDPs.action(policy, current_belief)
            if sconfig.debug >= 1
                println("s: $s")
                println("a: $a")
            end

            # for now just use the deterministic action outcome version for the "world"
            s, o, r = POMDPs.generate_sor(p, s, a, world_rng, true)

            push!(step_rewards, r)
            if sconfig.debug >= 1
                println("s': $s")
                println("o': $o")
                println("r: $r")
                println("=======================================")
            end

            # update belief
            POMDPs.update(bu, current_belief, a, o, updated_belief)
            current_belief = deepcopy(updated_belief)
            # show(updated_belief) #NOTE: don't show for now

            if LPDM.isterminal(p, current_belief.particles)
                if sconfig.debug >= 1
                    println("Terminal belief. Execution completed.")
                    show(current_belief)
                end
                break
            end

            if step ∈ vis
                t = LPDM.d3tree(solver,
                                detect_repeat=false,
                                title="Step $step",
                                init_expand=2)
                # # show(t)
                inchrome(t)
                # blink(t)
            end
            if sconfig.debug >= 1
                println("root actions: $(solver.root.action_space)")
            end
        end

        sim_steps[sim] = step
        sim_rewards[sim] = sum(step_rewards)
        sconfig.debug >= 0 && println("steps=$(sim_steps[sim]), reward=$(sim_rewards[sim])")
    end

    return mean(sim_steps), std(sim_steps), mean(sim_rewards), std(sim_rewards)

end

end #module

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
"""
A SemanticWrapper determines the semantics of data that it wraps.
Any subtype needs to implement function `unwrap(wrapper)` that returns the
wrapped data.
"""
abstract type SemanticWrapper end

"""
Mask wraps a segmentation mask.
"""
struct Mask{AT<:AbstractArray} <: SemanticWrapper
    img::AT
end

"""
    unwrap(sw::SemanticWrapper)

Returns the original object.
"""
unwrap(m::Mask) = m.img

"""
    shouldapply(op, wrapper)
    shouldapply(typeof(op), typeof(wrapper))

Determines if operation `op` should be applied to semantic wrapper `wrapper`.
"""
shouldapply(op::Operation, what::SemanticWrapper) = shouldapply(typeof(op), typeof(what))
shouldapply(::Type{<:ImageOperation}, ::Type{<:SemanticWrapper}) = Val(true)
shouldapply(::Type{<:ColorOperation}, ::Type{<:Mask})            = Val(false)
# By default any operation is applicable to any semantic wrapper. Add new
# methods to this function to define exceptions.

# Allows doing `unwrap.(augment(img, Mask(img2), pl))`
unwrap(A::AbstractArray) = A


text:
using PowerModelsSecurityConstrained
using PowerModels

using Memento

using Test

using Ipopt
using Cbc

Memento.setlevel!(Memento.getlogger(PowerModels), "error")
Memento.setlevel!(Memento.getlogger(PowerModelsSecurityConstrained), "error")

nlp_solver = with_optimizer(Ipopt.Optimizer, tol=1e-6, print_level=0)
#nlp_solver = with_optimizer(Ipopt.Optimizer, tol=1e-6)
lp_solver = with_optimizer(Cbc.Optimizer, logLevel=0)
#lp_solver = with_optimizer(Cbc.Optimizer)

ini_file = "../test/data/inputfiles.ini"
case_01 = parse_goc_files(ini_file, scenario_id="scenario_01")
case_02 = parse_goc_files(ini_file, scenario_id="scenario_02")
cases = [case_01, case_02]

networks = [build_pm_model(case) for case in cases]
solutions = [read_solution1(networks[i], output_dir=dirname(case.files["raw"])) for (i,case) in enumerate(cases)]

case_infeasible = parse_goc_files(ini_file, scenario_id="scenario_03")
network_infeasible = build_pm_model(case_infeasible)

@testset "PowerModelsSecurityConstrained" begin

    include("common.jl")

    include("opf.jl")

    include("contingency.jl")

    include("contingency-stage.jl")

    include("scopf.jl")

end


text:
module DotfileManager

    include("argparse.jl")

end # module


text:
using DataStructures

@testgroup "fixtures" begin

@testcase "testcase with run_options" for ro in run_options
    @test haskey(ro, :verbosity)
end

fx_with_run_options = @global_fixture for ro in run_options
    @produce ro[:verbosity]
end

@testcase "testcase with a fixture with run_options" for v in fx_with_run_options
    @test isa(v, Int)
end

@testcase "constant fixture with non-1D iterable" begin
    fx = Jute.ConstantFixture(1)
    @test Jute.setup(fx) == [Jute.LabeledValue(1)]

    fx = Jute.ConstantFixture([1 2; 3 4])
    @test Jute.setup(fx) == map(Jute.LabeledValue, [1, 3, 2, 4])

    fx = Jute.ConstantFixture([1 2; 3 4], ["one" "two"; "three" "four"])
    @test Jute.setup(fx) == map(Jute.LabeledValue, [1, 3, 2, 4], ["one", "three", "two", "four"])
end

@testcase "constant fixture with labels" begin
    fx = Jute.ConstantFixture([1], ["one"])
    @test Jute.setup(fx) == [Jute.LabeledValue(1, "one")]
end

@testcase "constant fixture checks for iterable" begin
    @test_throws ErrorException Jute.ConstantFixture(:a)
    @test_throws ErrorException Jute.ConstantFixture([1], :a)
end

@testcase "global fixture checks for callable" begin
    @test_throws ErrorException Jute.GlobalFixture(1)
end

@testcase "local fixture checks for callable" begin
    @test_throws ErrorException Jute.LocalFixture(1)
end

@testcase "constant fixture from a pair" begin
    fx = @global_fixture for x in ([1, 2] => ["one", "two"])
        @produce x
    end
    constant_fx = Jute.parameters(fx)[1]
    @test Jute.setup(constant_fx) == map(Jute.LabeledValue, [1, 2], ["one", "two"])
end

testcases_for_custom_label_test = Jute.collect_testobjs() do
    fx_with_custom_label = @global_fixture for x in [1, 2]
        # Try a string and non-string label as a regression test
        # (@produce used to hang if passed a non-string label).
        label = x == 1 ? "one" : 2
        @produce x label
    end

    @testcase "custom label" for x in fx_with_custom_label
        @test 1 == 1
    end
end

@testcase "produce with a custom label" begin
    exitcode, output = nested_run_with_output(
        testcases_for_custom_label_test, Dict(:verbosity => 2))
    @test exitcode == 0

    template = """
        Collecting testcases...
        Using 1 out of 1 testcase definitions...
        ================================================================================
        Platform: <<<platform>>>, Julia <<<julia_version>>>, Jute <<<jute_version>>>
        --------------------------------------------------------------------------------
        custom label [one] (<<<time>>>) [PASS]
        custom label [2] (<<<time>>>) [PASS]
        --------------------------------------------------------------------------------
        2 tests passed, 0 failed, 0 errored in <<<full_time>>> (total test time <<<test_time>>>)
    """

    test_match_text(template, output)
end

@testcase "fixture dependencies" begin
    fx1 = @global_fixture begin
        @produce 1
    end

    fx2 = @global_fixture for x in fx1
        @produce x
    end

    fx3 = @global_fixture for x in fx1
        @produce x
    end

    fx4 = @local_fixture for x in fx3, y in fx2
        @produce x + y
    end

    fx5 = @global_fixture for x in fx3, y in fx2
        @produce x + y
    end

    @test Jute.dependencies(fx4) == OrderedSet([fx1, fx3, fx2])
    @test Jute.parameters(fx4) == [fx3, fx2]

    @test Jute.dependencies(fx5) == OrderedSet([fx1, fx3, fx2])
    @test Jute.parameters(fx5) == [fx3, fx2]
end

end


text:
using Flux
using Random
using MinimalRLCore

"""
    DQNAgent

A DQN which implements the MinimalRLCore interface with support for managed RNGs.
Currently only supports a vanilla Experience Replay buffer (as implemented in ExperienceReplay) with
plans to extend to more general buffers with dispatch.
Plans to include all changes made by Rainbow taking advantage of dispatch.

# Arguments:
- `model`: A model which is compatable with the learning update and optimizer. Will call update! (see update_parameters).
- `target_network`: A target network, which can be None if using default QLearning updates or if you have specialized update!.
- `learning_update`: A learning update (ala QLearning or DoubleQLearning).
- `optimizer`: An optimizer compatable w/ learning update and the model type.
- `acting_policy::AbstractValuePolicy`: An acting policy. Used to determine actions.
- `replay_size`: The size of the replay buffer.
- `hist_length`: The length of the state history used for the input state to the network.
- `example_state`: An example of the state produced by the environment which is sent to the Agent.
- `batch_size`: The batch_size used for updating.
- `target_update_freq`: The number of agent steps between target network updates.
- `update_freq`: The number of agent steps between updating the model.
- `min_mem_size`: The size of the replay buffer before learning begins.

# Kwargs:
- `device = Val{:cpu}()`: The device the agent uses (either `Val{:cpu}()` or `Val{:gpu}()`).
- `state_preproc = identity`: A function for processing observations before storage in the state_buffer
- `state_postproc = identity`: A function for processing the observations after storage in the state_buffer.
- `hist_squeeze=false`: Whether to squeeze the observations onto a single dimension. (see `HistStateBuffer`). Can be either a boolean or Val{bool}
"""
Base.@kwdef mutable struct DQNAgent{ER<:AbstractReplay, SB, AP<:AbstractValuePolicy, LU, M, TN, O, Φ, SP1, SP2, RT, UC<:Val} <: AbstractAgent
    # Models
    model::M # Currently assumed to be a Flux model
    target_network::TN # Can be either Nothing or a Flux model

    # Learning
    learning_update::LU # Follows the implementation of a Learning Update.
    optimizer::O # Flux optimizer, but could wrap your own (see optimizers.jl)

    # Acting policy (Abstract Value Policy)
    acting_policy::AP # An AbstractValuePolicy. This can be ϵGreedy or something else.

    # State/experience processing.
    replay::ER # Experience Replay. Currently, we only support a vanilla replay but will change in the future based on dispatch.
    state_buffer::SB # A state buffer. Can be Nothing, but the current constructor uses StateBuffer of HistStateBuffer based on hist_length.
    state_preproc::SP1 # A function for processing observations before storage in the state_buffer
    state_postproc::SP2 # Another function for processing the observations after storage in the state_buffer.
    rew_transform::RT # A function for processing the rewards.
    prev_s::Φ # storing the prev_state for adding to the er buffer.

    # params
    batch_size::Int
    target_update_freq::Int
    update_freq::Int
    min_mem_size::Int

    # minor details
    action::Int = 0
    training_steps::Int = 0

    # extra
    device::UC = Val{:cpu}() # Determines where to send the data when doing updates. I feel like we should also send the model and target_network but currently don't.

end

function DQNAgent(model,
                  target_network,
                  learning_update,
                  optimizer,
                  acting_policy,
                  replay_size,
                  hist_length,
                  example_state,
                  batch_size,
                  target_update_freq,
                  update_freq,
                  min_mem_size;
                  device = Val{:cpu}(),
                  state_preproc = identity,
                  state_postproc = identity,
                  rew_transform = identity,
                  hist_squeeze=false)

    proc_state = state_preproc(example_state)
    @assert hist_length >= 1
    state_buffer = if hist_length == 1
        DeepRL.StateBuffer{eltype(proc_state)}(replay_size, length(proc_state))
    else
        DeepRL.HistStateBuffer{eltype(proc_state)}(replay_size, size(proc_state), hist_length, hist_squeeze)
    end

    prev_s = if state_buffer isa Nothing
        state_processor(example_state)
    elseif state_buffer isa StateBuffer
        0
    elseif state_buffer isa HistStateBuffer
        zeros(Int, state_buffer.hist_length)
    else
        throw("Unknown StateBuffer please use default DQN constructor.")
    end

    replay = ExperienceReplayDef(replay_size, length(prev_s), eltype(prev_s))

    DQNAgent(model = model,
             target_network = target_network,
             learning_update = learning_update,
             optimizer = optimizer,
             acting_policy = acting_policy,
             replay = replay,
             state_buffer = state_buffer,
             state_preproc = state_preproc,
             state_postproc = state_postproc,
             rew_transform = rew_transform,
             prev_s = prev_s,
             batch_size = batch_size,
             target_update_freq = target_update_freq,
             update_freq = update_freq,
             min_mem_size = min_mem_size,
             device = device)

end

function process_state(agent::DQNAgent, s; start=false)
    # this stores the state in the state buffer (if not Nothing) and preproccesses the state.
    if agent.state_buffer isa Nothing
        agent.state_preproc(s)
    else
        Base.push!(agent.state_buffer, agent.state_preproc(s); new_episode=start)
        laststate(agent.state_buffer)
    end
end

function process_input(agent::DQNAgent, s, r; start=false)
    proc_state = process_state(agent, s; start)
    proc_rew = agent.rew_transform(r) # Atari implementation returns float64s for the reward.
    proc_action = findfirst((a)->a==agent.action,
                            action_set(agent.acting_policy))::Int
    proc_state, proc_action, proc_rew
end

function get_action(agent::DQNAgent, s, rng)
    if agent.acting_policy isa ϵGreedyDecay
        sample(agent.acting_policy,
               to_host(agent.model(s)),
               agent.training_steps,
               rng)
    else
        sample(agent.acting_policy,
               to_host(agent.model(s)),
               rng)
    end
end

function get_state_from_buffer(agent::DQNAgent, s)
    if agent.state_buffer isa Nothing
        to_device(agent.device, agent.state_postproc(s))
    else
        to_device(agent.device, agent.state_postproc(agent.state_buffer[s]))
    end
end

function MinimalRLCore.start!(agent::DQNAgent,
                              env_s_tp1,
                              rng::AbstractRNG=Random.GLOBAL_RNG)

    agent.prev_s = process_state(agent, env_s_tp1; start=true)
    state = get_state_from_buffer(agent, agent.prev_s)
    resh_state = reshape(state, size(state)..., 1)

    agent.action = get_action(agent, resh_state, rng)

    return agent.action
end

function MinimalRLCore.step!(agent::DQNAgent,
                             env_s_tp1,
                             r,
                             terminal,
                             rng::AbstractRNG=Random.GLOBAL_RNG)

    # proc_state = process_state(agent, env_s_tp1)

    proc_state, a_t, proc_rew = process_input(agent, env_s_tp1, r)

    add_ret = Base.push!(agent.replay,
                         (agent.prev_s,
                          a_t,
                          proc_state,
                          proc_rew, # Atari implementation returns float64s for the reward.
                          terminal))

    update_params!(agent, rng)

    agent.prev_s = proc_state
    prev_s = get_state_from_buffer(agent, agent.prev_s)
    resh_prev_s = reshape(prev_s, size(prev_s)..., 1)

    agent.action = get_action(agent, resh_prev_s, rng)

    return agent.action
end

function update_params!(agent::DQNAgent, rng)

    us = if length(agent.replay) > agent.min_mem_size &&
        agent.training_steps % agent.update_freq == 0

        e = sample(rng,
                   agent.replay,
                   agent.batch_size)

        s = get_state_from_buffer(agent, e.s)
        r = to_device(agent.device, e.r)
        t = to_device(agent.device, e.t)
        sp = get_state_from_buffer(agent, e.sp)

        update!(agent.model,
                agent.learning_update,
                agent.optimizer,
                s,
                e.a,
                sp,
                r,
                t,
                agent.target_network)

    end

    # Target network updates
    if !(agent.target_network isa Nothing)
        if agent.training_steps%agent.target_update_freq == 0
            update_target_network(agent.model, agent.target_network)
        end
    end

    agent.training_steps += 1

    us
end


text:
struct Hexagonal end

capacity(::Hexagonal, d::Int) = 3 * d

@inline function eachdir(::Hexagonal, x::Site{2})
  return ((x[1] - 1, x[2]),
          (x[1] + 1, x[2]),
          (x[1], x[2] - 1),
          (x[1], x[2] + 1),
          (x[1] + 1, x[2] - 1),
          (x[1] - 1, x[2] + 1))
end

function distance(::Hexagonal, x::Site{2}, y::Site{2})
  d1 = x[1] - y[1]
  d2 = x[2] - y[2]

  return sign(d1) == sign(d2) ? abs(d1 + d2) : max(abs(d1), abs(d2))
end


text:
const FILTERED_ACTIVE_CMDS = keys((filter(check -> last(check) == true, ACTIVE_COMMANDS)))

function commander(c::Client, m::Message, ::Val{:source})
    startswith(m.content, COMMAND_PREFIX * "src") || return nothing

    regex = Regex(COMMAND_PREFIX * raw"src( help| [a-zA-Z]*)? *(.*)$")

    matches = match(regex, m.content)

    if matches === nothing || matches.captures[1] ∈ (" help", nothing)
        help_commander(c, m, :source)
    else
        command = lowercase(strip(matches.captures[1]))
        if Symbol(command) ∈ FILTERED_ACTIVE_CMDS
            msg = string(BOT_REPO_URL, "/blob/main/src/command/", command, ".jl")
            reply(c, m, msg)
        else
            reply(c, m, "No such blob or file. Please check out $(BOT_REPO_URL) instead")
        end
    end
    return nothing
end

function help_commander(c::Client, m::Message, ::Val{:source})
    list_commands = join(
        sort([string("src", " ", cmd, "\n") for cmd in FILTERED_ACTIVE_CMDS])
    )

    return reply(
        c,
        m,
        """
        Returns the source file of the command.
        Usage: `src <command>`
        List of available commands:
        ```
        $(list_commands)
        ```
        `src` or `src help` returns help
        """,
    )
end


text:
# # Using a custom model data structure

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/@__NAME__.ipynb)

# This notebooks shows how to utilize the generic accessors and modification
# functions in COBREXA.jl to run the analysis on any custom model type. We will
# create a simple dictionary-style structure that describes the model, allow
# COBREXA to run a FVA on it, and create a simple reaction-removing
# modification.

# First, let's define a very simple stoichiometry-only structure for the model:

using COBREXA

mutable struct MyReaction
    max_rate::Float64 # maximum absolute conversion rate
    stoi::Dict{String,Float64} # stoichimetry of the reaction

    MyReaction() = new(0.0, Dict{String,Float64}())
end

mutable struct MyModel <: MetabolicModel
    optimization_target::String # the "objective" reaction name
    reactions::Dict{String,MyReaction} # dictionary of reactions

    MyModel() = new("", Dict{String,MyReaction}())
    MyModel(o, r) = new(o, r)
end

# With this, we can start defining the accessors:

COBREXA.n_reactions(m::MyModel) = length(m.reactions)
COBREXA.reactions(m::MyModel) = sort(collect(keys(m.reactions)))

# Metabolites are defined only very implicitly, so let's just make a function
# that collects all names. [`n_metabolites`](@ref) can be left at the default
# definition that just measures the output of [`metabolites`](@ref).

function COBREXA.metabolites(m::MyModel)
    mets = Set{String}()
    for (_, r) in m.reactions
        for (m, _) in r.stoi
            push!(mets, m)
        end
    end
    return sort(collect(mets))
end

# Now, the extraction of the linear model. Remember the order of element in the
# vectors must match the order in the output of [`reactions`](@ref) and
# [`metabolites`](@ref).

using SparseArrays

function COBREXA.bounds(m::MyModel)
    max_rates = [m.reactions[r].max_rate for r in reactions(m)]
    (sparse(-max_rates), sparse(max_rates))
end

function COBREXA.objective(m::MyModel)
    if m.optimization_target in keys(m.reactions)
        c = spzeros(n_reactions(m))
        c[first(indexin([m.optimization_target], reactions(m)))] = 1.0
        c
    else
        throw(
            DomainError(
                m.optimization_target,
                "The target reaction for flux optimization not found",
            ),
        )
    end
end

function COBREXA.stoichiometry(m::MyModel)
    sparse([
        get(m.reactions[rxn].stoi, met, 0.0) for met in metabolites(m), rxn in reactions(m)
    ])
end

# Now the model is complete! We can generate a random one and see how it
# performs
import Random
Random.seed!(123)

rxn_names = ["Reaction $i" for i = 'A':'Z'];
metabolite_names = ["Metabolite $i" for i = 1:20];

m = MyModel();
for i in rxn_names
    m.reactions[i] = MyReaction()
end

for i = 1:50
    rxn = rand(rxn_names)
    met = rand(metabolite_names)
    m.reactions[rxn].stoi[met] = rand([-3, -2, -1, 1, 2, 3])
    m.reactions[rxn].max_rate = rand()
end

# Let's see what the model looks like now:
m

# We can run most of the standard function on the model data right away:
using Tulip
m.optimization_target = "Reaction A"
flux_balance_analysis_dict(m, Tulip.Optimizer)

# To be able to use the model conveniently in functions such as
# [`screen`](@ref), you usually want to be able to easily specify the
# modifications. In this example, we enable use of
# [`with_removed_reactions`](@ref) by overloading the internal
# [`remove_reactions`](@ref) for this specific model type:
#
# We need to make an as-shallow-as-possible copy of the model that allows us to
# remove the reactions without breaking the original model.

function COBREXA.remove_reactions(m::MyModel, rxns::AbstractVector{String})
    m = MyModel(m.optimization_target, copy(m.reactions))
    delete!.(Ref(m.reactions), rxns)
    return m
end

# The screening is ready now!

reactions_to_remove = ("Reaction $i" for i = 'B':'Z')

reactions_to_remove .=> screen_variants(
    m,
    [[with_removed_reactions([r])] for r in reactions_to_remove],
    m -> flux_balance_analysis_dict(m, Tulip.Optimizer)["Reaction A"],
)


text:
import Obj
import Ls

function bt()
    println("\nTesting Ls.bt...")
    maxiter = 100

    # init intractable
    obj = Obj.Objective(
        w -> w[1]^2,
        w -> [2w[1], 0]
    )
    bt = Ls.bt(obj, maxiter=maxiter)
    w = [0.0, 0.0]
    s = [1.0, 0.0]
    # assert
    α, iter1 = Ls.ls(bt, w, s)
    success1 = (iter1 == maxiter)

    # init tractable
    obj = Obj.Objective(
        w -> -w[1]^2,
        w -> [-2w[1], 0]
    )
    bt = Ls.bt(obj, maxiter=maxiter)
    w = [0.0, 0.0]
    s = [1.0, 0.0]
    # assert
    α, iter2 = Ls.ls(bt, w, s)
    success2 = (iter2 < maxiter)

    # print results
    if (success1 && success2)
        println("SUCCEEDED: Intractable LS required max number of iterations, tractable LS required less.")
    elseif (!success1)
        println("FAILED: Intractable LS required less than max number of iterations ($iter1)!")
    elseif (!success2)
        println("FAILED: Tractable LS required max number of iterations!")
    end

    return nothing
end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso