text:
"""
    fill_vertices!(pvector, parray)
"""
function fill_vertices!(pvector, parray)
    dim = size(parray[1], 1)
    # check dimensions
    if sum(length.(parray)) != length(pvector)
        error("pb of dimensions in fill_vertices $(sum(length.(parray))) != $(length(pvector)) ")
    end
    c = 1
    for k = 1:length(parray)
        for l = 1:dim # fill by column
            for m = 1:size(parray[k], 1)
                parray[k][m, l] = pvector[c]
                c += 1
            end
        end
    end
    return nothing
end
#-------------------------------------------------------------------------------
function costnormU(pvin, cellsb, t, rhs, boundary_condition, meshboundary)
    # copy coordinates
    pv = reshape(pvin, length(pvin) ÷ 2, 2)
    n_dofs, n_polys = size(pv, 1), 3 # method has 1 degree of freedom per vertex
    u = zeros(n_dofs) # degrees of freedom of the virtual element solution
    # call assemble function
    IK, JK, SK, IF, SF = assembKM_vemKsource(pv, cellsb, rhs,
                                             boundary_condition, meshboundary)
    K = sparse(IK, JK, SK, n_dofs, n_dofs)
    F = Vector(sparsevec(IF, SF, n_dofs))
    boundary_vals = boundary_condition(pv[meshboundary, :])
    # vertices which aren’t on the boundary
    internal_dofs =  setdiff(1:n_dofs, meshboundary)
    # solve
    u[internal_dofs] = K[internal_dofs, internal_dofs] \ F[internal_dofs]
    u[meshboundary] .= boundary_vals # set the boundary values
    # L² cost
    val = sum(u[k] ^ 2 for k = 1:n_dofs) / n_dofs
    return val
end
#-------------------------------------------------------------------------------
function ∇costnormU(pvin, cellsb, t, rhs, boundary_condition, meshboundary)
    # copy coordinates
    pv = reshape(pvin, length(pvin) ÷ 2, 2)
    n_dofs, n_polys = size(pv, 1), 3 # method has 1 degree of freedom per vertex
    u = zeros(n_dofs) # degrees of freedom of the virtual element solution
    # call assemble function
    IK, JK, SK, IF, SF = assembKM_vemKsource(pv, cellsb, rhs,
                                             boundary_condition, meshboundary)
    K = sparse(IK, JK, SK, n_dofs, n_dofs)
    F = Vector(sparsevec(IF, SF, n_dofs))
    boundary_vals = boundary_condition(pv[meshboundary, :])
    # vertices which aren’t on the boundary
    internal_dofs =  setdiff(1:n_dofs, meshboundary)
    # solve
    u[internal_dofs] = K[internal_dofs, internal_dofs] \ F[internal_dofs]
    u[meshboundary] .= boundary_vals # set the boundary values
    # cost gradient and adjoint function
    JFU = - 2 * u / n_dofs
    λ = K' \ JFU
    λ = zeros(length(u))
    λ[internal_dofs] = K'[internal_dofs, internal_dofs] \ JFU[internal_dofs]
    # compute full gradient by automatic differentiation
    fprod(x) = prodgrad(x, cellsb, t, rhs, boundary_condition, meshboundary,
                        u, λ)
    ∇prodgrad = ReverseDiff.gradient(fprod, pv)
    #= optimization tests
    ∇prodgrad2 = similar(∇prodgrad)
    f_tape = GradientTape(fprod, randn(size(pv)))
    compiled_f_tape = compile(f_tape)
    @time gradient!(∇prodgrad2, compiled_f_tape, pv)
    println("diff grads = $(norm(∇prodgrad2 - ∇prodgrad2))")
    =#
    # gradient contribution of boundary points
    fboundary(x) = sum(boundary_condition(x[meshboundary, :]) .^ 2) / n_dofs
    ∇Dirichlet = ReverseDiff.gradient(fboundary, pv)
    return ∇prodgrad + ∇Dirichlet
end
#-------------------------------------------------------------------------------
function prodgrad(pv, cellsb, t, rhs, boundary_condition, meshboundary, u, λ)
    # assemble matrices
    IK, JK, SK, IF, SF = assembKM_vemKsource(pv, cellsb, rhs,
                                             boundary_condition, meshboundary)
    val = 0.
    # consider only internal points
    Jinmesh = indexin(JK, meshboundary)
    Jboundary = findall(Jinmesh .!= nothing)
    SK[Jboundary] .= 0.
    for (i, j, s) ∈ zip(IK, JK, SK)
        val += λ[i] * s * u[j]
    end
    for (i, s) ∈ zip(IF, SF)
        val -= λ[i] * s
    end
    return val
end
#-------------------------------------------------------------------------------
function costnormU_elas2D(pvin, cellsb, t, rhs, boundary_condition,
                          meshboundary)
    # copy coordinates
    pv = reshape(pvin, length(pvin) ÷ 2, 2)
    # call assemble function
    IK, JK, SK, IF, SF = assembKM_vemKsource_elas2D(pv, cellsb, rhs,
                                                     boundary_condition,
                                                     meshboundary)
    # solve
    K, F, internal_dofs, u = solve_elas2D(IK, JK, SK, IF, SF, pv, meshboundary,
                                          boundary_condition)
    # L² cost
    n_dofs = length(u)
    val = sum(u[k] ^ 2 for k = 1:n_dofs) / n_dofs
    return val
end
#-------------------------------------------------------------------------------
function ∇costnormU_elas2D(pvin, cellsb, t, rhs, boundary_condition,
                           meshboundary)
    # copy coordinates
    pv = reshape(pvin, length(pvin) ÷ 2, 2)
    # call assemble function
    IK, JK, SK, IF, SF = assembKM_vemKsource_elas2D(pv, cellsb, rhs,
                                             boundary_condition, meshboundary)
    # solve
    K, F, internal_dofs, u = solve_elas2D(IK, JK, SK, IF, SF, pv, meshboundary,
                                          boundary_condition)
    n_dofs = length(u)
    # cost gradient and adjoint function
    JFU = - 2 * u / n_dofs
    λ = K' \ JFU
    λ = zeros(n_dofs)
    λ[internal_dofs] = K'[internal_dofs, internal_dofs] \ JFU[internal_dofs]
    # compute full gradient by automatic differentiation
    fprod(x) = prodgrad_elas2D(x, cellsb, t, rhs, boundary_condition,
                               meshboundary, u, λ)
    ∇prodgrad = ReverseDiff.gradient(fprod, pv)
    # gradient contribution of boundary points
    fboundary(x) = sum(boundary_condition(x[meshboundary, :]) .^ 2) / n_dofs
    ∇Dirichlet = ReverseDiff.gradient(fboundary, pv)
    return ∇prodgrad #+ ∇Dirichlet
end
#-------------------------------------------------------------------------------
function prodgrad_elas2D(pv, cellsb, t, rhs, boundary_condition,
                         meshboundary, u, λ)
    # assemble matrices
    IK, JK, SK, IF, SF = assembKM_vemKsource_elas2D(pv, cellsb, rhs,
                                                    boundary_condition,
                                                    meshboundary)
    val = 0.
    # boundary dofs
    #=
    Vmeshboundary = vcat(2 * meshboundary .- 1, 2 * meshboundary)
    Jinmesh = indexin(JK, Vmeshboundary)
    Jboundary = findall(Jinmesh .!= nothing)
    SK[Jboundary] .= 0.
    =#
    for (i, j, s) ∈ zip(IK, JK, SK)
        val += λ[i] * s * u[j]
    end
    for (i, s) ∈ zip(IF, SF)
        val -= λ[i] * s
    end
    return val
end


text:
using Py2Jl
py2jl"""
def sum_by(f, seq):
    s = 0
    for e in seq:
        s = s + f(e)
    return s
result = sum_by(lambda x: x * 10, [1, 2, 3])
"""
@info :print_result result

@info :interops sum_by(x -> 2x, [1, 2, 3])

@info :transpiled_to_julia py2jl("mymod",
"""
def sum_by(f, seq):
    s = 0
    for e in seq:
        s = s + f(e)
    return s
""")

ext: jl


text:
# distritbutions
@rule typeof(+)(:out, Marginalisation) (m_in1::Distribution, m_in2::Distribution) = begin
    return convolve(m_in1, m_in2)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_var(m_in1)
    min2, vin2 = mean_var(m_in2)
    return NormalMeanVariance(min1 + min2, vin1 + vin2)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::MultivariateNormalDistributionsFamily, m_in2::MultivariateNormalDistributionsFamily) = begin
    min1, vin1 = mean_cov(m_in1)
    min2, vin2 = mean_cov(m_in2)
    return MvNormalMeanCovariance(min1 + min2, vin1 + vin2)
end

# PointMass
@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::PointMass) = PointMass(mean(m_in1) + mean(m_in2))

@rule typeof(+)(:out, Marginalisation) (m_in1::NormalMeanPrecision, m_in2::PointMass) = begin
    min1, win1 = mean_precision(m_in1)
    return NormalMeanPrecision(min1 + mean(m_in2), win1)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::NormalMeanPrecision) = begin
    return @call_rule typeof(+)(:out, Marginalisation) (m_in1 = m_in2, m_in2 = m_in1, meta = meta)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::MvNormalMeanPrecision, m_in2::PointMass) = begin
    min1, win1 = mean_precision(m_in1)
    return MvNormalMeanPrecision(min1 + mean(m_in2), win1)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::MvNormalMeanPrecision) = begin
    return @call_rule typeof(+)(:out, Marginalisation) (m_in1 = m_in2, m_in2 = m_in1, meta = meta)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::UnivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    min1, vin1 = mean_var(m_in1)
    return NormalMeanVariance(min1 + mean(m_in2), vin1)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::UnivariateNormalDistributionsFamily) = begin
    return @call_rule typeof(+)(:out, Marginalisation) (m_in1 = m_in2, m_in2 = m_in1, meta = meta)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::MultivariateNormalDistributionsFamily, m_in2::PointMass) = begin
    min1, vin1 = mean_cov(m_in1)
    return MvNormalMeanCovariance(mean(m_in2) + min1, vin1)
end

@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::MultivariateNormalDistributionsFamily) = begin
    return @call_rule typeof(+)(:out, Marginalisation) (m_in1 = m_in2, m_in2 = m_in1, meta = meta)
end

# specialized
@rule typeof(+)(:out, Marginalisation) (m_in1::MvNormalWeightedMeanPrecision{T1}, m_in2::MvNormalWeightedMeanPrecision{T2}) where { T1 <: LinearAlgebra.BlasFloat, T2 <: LinearAlgebra.BlasFloat } = begin

    min2, vin2 = mean_cov(m_in2)
    vin1 = cov(m_in1)
    T = promote_type(T1, T2)
    BLAS.gemv!('N', one(T), vin1, weightedmean(m_in1), one(T), min2)
    vin2 .+= vin1
    return MvNormalMeanCovariance(min2, vin2)

end

@rule typeof(+)(:out, Marginalisation) (m_in1::MvNormalWeightedMeanPrecision, m_in2::PointMass) = begin

    ξin1, win1 = weightedmean_precision(m_in1)
    ξout = win1*mean(m_in2)
    ξout .+= ξin1
    return MvNormalWeightedMeanPrecision(ξout, win1)

end

@rule typeof(+)(:out, Marginalisation) (m_in1::PointMass, m_in2::MvNormalWeightedMeanPrecision) = begin

    ξin2, win2 = weightedmean_precision(m_in2)
    ξout = win2*mean(m_in1)
    ξout .+= ξin2
    return MvNormalWeightedMeanPrecision(ξout, win2)

end


text:
using BioSequences, GLM, DataFrames, CSV

####################################
## Following Code is adapted from ##
## La Fleur et al. 2022           ##
####################################
# k and BETA for La Fleur dataset
LOGK   = -2.80271176
BETA    = 0.81632623

# Dictionaries
groove_access = Dict(
    "CG" => 43,
    "CA"=>  42, "TG" => 42, #revcomp
    "GG"=>  42, "CC" => 42,  #revcomp
    "GC"=>  25,
    "GA"=>  22, "TC" => 22, #revcomp
    "TA"=>  14,
    "AG"=>  9, "CT" => 9, #revcomp
    "AA"=>  5, "TT" => 5, #revcomp
    "AC"=>  4, "GT" => 4, #revcomp
    "AT"=>  0
)

stacking_dict = Dict(
    "AA" =>-2.41,
    "AC"=>-2.06, "CA" => -2.06, #reverse
    "AG"=>-2.54, "GA" => -2.54, #reverse
    "AT"=>-2.25, "TA" => -2.25, #reverse
    "CC"=>-1.93,
    "CG"=>-2.24, "GC" => -2.24, #reverse
    "CT"=>-2.01, "TC" => -2.01, #reverse
    "GG"=>-2.66,
    "GT"=>-2.46,
    "TG" => -2.46, #reverse
    "TT"=>-1.82
)

average_twist = Dict(
    "AA" => 35.6, "TT" => 35.6, #revcomp
    "AT"=>  32.1, "AT" => 32.1, #revcomp
    "TA"=>  35.3, "TA" => 35.3, #revcomp
    "GG"=>  33.65, "CC" => 33.65, #revcomp
    "GC"=>  40.2,
    "CG"=>  29.9,
    "AG"=>  27.9, "CT" => 27.9, #revcomp
    "GA"=>  36.8, "TC" => 36.8, #revcomp
    "AC"=>  34.0, "GT" => 34.0, #revcomp
    "CA"=>  34.4, "TG" => 34.4 #revcomp
)

# Sequence dependence of DNA bending rigidity, Geggier and Vologodskii
persistence = Dict(
    "AA" => 50.4, "TT" => 50.4, #revcomp
    "AC"=>  55.4, "GT" => 55.4, #revcomp
    "AG"=>  51.0, "CT" => 51.0, #revcomp
    "AT"=>  40.9, "AT" => 40.9, #revcomp
    "CA"=>  46.7, "TG" => 46.7, #revcomp
    "CC"=>  41.7, "GG" => 41.7, #revcomp
    "CG"=>  56.0,
    "GA"=>  54.4, "TC" => 54.4, #revcomp
    "GC"=>  44.6,
    "TA"=>  44.7,
)

RNA_DNA_hybrids = Dict(
    "TT" =>  -1.0, "AA" => -1.0,  #revcomp
    "TG"=>   -2.1, "CA" => -2.1,  #revcomp
    "TC"=>   -1.8, "GA" => -1.8,  #revcomp
    "TA"=>   -0.9,
    "GT"=>   -0.9, "AC" => -0.9,  #revcomp
    "GG"=>   -2.1, "CC" => -2.1,  #revcomp
    "GC"=>   -1.7,
    "GA"=>   -0.9, "TC" => -0.9,  #revcomp
    "CT"=>   -1.3, "AG" => -1.3,  #revcomp
    "CG"=>   -2.7,
    "CC"=>   -2.9, "GG" => -2.9,  #revcomp
    "CA"=>   -1.1, "TG" => -1.1,  #revcomp
    "AT"=>   -0.6,
    "AG"=>   -1.5, "CT" => -1.5,  #revcomp
    "AC"=>   -1.6, "GT" => -1.6,  #revcomp
    "AA"=>   -0.2, "TT" => -0.2   #revcomp
)

DNA_DNA_hybrids = Dict(
    "AA"=>   -1.00, "TT" => -1.00, #revcomp
    "AT"=>   -0.88,
    "TA"=>   -0.58,
    "CA"=>   -1.45, "TG" => -1.45, #revcomp
    "GT"=>   -1.44, "AC" => -1.44, #revcomp
    "CT"=>   -1.28, "AG" => -1.28, #revcomp
    "GA"=>   -1.30, "TC" => -1.30, #revcomp
    "CG"=>   -2.17,
    "GC"=>   -2.24,
    "GG"=>   -1.42, "CC" => -1.42  #revcomp
)

function get_model_params()
    dir = @__DIR__
    home_dir = joinpath(split(dir, '/')[1:end-1]...)
    # Initialize model and matrices
    layer1 = CSV.read("/$home_dir/data/salis/free_energy_coeffs.txt", DataFrame, header=false)[!, "Column1"]

    inters = open("/$home_dir/data/salis/model_intercept.txt") do file
        read(file, String)
    end
    inters = parse(Float64, inters)
    return layer1, inters

end

all_kmers(k) = vcat(join.(collect(Iterators.product(fill([DNA_A, DNA_C, DNA_G, DNA_T], k)...)))...) |> sort

### This function calcs the hamming dist between 2 sequences of = len
function hamming(seq1, seq2)
    if length(seq1) != length(seq2)
        throw(ArgumentError("Sequences are not of equal length! Sequence 1: $(length(seq1)), Sequence 2: $(length(seq2))"))
    end
    distance = sum([x1 != x2 for (x1, x2) in zip(collect(seq1), collect(seq2))])
    return distance
end

function get_rsquared(TXpred, TXempirical)
    ols = lm(@formula(A  ~ B), DataFrame(A=TXpred, B=TXempirical))
    return r2(ols)
end

function AT_content(seq)
    ATcont = count(x -> x == DNA_A, collect(seq)) + count(x -> x == DNA_T, collect(seq))

    return ATcont / length(seq) * 100
end

# Functions to calculate DNA or RNA properties
function calc_stacking_and_twist_energy(seq)
    #find stacking free enegy of spacer
    spacer_sfe  = 0
    spacer_twist = 0
    for j in range(1, stop=length(seq), step=2)
        mer = seq[j:j+1]
        spacer_sfe += stacking_dict[string(mer)]
        spacer_twist += average_twist[string(mer)]
    end
    return spacer_sfe, spacer_twist
end

# local DNA element rigidity
function calc_rigidity(seq)
    rigidity = 0
    for m in range(1, stop=length(seq)-1, step=2)
        mer = seq[m:m+1]
        rigidity += persistence[string(mer)]
    end
    rigidity = rigidity / length(seq)
    return rigidity
end

function calc_groove_width(seq)
    groove_width = 0
    for l in range(1, stop=length(seq)-1, step=2)
        mer1 = seq[l:l+1]
        groove_width += groove_access[string(mer1)]
    end
    return groove_width
end

function calc_DNA_RNA_hybrid_energy(seq)
    dg_dna = 0
    dg_rna = 0
    for j in range(1, stop=length(seq[1:15])-1, step=2)
         # 15 is the length of long abortive transcripts
        mer = seq[j:j+1]
        dg_dna += DNA_DNA_hybrids[string(mer)]
        dg_rna += RNA_DNA_hybrids[string(mer)]
    end
    dg_hybrid = dg_dna - dg_rna
    return dg_dna, dg_rna, dg_hybrid
end

function get_matrices(coeffs)
    three_mer = LongDNA{4}.(all_kmers(3))
    two_mer = LongDNA{4}.(all_kmers(2))
    #Extract dG values from model coefficients
    dg10_0 = Dict(three_mer .=> coeffs[1:64])
    dg10_3 = Dict(three_mer .=> coeffs[65:128])
    dg35_0 = Dict(three_mer .=> coeffs[129:192])
    dg35_3 = Dict(three_mer .=> coeffs[193:256])
    dmers = Dict(three_mer .=> coeffs[257:256+64])
    x10mers = Dict(two_mer .=> coeffs[256+64+1:256+64+16])
    spacers = Dict(string.(collect(16:18)) .=> [256+64+16+1:256+64+16+3])

    return dg10_0, dg10_3, dg35_0, dg35_3, dmers, x10mers, spacers
end

function linear_free_energy_model(
    UP,
    h35,
    spacer,
    h10,
    disc,
    ITR,
    dg10_0,
    dg10_3,
    dg35_0,
    dg35_3,
    dmers,
    x10mers,
    coeffs,
    inters
    )

    prox_UP = UP[end-Int(floor(length(UP)/2))+1:end]
    dist_UP = UP[1:Int(floor(length(UP)/2))]

    # CATEGORICAL FEATURES
    ext10           = spacer[end-2:end-1] # TGN motif, contacts sigma
    hex10_0         = h10[1:3]
    hex10_3         = h10[4:end]
    hex35_0         = h35[1:3]
    hex35_3         = h35[4:end]
    disc_first_3mer = disc[1:3]
    spacer_length   = string(length(spacer))

    # NUMERICAL FEATURES
    dg_dna, dg_rna, dg_ITR     = calc_DNA_RNA_hybrid_energy(ITR) # calc R-loop strength
    rigidity                   = calc_rigidity(join([UP, h35, spacer[1:14]]))

    width_proxy_prox = calc_groove_width(prox_UP)
    width_proxy_dist = calc_groove_width(dist_UP)

    # NORMALIZE NUMERICAL FEATURES BY MAX IN TRAINING SET
    numericals         = [width_proxy_dist, width_proxy_prox, dg_ITR, rigidity]
    normalizing_values = [256.0, 255.0, 4.300000000000002, 25.780434782608694]
    numerical_coefs    = coeffs[end-3:end]
    normald            = numericals ./ normalizing_values
    dg_numerical       = normald .* numerical_coefs

    dg10      = dg10_0[hex10_0] + dg10_3[hex10_3]
    dg35      = dg35_0[hex35_0] + dg35_3[hex35_3]
    dg_disc   = dmers[disc_first_3mer]
    dg_ITR    = dg_numerical[end-1]
    dg_ext10  = x10mers[ext10]
    x = parse(Float64, spacer_length)
    dg_spacer = 0.1463*x^2 - 4.9113*x + 41.119

    dg_UP        = dg_numerical[1] + dg_numerical[2] + dg_numerical[end]
    dG_apparent  = (dg10 + dg35 + dg_disc + dg_ITR + dg_ext10 + dg_spacer + dg_UP + inters[1] - LOGK)/BETA
    dG_total     = dg10 + dg35 + dg_disc + dg_ITR + dg_ext10 + dg_spacer + dg_UP + inters[1]

    return dG_total, dG_apparent, dg10, dg35, dg_disc, dg_ITR, dg_ext10, dg_spacer, dg_UP
end

struct Promoter_Calculator
    model
    inters
    dg10_0
    dg10_3
    dg35_0
    dg35_3
    dmers
    x10mers
    spacers
    K
    BETA
    function Promoter_Calculator()
        layer1, inters = get_model_params()
        dg10_0, dg10_3, dg35_0, dg35_3, dmers, x10mers, spacers = get_matrices(layer1)
        K = 42.00000
        BETA = 1.636217004872062
        return new(layer1, inters, dg10_0, dg10_3, dg35_0, dg35_3, dmers, x10mers, spacers, K, BETA)
    end
end

function (p::Promoter_Calculator)(sequence::BioSequences.LongDNA{4}, TSS_range::Tuple{Int64, Int64}=(-1, -1))

    if TSS_range == (-1, -1)
        TSS_range = (1, length(sequence))
    end

    TSS_range_rev = (length(sequence) - TSS_range[2]+1, length(sequence) - TSS_range[1]+1)
    fwd_sequence = sequence
    rev_sequence = reverse_complement(sequence)
    Forward_Min_States, Forward_All_States = predict(p, fwd_sequence, TSS_range)
    Reverse_Min_States_Temp, Reverse_All_States_Temp = predict(p, rev_sequence, TSS_range_rev)

    Reverse_Min_States = Dict()
    Reverse_All_States = Dict()

    for TSS in Reverse_Min_States_Temp |> keys |> collect
        Reverse_Min_States[length(sequence) - TSS] = Reverse_Min_States_Temp[TSS]
        Reverse_All_States[length(sequence) - TSS] = Reverse_All_States_Temp[TSS]
    end
    Forward_Predictions_per_TSS = Forward_Min_States
    Reverse_Predictions_per_TSS = Reverse_Min_States
    output = Dict(
        "K" => p.K,
        "beta" => p.BETA,
        "sequence" => sequence,
        "TSS_range" => TSS_range,
        "Forward_Predictions_per_TSS" => Forward_Predictions_per_TSS,
        "Reverse_Predictions_per_TSS" => Reverse_Predictions_per_TSS
    )
    return output
end

function predict(p::Promoter_Calculator, sequence::BioSequences.LongDNA{4}, TSS_range::Tuple{Int64, Int64})
    UPS_length = 24
    HEX35_length = 6
    UPS_HEX35_SPACER = 1
    SPACER_length_range = [15, 20]
    HEX10_length = 6
    DISC_length_range = [6, 10]
    ITR_length = 20

    MinPromoterSize = UPS_length + UPS_HEX35_SPACER + HEX35_length + SPACER_length_range[1] + HEX10_length + DISC_length_range[1] + ITR_length
    MaxPromoterSize = UPS_length + UPS_HEX35_SPACER + HEX35_length + SPACER_length_range[2] + HEX10_length + DISC_length_range[2] + ITR_length
    MinimumTSS = UPS_length + UPS_HEX35_SPACER + HEX35_length + SPACER_length_range[1] + HEX10_length + DISC_length_range[1]

    All_States = Dict()
    Min_States = Dict()
    for TSS in TSS_range[1]:TSS_range[2]
        All_States[TSS] = Dict()
        for DISC_length in DISC_length_range[1]:DISC_length_range[2]
            if (TSS - DISC_length >= 0) && (TSS + ITR_length <= length(sequence))
                tempdisc = sequence[TSS - DISC_length + 1: TSS]
                tempITR  = sequence[TSS + 1 : TSS + ITR_length]

                for SPACER_length in SPACER_length_range[1]:SPACER_length_range[2]
                    if TSS - DISC_length - HEX10_length - SPACER_length - HEX35_length - UPS_length - UPS_HEX35_SPACER >= 0

                        temp10     = sequence[ TSS - DISC_length - HEX10_length + 1 : TSS - DISC_length ]
                        tempspacer = sequence[ TSS - DISC_length - HEX10_length - SPACER_length + 1 : TSS - DISC_length - HEX10_length]
                        temp35     = sequence[( TSS - DISC_length - HEX10_length - SPACER_length - HEX35_length +1): TSS - DISC_length - HEX10_length - SPACER_length]

                        tempUP     = sequence[ TSS - DISC_length - HEX10_length - SPACER_length - HEX35_length - UPS_length - UPS_HEX35_SPACER + 1:  TSS - DISC_length - HEX10_length - SPACER_length - HEX35_length - UPS_HEX35_SPACER]

                        dG_total, dG_apparent, dG_10, dG_35, dG_disc, dG_ITR, dG_ext10, dG_spacer, dG_UP = linear_free_energy_model(tempUP, temp35, tempspacer, temp10, tempdisc, tempITR, p.dg10_0, p.dg10_3, p.dg35_0, p.dg35_3, p.dmers, p.x10mers, p.model, p.inters)

                        dG_bind  = dG_10 + dG_35 + dG_spacer + dG_ext10 + dG_UP

                        Tx_rate = p.K * exp(- p.BETA * dG_total)

                        result = Dict(
                            "promoter_sequence" => sequence[TSS - DISC_length - HEX10_length - SPACER_length - HEX35_length - UPS_length - UPS_HEX35_SPACER + 1: TSS + ITR_length],
                            "TSS" => TSS,
                            "UP" => tempUP,
                            "hex35" => temp35,
                            "spacer" => tempspacer,
                            "hex10" => temp10,
                            "disc" => tempdisc,
                            "ITR" => tempITR,
                            "dG_total" => dG_total,
                            "dG_10" => dG_10,
                            "dG_35" => dG_35,
                            "dG_disc" => dG_disc,
                            "dG_ITR" => dG_ITR,
                            "dG_ext10"=> dG_ext10,
                            "dG_spacer" => dG_spacer,
                            "dG_UP" => dG_UP,
                            "dG_bind" => dG_bind,
                            "Tx_rate" => Tx_rate,
                            "UP_position" => [TSS - DISC_length - HEX10_length - SPACER_length - HEX35_length - UPS_HEX35_SPACER - UPS_length + 1,
                                                       TSS - DISC_length - HEX10_length - SPACER_length - HEX35_length - UPS_HEX35_SPACER],
                            "hex35_position" => [TSS - DISC_length - HEX10_length - SPACER_length - HEX35_length + 1, TSS - DISC_length - HEX10_length - SPACER_length],
                            "spacer_position" => [TSS - DISC_length - HEX10_length - SPACER_length + 1, TSS - DISC_length - HEX10_length],
                            "hex10_position" => [TSS - DISC_length - HEX10_length + 1, TSS - DISC_length],
                            "disc_position" => [TSS - DISC_length + 1, TSS]
                        )

                        All_States[TSS][(DISC_length, SPACER_length)] = result
                        if TSS in Min_States |> keys |> collect
                            if result["dG_total"] < Min_States[TSS]["dG_total"]
                                Min_States[TSS] = result
                            end
                        else
                            Min_States[TSS] = result
                        end
                    end
                end
            end
        end
    end
    return (Min_States, All_States)
end


text:
"""
    triangle_range(a,b)

Find all (even) `k` such that `|a-b| ≤ k ≤ a + b`. This is useful when
expanding matrix elements of tensors between angular momenta `a` and
`b` in multipoles `k`; `triangle_range` can then be used to decided
which multipole terms are required.
"""
function triangle_range(a,b)
    kmin = abs(a-b)
    if !iseven(kmin + a + b)
        kmin += 1
    end
    kmin:2:(a+b)
end

"""
    ∏(ℓs...)

Calculates `√((2ℓ₁+1)(2ℓ₂+1)...(2ℓₙ+1))`, which is a common factor in
angular momentum algebra.
"""
function ∏(ℓs...)
    v = 1
    for ℓ in ℓs
        v *= 2ℓ + 1
    end
    √v
end

"""
    powneg1(k) = (-)ᵏ

Calculates powers of negative unity for integer `k`.
"""
powneg1(k::Integer) = isodd(k) ? -1 : 1

"""
    jmⱼ(o::SpinOrbital)

Return the angular momentum and its projection on the z axis of the
spin-orbital `o`.
"""
jmⱼ(o::SpinOrbital{<:Orbital}) = o.orb.ℓ, o.m[1]
jmⱼ(o::SpinOrbital{<:RelativisticOrbital}) = o.orb.j, o.m[1]

"""
    spin(o::SpinOrbital)

Return the spin of the spin-orbital `o`.
"""
spin(o::SpinOrbital{<:Orbital}) = o.m[2]

"""
    @δ (a,b)[, (c,d) ...]

Kronecker ``\\delta_{ab}\\delta_{cd}...`` that tests each pair of
values for equality and quick-returns `0` at the first inequality.
"""
macro δ(vars...)
    code = map(vars) do v
        :(isequal($(esc(v))...) || return 0)
    end
    Expr(:block, code..., :(1))
end

"""
    couplings(j₁, m₁, j₂, m₂)

Return all `j`, `m` for which the Clebsch–Gordan coefficient
`⟨j₁m₁ j₂m₂|jm⟩` is non-zero.
"""
function couplings(j₁, m₁, j₂, m₂)
    for (jᵢ,mᵢ) in ((j₁, m₁), (j₂, m₂))
        WignerSymbols.ϵ(jᵢ, mᵢ) || throw(DomainError((jᵢ, mᵢ), "invalid combination (jᵢ, mᵢ)"))
    end
    m = m₁ + m₂
    a = max(abs(j₁-j₂),abs(m))
    b = j₁ + j₂
    js = if iszero(m₁) && iszero(m₂) && iszero(m)
        (iseven(a) ? a : a + 1):2:b
    else
        a:b
    end
    js, m
end


text:
"""
    Ensemble

A supertype for all kinds of ensembles.
"""
abstract type Ensemble end

"""
    EnsembleOptions{T,E}

A type to hold necessary information for each kind of `Ensemble`.

# Fields
- `ensemble::E`: the ensemble that will specify this type; normally a sub-type of
    `Ensemble`.
- `nattempt::T`: the number of attempts that each ensemble has been called.
- `naccept::T`: the number of accepted attempts for each ensemble called.
"""
mutable struct EnsembleOptions{T,E}
    ensemble::E
    nattempt::T
    naccept::T
end

EnsembleOptions(e::Ensemble) = EnsembleOptions(e, 0, 1)


text:
using Test
using jInvVis
using jInv.Mesh

# tests for regular mesh
domain = [0 1.1 0 1.0 0  1.1]
n      = [8 5 3]
Mr    = getRegularMesh(domain,n)
Mt    = getTensorMesh3D(rand(n[1]),rand(n[2]),rand(n[3]))
I     = randn(tuple(n...))
println("=== test viewOrthoSlices  ===")
figure(2)
viewOrthoSlices2D(I,Mr)
figure(1)
viewOrthoSlices2D(I,Mt)
println("==== all passed === ")


text:
using ReactionNetworkImporters
using DiffEqBase, DiffEqBiological, DiffEqJump, OrdinaryDiffEq, Plots

using TimerOutputs

plotlyjs()

# Create a TimerOutput, this is the main type that keeps track of everything.
const to = TimerOutput()

# parameters
method = RSSA()
nsteps = 1e5                  # time to simulate to
networkname = "tester"
tf = 10.
#speciesf = path to species initial population file
#rxsf = path to reaction network file

# get the RSSA reaction network
reset_timer!(to)
@timeit to "netgen" prn = loadrxnetwork(RSSANetwork(), networkname, speciesf, rxsf; printrxs = false)
rn = prn.rn
initialpop = prn.u₀
println("network parsed")

# get the BioNetGen reaction network
@timeit to "bionetgen" prnbng = loadrxnetwork(BNGNetwork(),string(networkname,"bng"), bngfname)
rnbng = prnbng.rn; u₀ = prnbng.u₀; p = prnbng.p; shortsymstosyms = prnbng.symstonames;

u0 = round.(Int, u₀)
println(typeof(u0))

# one simulation
@timeit to "addjumps" addjumps!(rn,build_regular_jumps=false, minimal_jumps=true)
println("added jumps")
@timeit to "dprob" rprob = DiscreteProblem(rn, u0, (0.,tf))
@timeit to "jprob" rjprob = JumpProblem(rprob, method, rn; save_positions=(false,false))
@timeit to "solve" rsol = solve(rjprob, SSAStepper(), saveat=tf/1000)
show(to)

# bng file
reset_timer!(to)
@timeit to "addjumps" addjumps!(rnbng,build_regular_jumps=false, minimal_jumps=true)
@timeit to "dprob" bprob = DiscreteProblem(rnbng, u0, (0.,tf), p)
@timeit to "jprob" bjprob = JumpProblem(bprob, method, rnbng; save_positions=(false,false))
@timeit to "solve" bsol = solve(bjprob, SSAStepper(), saveat=tf/1000)
show(to)

# plot
#labs = [String(sym) for i=1:1, sym in rn.syms];
#plot(sol.t, sol[:,:]', linetype=:steppost, labels=labs)
#plot(sol.t, sol[:,:]', labels=labs)

rasyk = sum(rsol[prnbng.groupstoids[:Activated_Syk],:], dims=1)
basyk = sum(bsol[prnbng.groupstoids[:Activated_Syk],:], dims=1)
plot(rsol.t, rasyk', label=:rsol)
plot!(bsol.t, basyk', label=:bsol)


text:
"""
    update_best_bound!(com::CS.CoM)

Iterate over all backtrack objects to set the new best bound for the whole search tree
"""
function update_best_bound!(com::CS.CoM)
    if !isempty(com.backtrack_pq)
        if com.sense == MOI.MIN_SENSE
            max_val = typemax(com.best_bound)
            com.best_bound = minimum(
                bo.status == :Open ? bo.best_bound : max_val for bo in com.backtrack_vec
            )
        elseif com.sense == MOI.MAX_SENSE
            min_val = typemin(com.best_bound)
            com.best_bound = maximum(
                bo.status == :Open ? bo.best_bound : min_val for bo in com.backtrack_vec
            )
        end # otherwise no update is needed
    end
end

"""
    get_best_bound(com::CS.CoM, backtrack_obj::CS.BacktrackObj, obj_fct::SingleVariableObjective, vidx::Int, lb::Int, ub::Int)

Compute the best bound if we have a `SingleVariableObjective` and limit `vidx` using
    `lb <= var[vidx] <= ub` if `vidx != 0`.
Return a best bound given the constraints on `vidx`
"""
function get_best_bound(
    com::CS.CoM,
    backtrack_obj::CS.BacktrackObj,
    obj_fct::SingleVariableObjective,
    vidx::Int,
    lb::Int,
    ub::Int,
)
    if obj_fct.vidx != vidx
        if com.sense == MOI.MIN_SENSE
            return com.search_space[obj_fct.vidx].min
        else # MAX
            return com.search_space[obj_fct.vidx].max
        end
    else
        if com.sense == MOI.MIN_SENSE
            return lb
        else
            return ub
        end
    end
end

"""
    get_best_bound(com::CS.CoM, backtrack_obj::CS.BacktrackObj, obj_fct::LinearCombinationObjective, vidx::Int, lb::Int, ub::Int)

Compute the best bound if we have a `LinearCombinationObjective` and limit `vidx` using
    `lb <= var[vidx] <= ub` if `vidx != 0`.
Return a best bound given the constraints on `vidx`
"""
function get_best_bound(
    com::CS.CoM,
    backtrack_obj::CS.BacktrackObj,
    obj_fct::LinearCombinationObjective,
    vidx::Int,
    lb::Int,
    ub::Int,
)
    indices = obj_fct.lc.indices
    coeffs = obj_fct.lc.coeffs
    objval = obj_fct.constant
    if com.sense == MOI.MIN_SENSE
        for i in 1:length(indices)
            if indices[i] == vidx
                objval += min(coeffs[i] * lb, coeffs[i] * ub)
                continue
            end
            objval += min(
                coeffs[i] * com.search_space[indices[i]].min,
                coeffs[i] * com.search_space[indices[i]].max,
            )
        end
    else # MAX Sense
        for i in 1:length(indices)
            if indices[i] == vidx
                objval += max(coeffs[i] * lb, coeffs[i] * ub)
                continue
            end
            objval += max(
                coeffs[i] * com.search_space[indices[i]].min,
                coeffs[i] * com.search_space[indices[i]].max,
            )
        end
    end

    return objval
end

function get_best_bound_lp(com, backtrack_obj, vidx, lb, ub)
    # check if last best_bound is affected
    # check that we have a parent node to maybe use the bound of the parent
    if backtrack_obj.parent_idx != 0 && vidx == 0
        use_last = true
        for variable in com.search_space
            lb = com.search_space[variable.idx].min
            ub = com.search_space[variable.idx].max
            if lb > backtrack_obj.primal_start[variable.idx] ||
               ub < backtrack_obj.primal_start[variable.idx]
                use_last = false
            end
        end
        # best_bound is the best_bound from parent
        if use_last
            backtrack_obj.solution = backtrack_obj.primal_start
            return backtrack_obj.best_bound
        end
    end

    # compute bound using the lp optimizer
    # setting all bounds
    for variable in com.search_space
        if variable.idx == vidx
            set_lower_bound(com.lp_x[vidx], lb)
            set_upper_bound(com.lp_x[vidx], ub)
        else
            set_lower_bound(com.lp_x[variable.idx], com.search_space[variable.idx].min)
            set_upper_bound(com.lp_x[variable.idx], com.search_space[variable.idx].max)
        end
    end
    lp_backend = backend(com.lp_model)
    if MOI.supports(lp_backend, MOI.VariablePrimalStart(), MOI.VariableIndex)
        for variable in com.search_space
            v_idx = variable.idx
            MOI.set(
                lp_backend,
                MOI.VariablePrimalStart(),
                MOI.VariableIndex(v_idx),
                backtrack_obj.primal_start[v_idx],
            )
        end
    end

    # update bounds by constraints
    # check each constraint which has `update_best_bound = true` for a better bound
    for constraint in com.constraints
        update_best_bound_constraint!(
            com,
            constraint,
            constraint.fct,
            constraint.set,
            vidx,
            lb,
            ub,
        )
        for bound in constraint.bound_rhs
            set_lower_bound(com.lp_x[bound.idx], bound.lb)
            set_upper_bound(com.lp_x[bound.idx], bound.ub)
        end
    end

    optimize!(com.lp_model)
    if termination_status(com.lp_model) == MOI.OPTIMAL
        backtrack_obj.solution = JuMP.value.(com.lp_x)
        return objective_value(com.lp_model)
    end
    if com.sense == MOI.MIN_SENSE
        return typemax(com.best_bound)
    end
    return typemin(com.best_bound)
end


text:
@testset "default solvers" begin
    prob = MVP.UnconstrainedProblems.examples["Powell"]
    f = objective(prob)
    g! = gradient(prob)
    h! = hessian(prob)
    @test summary(optimize(f, prob.initial_x)) == summary(NelderMead())
    @test summary(optimize(f, g!, prob.initial_x)) == summary(LBFGS())
    @test summary(optimize(f, g!, h!, prob.initial_x)) == summary(Newton())
end


text:
module ParrotHelperFlagBerryMoon

using ..Ahorn, Maple

@mapdef Entity "ParrotHelper/FlagBerryMoon" FlagBerryMoon(x::Integer, y::Integer, flag::String="", set::Bool=true, nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const placements = Ahorn.PlacementDict(
    "FlagBerry (Moon, ParrotHelper)" => Ahorn.EntityPlacement(
        FlagBerryMoon
    )
)

Ahorn.nodeLimits(entity::FlagBerryMoon) = 0, -1

sprite = "collectables/strawberry/normal00"

function Ahorn.selection(entity::FlagBerryMoon)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ())
    winged = get(entity.data, "winged", false)
    if winged
        sprite = "collectables/strawberry/wings01"
    else
        sprite = "collectables/moonBerry/normal00"
    end
    hasPips = length(nodes) > 0
    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]

    if hasPips
        for node in nodes
            nx, ny = node

            push!(res, Ahorn.getSpriteRectangle("collectables/strawberry/seed00", nx, ny))
        end
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::FlagBerryMoon)
    x, y = Ahorn.position(entity)

    for node in get(entity.data, "nodes", ())
        nx, ny = node

        Ahorn.drawLines(ctx, Tuple{Number, Number}[(x, y), (nx, ny)], Ahorn.colors.selection_selected_fc)
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FlagBerryMoon, room::Maple.Room)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ())
    winged = get(entity.data, "winged", false)
    hasPips = length(nodes) > 0

    if winged
        sprite = "collectables/strawberry/wings01"
    else
        sprite = "collectables/moonBerry/normal00"
    end

    for node in nodes
        nx, ny = node

        Ahorn.drawSprite(ctx, "collectables/strawberry/seed00", nx, ny)
    end

    Ahorn.drawSprite(ctx, sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagBerryMoon, room::Maple.Room)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())
    winged = get(entity.data, "winged", false)
    hasPips = length(nodes) > 0

    if winged
        sprite = "collectables/strawberry/wings01"
    else
        sprite = "collectables/moonBerry/normal00"
    end
    if hasPips
        for node in nodes
            nx, ny = node
            Ahorn.drawSprite(ctx, "collectables/strawberry/seed00", nx, ny)
        end
    end

    Ahorn.drawSprite(ctx, sprite, x, y)
end

end

max_stars_repo_name: Parrotapocalypse/ParrotHelper

max_issues_repo_name: Parrotapocalypse/ParrotHelper

max_forks_repo_name: Parrotapocalypse/ParrotHelper


text:
#=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=#
# Problem set 1 solutions
# Written by Tyler Ransom
# Commented by Giuseppe Grasso
#=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=#

using DataFrames
using FreqTables
using CSV
using JLD
using Random
using LinearAlgebra
using Statistics
using Distributions
cd("/Users/peppegrass/Documents/GitHub/fall-2020/ProblemSets/PS1-julia-intro/") # GG: sets working directory
pwd() ## GG: prints working directory
readdir() # GG: equivalent to -ls- to see elements of working directory
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1
# Initializing variables and practice with basic matrix operations
#:::::::::::::::::::::::::::::::::::::::::::::::::::

# a. Set the seed
Random.seed!(1234)
#=
GG: It seems to me that Random.seed!() only applies to the next command.
Hence the use of function q1() so that it applies to everything.
But I'm not sure about it.
=#

#=
GG: The content of the function cannot be called from the .jl
You either copy and paste it in the REPL or call the functions (as done in the end of this script)
=#
function q1()
    ## a. i.-iv. create matrices

    # A 10×7- random numbers distributed U [−5,10]
    A = 15*rand(10,7).-5 # GG: same result as using Distributions A=rand(Uniform(-5,10),10,7). Not sure why you have 15 instead of 10 here, but it works.

    # B 10×7- random numbers distributed N (−2,15) [st dev is 15]
    B = 15*randn(10,7).-2 # GG: same result as using Distributions B=rand(Normal(-2,15),10,7)

    # C 5×7- the ﬁrst 5 rows and ﬁrst 5 columns of A and the last two columns and ﬁrst 5 rows of B
    C = cat(A[1:5,1:5],B[1:5,end-1:end]; dims=2) # GG: notice end arguments

    # D 10×7- where D i,j= A i,jif A i,j≤ 0, or 0 otherwise
    D = A.*(A.<=0) # GG: D contains only negative elements of A; rest filled with 0s
    # GG: very elegant subsetting; had no idea how to do it

    ## b. Use a built-in Julia function to list the number of elements of A
    println(length(A))

    ## c. Use a series of built-in Julia functions to list the number of unique elements of D
    println(length(unique(A)))

    ## d. Using the reshape() function, create a new matrix called E which is the ‘vec’ operator applied to B.
    ## Can you ﬁnd an easier way to accomplish this?

    # one way:
    E = reshape(B,length(B)) # GG: piles columns of B one above the other
    # easy way:
    E = B[:] ## GG: simple as that

    ## e. Create a new array called F which is 3-dimensional and contains
    ## A in the ﬁrst column of the third dimension and
    ## B in the second column of the third dimension
    F = cat(A,B; dims=3) # GG: making a cuboid (3-dimensional tensor) out of matrices A and B; see https://improbable-research.github.io/keanu/docs/tensors/

    ## f. Use the permutedims() function to twist F so that it is now F 2×10×7 instead of F 10×7×2 . Save this new matrix as F.
    F = permutedims(F, [3 1 2]) # GG: now F has 7 third dimensions, each containing a 2x10 matrix

    ## g. Create a matrix G which is equal to B ⊗ C (the Kronecker product of B and C). What happens when you try C ⊗ F?
    G = kron(B,C) # GG: A (5x10=50)x(7x7=49) matrix; see https://en.wikipedia.org/wiki/Kronecker_product
    #kron(C,F) # it does work; not sure whether it is that tensor of different dimension cannot be kronecher-multiplied (don't think so) or that this Julia command isn't suited for such operation

    ## h. Save the matrices A, B, C, D, E, F and G as a .jld ﬁle named matrixpractice.
    save("matrixpractice.jld","A",A,"B",B,"C",C,"D",D,"E",E,"F",F,"G",G)
        #=
        GG: "matrixpractice.jld" can be read back with
        d = load("matrixpractice.jld")
        where d is a dictionary
        Other examples of .jld loading are provided in Exercise 4
        =#

    ## i. Save only the matrices A, B, C, and D as a .jld ﬁle called firstmatrix.
    save("firstmatrix.jld","A",A,"B",B,"C",C,"D",D)

    ## j. Export C as a .csv ﬁle called Cmatrix. You will ﬁrst need to transform C into a DataFrame.
    CSV.write("Cmatrix.csv",DataFrame(C))

    ## k. Export D as a tab-delimited .dat ﬁle called Dmatrix. You will ﬁrst need to transform D into a DataFrame.
    CSV.write("Dmatrix.dat",DataFrame(D); delim="\t")
    return A,B,C,D
end
#=
l. Wrap a function deﬁnition around all of the code for question 1. Call the function q1().
The function should have 0 inputs and should output the arrays A, B, C and D.
At the very bottom of your script you should add the code A,B,C,D = q1().
=#

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2
# Practice with loops and comprehensions
#:::::::::::::::::::::::::::::::::::::::::::::::::::

function q2(A,B,C)
    ## a. Write a loop or use a comprehension that computes the element-by-element product of A and B. Name the new matrix AB.
    ## Create a matrix called AB2 that accomplishes this task without a loop or comprehension.

    #=
    GG: Using a COMPREHENSION to compute the element-by-element product of A and B; which has same dimension 10x7
    What is a comprehension? It allows the user in 1 line of code to create an object that could be a complex formula.
    e.g. computing a present value by PV = sum([β^t*Y[t] for t=1:T])
    or a shift-share instrument for what matters
    =#
    AB = [A[i,j]*B[i,j] for i=1:size(A,1),j=1:size(A,2)] # GG: could have been size(B,1) and size(B,2); what matter are the 1 and 2 to indicate dimensions

    AB2 = A.*B ## GG: doing it with the ELEMENT-BY-ELEMENT multiplication oparator .*

    ## b. Write a loop that creates a column vector called Cprime which contains only the elements of C that are between -5 and 5 (inclusive).
    ## Create a vector called Cprime2 which does this calculation without a loop.
    Cprime = [] # GG: creating empty array
    for j=1:size(C,2)
        for i=1:size(C,1)
            if C[i,j]>=-5 && C[i,j]<=5 ## GG: declaring the IF condition here; note usage of && coz boolean here
                push!(Cprime,C[i,j]) ## GG: sounds like some sort of "copy"; "pushes" only element satisfying IF clause
            end
        end
    end
    Cprime2 = C[(C.>=-5) .& (C.<=5)]

    ## c. Using loops or comprehensions, create a 3-dimensional array called X that is of dimension N × K × T where
    ## N = 15,169, K = 6, and T = 5. For all t, the columns of X should be (in order):
    ## - an intercept (i.e. vector of ones)
    ## - a dummy variable that is 1 with probability .75∗(6−t)/5
    ## - a continuous variable distributed normal with mean 15+t−1 and standard deviation 5(t−1)
        ## GG: apparently (mean).+(sd).*randn(N,1) generates a normal rv with mean and sd
    ## - a continuous variable distributed normal with mean π(6−t)/3 and standard deviation 1/e
    ## - a discrete variable distributed “discrete normal” with mean 12 and standard deviation 2.19. (A discrete normal random variable is properly called a binomial random variable. The distribution described above can be implemented by choosing binomial parameters n and p where n = 20 and p = 0.6. Use the following code (after loading Julia’s Distributions package) to generate this vector of X: rand(Binomial(20,0.6),N), where N is the length of the vector
    ## - a discrete variable distributed binomial with n = 20 and p = 0.5
    N = 15169
    K = 6
    T = 5
    #=
    GG: cat is for CONCATENATION
    WOW this is basically generating a dataset in one line :O
    The dataset has N=15169 observations, K=6 variables, and T=5 time periods
    =#
    X = cat([cat([ones(N,1) rand(N,1).<=(0.75*(6-t)/5) (15+t-1).+(5*(t-1)).*randn(N,1) (π*(6-t)/3).+(1/exp(1)).*randn(N,1) rand(Binomial(20,0.6),N) rand(Binomial(20,0.5),N)];dims=3) for t=1:T]...;dims=3) # discrete_normal binomial

    ## d. Use comprehensions to create a matrix β which is K × T and whose elements evolve across time in the following fashion:
    ## - 1,1.25,1.5,...
    ## - ln(t)
    ## - −√t
    ## - exp{t}-exp{t+1}
    ## - t
    ## - t/3
    β = vcat([cat([1+0.25*(t-1) log(t) -sqrt(t) exp(t)-exp(t+1) t t/3];dims=1) for t=1:T]...)' ## GG: using vertical concatenation here; mind the transpose

    ## e. Use comprehensions to create a matrix Y which is N × T deﬁned by Y_t= X_t β_t + ε_t , where ε_t ∼ iid N (0,σ=.36)
    Y = hcat([cat(X[:,:,t]*β[:,t] + .36*randn(N,1);dims=2) for t=1:T]...) # GG: 3-dots at the end have to do with the layout given by hcat
    return nothing
end

#=
f. Wrap a function deﬁnition around all of the code for question 2. Call the function q2().
The function should have take as inputs the arrays A, B and C. It should return nothing.
At the very bottom of your script you should add the code q2(A,B,C). Make sure q2() gets called after q1()!
=#

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3
# Reading in Data and calculating summary statistics
#:::::::::::::::::::::::::::::::::::::::::::::::::::

function q3()
    ## a. Clear the workspace and import the ﬁle nlsw88.csv into Julia as a DataFrame.
    ## Make sure you appropriately convert missing values and variable names. Save the result as nlsw88.jld.
    nlsw88 = CSV.read("nlsw88.csv") # GG: kinda -import csv-
    save("nlsw88.jld","nlsw88",nlsw88) # GG: kinda -save-

    ## b. What percentage of the sample has never been married? What percentage are college graduates?
    mean(nlsw88.never_married) # GG: kinda -summarize- but just returning the mean
    mean(nlsw88.collgrad)

    ## c. Use the tabulate command to report what percentage of the sample is in each race category
    freqtable(nlsw88, :race) # GG: oneway -tabulate race-; gotta have that colon to work though

    ## d. Use the describe() function to create a matrix called summarystats which lists the mean, median,
    ## standard deviation, min, max, number of unique elements, and interquartile range (75th percentile minus
    ## 25th percentile) of the data frame. How many grade observations are missing?
    summarystats = describe(nlsw88)
    #=
    GG: to see how many grade observation are missing you gotta visualize all rows using
    show(summarystats, allrows=true)
    see https://juliadata.github.io/DataFrames.jl/stable/man/getting_started/
    =#

    ## e. Show the joint distribution of industry and occupation using a cross-tabulation.
    freqtable(nlsw88, :industry, :occupation) # GG: twoway tabulate

    ## f. Tabulate the mean wage over industry and occupation categories.
    ## Hint: you should ﬁrst subset the data frame to only include the columns industry, occupation and wage.
    ## You should then follow the “split-apply-combine” directions here.
    wageonly = nlsw88[:,[:industry,:occupation,:wage]] # GG: [SPLIT] creating a sub-dataframe from nlsw88
    grouper = groupby(wageonly, [:industry,:occupation]) # GG: [APPLY] created a GroupedDataFrame whereby observations are grouped, e.g. 12 observations/rows with industry=5 and occupation=6
    combine(grouper, valuecols(grouper) .=> mean) # GG: [COMBINE] essentially "collapsing" to means by groups defined in previous step
    return nothing
end

#=
g. Wrap a function deﬁnition around all of the code for question 3. Call the function q3().
The function should have no inputs and no outputs. At the very bottom of your script you should add the code q3().
=#

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 4
# Practice with functions
#:::::::::::::::::::::::::::::::::::::::::::::::::::

function q4()
    ## a. Load firstmatrix.jld.
    mats = load("firstmatrix.jld") # GG: loads .jld file in dictionary "mats"
    A = mats["A"] # GG: creating/extracting array from "mats"
    B = mats["B"]
    C = mats["C"]
    D = mats["D"]

    ## b. Write a function called matrixops that takes as inputs the matrices A and B from question (a) of problem 1 and has
    ## three outputs: (i) the element-by-element product of the inputs, (ii) the product A ′ B, and (iii) the sum of all the
    ## elements of A+B.
    function matrixops(m1,m2) # GG: defining function matrixops to have two argumets/inputs
        ## c. Starting on line 2 of the function, write a comment that explains what matrixops does.
        ## this function computes the following matrix formulas: i) element-wise multiplication; ii) transpose multiplication; iii) sum of the elementwise product

        ## e. Just before the ﬁrst executable line of matrixops.m (i.e. right after the ﬁrst-line comments),
        ## write an if statement which gives an error if the two inputs are not the same size. Have the error say “inputs must have the same size.”
        if size(m1)!=size(m2)
            error("inputs must have the same size.") # GG: creating an error which stops function if matrices aren't the same size (needed for element-wise multiplication)
        end
        ret1 = m1.*m2
        ret2 = m1'*m2 # GG: '* is operator for transpose multiplication
        ret3 = sum(m1+m2)
        return ret1,ret2,ret3
    end

    ## d. Evaluate [GG: the returns one-ny-one in the REPL and the function itself] matrixops() using A and B from question (a) of problem 1
    matrixops(A,B) # GG: everything goes fine

    ## f. Evaluate matrixops.m using C and D from question (a) of problem 1. What happens?
    # matrixops(C,D) # GG: since C and D aren't the same size the programmed error is raised

    ## g. Now evaluate matrixops.m using ttl_exp and wage from nlsw88.jld.
    ## Hint: before doing this, you will need to convert the data frame columns to Arrays.
    ## e.g. convert(Array,nlsw88.ttl_exp), depending on what you called the data frame object [I called it nlsw88].
    #= GG: Two lines below actually not working. Apparently, it has to do with the conversion that are purported to do next...
    mat1 = load("nlsw88.jld")
    nlsw88 = mat1["nlsw88"]
    ... Instead a re-import from the CSV and create the DataFrame from scratch
    =#
    nlsw88 = CSV.read("nlsw88.csv",DataFrame;) # GG: nlsw88 = CSV.read("nlsw88.csv") would work to but appears to be deprecated in favour of the one used
    matrixops(convert(Array,nlsw88.ttl_exp),convert(Array,nlsw88.wage)) # GG: evaluating matrixops() on the converted versions of nlsw88.ttl_exp and nlsw88.wage
    return nothing
end

#=
h. Wrap a function deﬁnition around all of the code for question 4. Call the function q4().
The function should have no inputs or outputs. At the very bottom of your script you should add the code q4().
=#

# Call the functions defined above
A,B,C,D = q1()
q2(A,B,C)
q3()
q4()


text:
using Iris
using Test

@testset "Lattices" begin
    @testset "1D" begin
        @testset "Cartesian" begin
            lat = @test_nowarn Lattices.Cartesian(.313)
            @test lat[latticeindex(lat,π)] ≈ Point(π)
            @test all(lat[0:2] .≈ [Point(0),Point(.313),Point(.626)])
            lat = @test_nowarn Lattices.Cartesian(.1; origin = Point(2))
            @test lat[latticeindex(lat,π)] ≈ Point(π)
            @test_throws ArgumentError Lattices.Cartesian(.1; origin = Point(2), angles = .1)
        end
        @testset "Polar" begin
            lat = @test_nowarn Lattices.Polar(.313)
            @test lat[latticeindex(lat,π)] ≈ Point(π)
            @test all(lat[0:2] .≈ [Point(lat.dr/2),Point(3lat.dr/2),Point(5lat.dr/2)])
            lat = @test_nowarn Lattices.Polar(.1; origin = Point(2))
            @test lat[latticeindex(lat,π)] ≈ Point(π)
            @test_throws ArgumentError Lattices.Polar(.1; origin = Point(2), angles = .1)
        end
        @testset "Spherical" begin
            lat = @test_nowarn Lattices.Spherical(.313)
            @test lat[latticeindex(lat,π)] ≈ Point(π)
            @test all(lat[0:2] .≈ [Point(lat.dr/2),Point(3lat.dr/2),Point(5lat.dr/2)])
            lat = @test_nowarn Lattices.Spherical(.1; origin = Point(2))
            @test lat[latticeindex(lat,π)] ≈ Point(π)
            @test_throws ArgumentError Lattices.Spherical(.1; origin = Point(2), angles = .1)
        end
    end
    @testset "2D" begin
        @testset "Cartesian" begin
            lat = @test_nowarn Lattices.Cartesian(.313,.213)
            @test lat[latticeindex(lat,Point(π,2))] ≈ Point(π,2)
            @test all(lat[0:2,0:2] .≈ [Point(0,0),Point(.313,.213),Point(.626,.426)])
            @test all([lat[i,j] for i ∈ 0:1, j ∈ 0:1] .≈ [Point(0,0) Point(0,.213);Point(.313,0) Point(.313,.213)])
            lat = @test_nowarn Lattices.Cartesian(.1,.2; origin = Point(2,3))
            @test lat[latticeindex(lat,Point(π,sqrt(5)))] ≈ Point(π,sqrt(5))
            @test_nowarn Lattices.Cartesian(.1,.2; origin = Point(2,3), angles = .1)
        end
        @testset "Polar" begin
            lat = @test_nowarn Lattices.Polar(.1,.1)
            @test lat[latticeindex(lat,Point{Polar}(1,π))] ≈ Point(-1,0)
            @test all(lat[0:2,0:2] .≈ [Point{Polar}(.05,0),Point{Polar}(0.15,.1),Point{Polar}(.25,.2)])
            lat = @test_nowarn Lattices.Polar(.1; origin = Point(2))
            @test lat[latticeindex(lat,π)] ≈ Point(π)
            @test_throws ArgumentError Lattices.Polar(.1; origin = Point(2), angles = .1)
        end
        @testset "Spherical" begin
            lat = @test_nowarn Lattices.Spherical(1,π/4)
            @test lat[latticeindex(lat,1,π)] ≈ Point(0,-1)
            @test all(lat[0:2,0:2] .≈ [Point(0,1/2),Point(1.5/sqrt(2),1.5/sqrt(2)),Point(2.5,0)])
            lat = @test_nowarn Lattices.Spherical(.1,.2; origin = Point(2,1))
            @test lat[latticeindex(lat,π,π)] ≈ Point(0,-π)
            @test_throws ArgumentError Lattices.Spherical(.1; origin = Point(2), angles = .1)
        end
    end
end


text:
using Test
using Hello

@testset "hello" begin
    @test Hello.command_main(["-h"]) == 0
    @test Hello.command_main(["--help"]) == 0
    @test Hello.command_main(["2", "--opt1=3", "--opt2", "5", "-f"]) == 0
    @test Hello.command_main(["2", "--opt1=3", "--opt2", "5", "--flag"]) == 0
end


text:
"""Optimize a quantum control problem.

```julia
opt_result = optimize(problem; method=<method>, kwargs...)
```

optimizes towards a solution of given [`problem`](@ref ControlProblem) with the
given optimization `method`. Any keyword argument temporarily overrides the
corresponding keyword argument in `problem`.
"""
function optimize(problem::ControlProblem; method, kwargs...)
    if length(kwargs) > 0
        tempproblem = copy(problem)
        merge!(tempproblem.kwargs, kwargs)
        return optimize(tempproblem, method)
    else
        return optimize(problem, method)
    end
end

optimize(problem::ControlProblem, method::Symbol) = optimize(problem, Val(method))
#
# Note: Methods *must* be defined in the various optimization packages as e.g.
#
#   optimize(problem, method::Val{:krotov}) = optimize_krotov(problem)
#

import DrWatson

const _SAVENAME_AVAILABLE_KEYS = Set([
    :accesses,
    :allowedtypes,
    :connector,
    :digits,
    :equals,
    :expand,
    :ignores,
    :sigdigits,
    :sort,
    :val_to_string,
])
const DEFAULT_OPTIMIZATION_SAVENAME_KWARGS = Dict{Symbol,Any}()

# See @optimize_or_load for documentation –
# Only the macro version should be public!
function optimize_or_load(
    _filter,
    path,
    problem;
    method,
    filename::Union{Nothing,AbstractString}=nothing,
    suffix=((isnothing(filename) ? "jld2" : replace(splitext(filename)[2], "." => ""))),
    prefix="",
    tag::Bool=DrWatson.readenv("DRWATSON_TAG", DrWatson.istaggable(suffix)),
    gitpath=DrWatson.projectdir(),
    storepatch::Bool=DrWatson.readenv("DRWATSON_STOREPATCH", false),
    force=false,
    verbose=false,
    wsave_kwargs=Dict(),
    savename_kwargs=DEFAULT_OPTIMIZATION_SAVENAME_KWARGS,
    metadata=nothing,
    dry_run=false,
    kwargs...
)

    # _filter is only for attaching metadata to the result

    if isnothing(filename)
        filename = optimization_savename(
            problem;
            method=method,
            suffix=suffix,
            prefix=prefix,
            savename_kwargs=savename_kwargs,
            kwargs...
        )
    end

    file = joinpath(path, filename)
    if dry_run
        if verbose
            if isfile(file) && !force
                @info "Would load result from $file"
            else
                @info "Would optimize and store in $file"
            end
        end
        return nothing, file
    end
    if isfile(file) && verbose
        if force
            @info "Ignoring existing $file (force)"
        else
            @info "Loading result from $file"
        end
    end

    data, file = DrWatson.produce_or_load(
        path,
        Dict();
        filename=filename,
        suffix=suffix,
        prefix=prefix,
        tag=tag,
        gitpath=gitpath,
        loadfile=true,
        storepatch=storepatch,
        force=force,
        verbose=verbose,
        wsave_kwargs=wsave_kwargs,
        savename_kwargs...
    ) do _
        result = optimize(problem; method=method, verbose=verbose, kwargs...)
        data = Dict("result" => result)
        if !isnothing(_filter)
            data = _filter(data)
        end
        !isnothing(metadata) && merge!(data, metadata)
        return data
    end

    return data["result"], file

end

optimize_or_load(problem; kwargs...) = optimize_or_load(nothing, "", problem; kwargs...)
optimize_or_load(path, problem; kwargs...) =
    optimize_or_load(nothing, path, problem; kwargs...)

"""
Run [`optimize`](@ref) and store the result, or load the result if it exists.

```julia
result, file = @optimize_or_load(
    path="",
    problem;
    method=<method>,
    filename=nothing,
    suffix="jld2",
    prefix=DrWatson.default_prefix(config),
    tag=DrWatson.readenv("DRWATSON_TAG", true),
    gitpath=DrWatson.projectdir(),
    storepatch::Bool=DrWatson.readenv("DRWATSON_STOREPATCH", false),
    force=false,
    verbose=true,
    wsave_kwargs=Dict(),
    savename_kwargs=DEFAULT_OPTIMIZATION_SAVENAME_KWARGS,
    metadata=nothing,
    dry_run=false,
    kwargs...
)
```

runs `result = optimize(problem; method=<method>, kwargs...)` and stores
`result` in an automatically named file inside `path`. The
automatic file name is determined by [`optimization_savename`](@ref) and can be
overriden by passing an explicit `filename`. The full path to the output file
(`joinpath(path, filename)`) is returned as `file`.

In addition to the `result`, the data in the output `file` may also contain
some metadata, e.g. (automatically) "gitcommit" containing the git commit hash
of the project the produced the file, and "script" with the file name and line
number where `@optimize_or_load` was called, see [`load_optimization`](@ref).
If `metadata` is given as a dict on input, the data it contains will be
included in the output file.

If `file` already exists (and `force=false`), load the `result` from that file
instead of running the optimization.

If `dry_run=true`, return `(nothing, file)`. Depending on `verbose`, this will
print information about whether `file` would be loaded or generated.

The `@optimize_or_load` macro is intended to integrate well with the
[`DrWatson`](https://juliadynamics.github.io/DrWatson.jl/stable/) framework
for scientific projects and utilizes several configuration options and utility
functions from `DrWatson`, see below. Note that even though `DrWatson` is
recomended, you are not *required* to use if for your projects in order to use
`@optimize_or_load` or any other part of `QuantumControl`.

## I/O Keywords

The following keyword arguments determine where the result is stored and in
which format.

* `filename`: A file name to override the automatic file name. The `filename`
   should not contain slashes: use `path` for the folder where `filename`
   should be created.
* `suffix`, `prefix`, `savename_kwargs`: Parameters for
  [`optimization_savename`](@ref), which determines the automatic file name
* `tag`: Whether to record the current "gitcommit" as metadata alongside the
   optimization result, via
   [`DrWatson.tagsave`](https://juliadynamics.github.io/DrWatson.jl/stable/save/#DrWatson.tagsave).
   If not given explicitly, determine automatically from `suffix` or the
   extension of `filename`.
* `gitpath`, `storepatch`: Passed to `DrWatson.tagsave` if `tag` is `true`.
* `force`: If `true`, run and store the optimization regardless of whether
  `file` already exists.
* `verbose`: If `true`, print info about the process, if `file` does not exist.
* `wsave_kwargs`: Additional keyword arguments to pass to
  [`DrWatson.wsave`](https://juliadynamics.github.io/DrWatson.jl/stable/save/#Saving-Tools-1),
  e.g., to enable compression

All other keyword arguments are passed directly to [`optimize`](@ref).

## Related Functions

* [`optimization_savename`](@ref): Function that determines the automatic
  filename
* [`DrWatson.@produce_or_load`](https://juliadynamics.github.io/DrWatson.jl/stable/save/#DrWatson.@produce_or_load):
  The lower-level backend implementing the functionality of
  `@optimize_or_load`.
* [`load_optimization`](@ref): Function to load a file produced by
  `@optimize_or_load`
"""
macro optimize_or_load(path, problem, args...)
    args = Any[args...]
    # Keywords added after a `;` are moved to the front of the expression
    # that is passed to the macro. So instead of getting the path string
    # an Expr is passed.
    if path isa Expr && path.head == :parameters
        length(args) > 0 || return :(throw(
            MethodError(@optimize_or_load, $(esc(path)), $(esc(problem)), $(esc.(args)...)),
        ))
        extra_kw_def = path.args
        path = problem
        problem = popfirst!(args)
        append!(args, extra_kw_def)
    end
    # Save the source file name and line number of the calling line.
    s = QuoteNode(__source__)
    # Wrap the function f, such that the source can be saved in the data Dict.
    return quote
        optimize_or_load(
            $(esc(path)),
            $(esc(problem));
            $(esc.(DrWatson.convert_to_kw.(args))...)
        ) do data  # _filter
            # Extract the `gitpath` kw arg if it's there
            kws = ((; kwargs...) -> Dict(kwargs...))(
                $(esc.(DrWatson.convert_to_kw.(args))...),
            )
            gitpath = get(kws, :gitpath, DrWatson.projectdir())
            # Include the script tag with checking for the type of dict keys, etc.
            data = DrWatson.scripttag!(data, $s; gitpath=gitpath)
            return data
        end
    end
end

"""Determine an automatic filename for storing an optimization result.

```julia
file = optimization_savename(
    path="",
    problem;
    method=<method>,
    suffix="jld2",
    prefix="",
    savename_kwargs=DEFAULT_OPTIMIZATION_SAVENAME_KWARGS,
    kwargs...,
)
```

finds an appropriate automatic filename for the result of
`optimize(problem; method=<method>, kwargs...)`.

By default, the `file` has the structure
`<path>/<prefix>_<key1>=<value1>_..._<keyN>=<valueN>_method=<method>.jld2`
where the key-value pairs are a subset of the keyword arguments used to
instantiate `problem`, respectively the keyword arguments in `kwargs`. The
`prefix` is best used as a "name" for the optimization problem to ensure a
unique file name.

Which key-value pairs that are taken into account and the way they are
formatted can be customized via `savename_kwargs`.
See [`default_optimization_savename_kwargs`](@ref) for the supported options.
"""
function optimization_savename(
    path::AbstractString,
    problem::ControlProblem;
    method,
    suffix="jld2",
    prefix="",
    savename_kwargs=DEFAULT_OPTIMIZATION_SAVENAME_KWARGS,
    kwargs...
)
    for k ∈ keys(savename_kwargs)
        if k ∉ _SAVENAME_AVAILABLE_KEYS
            throw(
                ArgumentError(
                    "'$k' is not a valid keyword argument for savename. Use one of $(join(_SAVENAME_AVAILABLE_KEYS, ", "))",
                ),
            )
        end
    end
    c = Saving.OptimizationConfig(problem, method, kwargs)
    filename = DrWatson.savename(prefix, c, suffix; savename_kwargs...)
    return joinpath(path, filename)
end

optimization_savename(problem; kwargs...) = optimization_savename("", problem; kwargs...)

"""Set the default `savename_kwargs` for [`optimization_savename`](@ref).

```julia
savename_kwargs = default_optimization_savename_kwargs(;kwargs...)
```

sets entries in the `DEFAULT_OPTIMIZATION_SAVENAME_KWARGS` used in
`optimization_savename` and thus determines the automatic name used to store
optimization results.

Use

```julia
default_optimization_savename_kwargs(reset=true)
```

to clear the settings from any previous call to
`default_optimization_savename_kwargs`.

The following keyword arguments are supported, cf. [`DrWatson.savename`]
(https://juliadynamics.github.io/DrWatson.jl/dev/name/#DrWatson.savename):

* `accesses` - List of strings indicating which fields (keys in `kwargs` of
  [`ControlProblem`](@ref) or `kwargs` of
  [`optimize`](@ref)/[`@optimize_or_load`](@ref)) can be included in the
  output filename. By default, all fields with values matching `allowedtypes`
  (excluding anonymous functions) are used.
* `allowedtypes` - List of types of values eligible to be included in the
  filename. Defaults to `[Real, String, Symbol, TimeType, Function]`
* `connector` - String used to separate key-value pairs in the output filename.
  Defaults to `"_"`.
* `digits` - Used in `round` when formatting numbers, if no custom
  `val_to_string`.
* `equals` - String used between keys and values. Defaults to `"="`.
* `ignores` - List of strings indicating which fields should be ignored.
* `sigdigits` - Used in `round` when formatting numbers, if no custom
  `val_to_string`.
* `sort` - Whether to sort the fields alphabetically (default). If `false`, the
  resulting file name my not be stable.
* `val_to_string` - Function to convert values to string.
"""
function default_optimization_savename_kwargs(; reset=false, kwargs...)
    # Note: default values are effectively set by the methods defined in
    # src/saving.jl
    if reset
        empty!(DEFAULT_OPTIMIZATION_SAVENAME_KWARGS)
    end
    for (k, v) ∈ kwargs
        if k ∉ _SAVENAME_AVAILABLE_KEYS
            throw(
                ArgumentError(
                    "'$k' is not a valid keyword argument for savename. Use one of $(join(_SAVENAME_AVAILABLE_KEYS, ", "))",
                ),
            )
        end
        DEFAULT_OPTIMIZATION_SAVENAME_KWARGS[k] = v
    end
    return copy(DEFAULT_OPTIMIZATION_SAVENAME_KWARGS)
end

"""Load a previously stored optimization.

```julia
result = load_optimization(filename; verbose=true, kwargs...)
```

recovers a `result` previously stored by [`@optimize_or_load`](@ref).

```julia
result, metadata = load_optimization(filename; return_metadata=true, kwargs...)
```

also obtains a metadata dict containing e.g., "gitcommit" or "script" depending
on the options to [`@optimize_or_load`](@ref).

Calling `load_optimization` with `verbose=true` (default) will show the
metadata after loading the file.
"""
function load_optimization(file; return_metadata=false, verbose=true, kwargs...)
    data = DrWatson.wload(file)
    result = data["result"]
    metadata = filter(kv -> (kv[1] != "result"), data)
    if verbose
        metadata_str = join(["  $key: $val" for (key, val) ∈ metadata], "\n")
        @info ("Loaded optimization result from $file\n" * metadata_str)
    end
    if return_metadata
        return result, metadata
    else
        return result
    end
end


text:
import Base.copy

#*********************************#
#       COST FUNCTION CLASS       #
#*********************************#

abstract type CostFunction end
CostTrajectory = Vector{C} where C <: CostFunction

function cost(c::CostTrajectory, X::VectorTrajectory{T}, U::VectorTrajectory{T},dt::Vector{T})::T where T
    N = length(X)
    J = 0.0
    for k = 1:N-1
        J += stage_cost(c[k],X[k],U[k],dt[k])
    end
    J += stage_cost(c[N],X[N])
    return J
end

"$(TYPEDEF) Expansion of cost function"
struct Expansion{T<:AbstractFloat,Q<:AbstractMatrix,R<:AbstractMatrix}
    x::Vector{T}
    u::Vector{T}
    xx::Q
    uu::R
    ux::Matrix{T}
end

function Expansion{T,Q,R}(n::Int, m::Int) where {T,Q,R}
    x = zeros(T,n)
    u = zeros(T,m)
    xx = Q(zeros(T,n,n))
    uu = R(zeros(T,m,m))
    ux = zeros(T,m,n)
    Expansion(x, u, xx, uu, ux)
end
Expansion{T}(n::Int, m::Int) where T = Expansion{T,Matrix{T},Matrix{T}}(n,m)

import Base./, Base.*
function *(e::Expansion, a::Real)
    e.x .*= a
    e.u .*= a
    e.xx .*= a
    e.uu .*= a
    e.ux .*= a
    return nothing
end
*(a::Real, e::Expansion) = e*a

function /(e::Expansion,a::Real)
    e.x ./= a
    e.u ./= a
    e.xx ./= a
    e.uu ./= a
    e.ux ./= a
    return nothing
end

function copy(e::Expansion)
    Expansion(copy(e.x),copy(e.u),copy(e.xx),copy(e.uu),copy(e.ux))
end

function reset!(e::Expansion)
    !isempty(e.x) ? e.x .= zero(e.x) : nothing
    !isempty(e.u) ? e.u .= zero(e.u) : nothing
    !isempty(e.xx) ? e.xx .= zero(e.xx) : nothing
    !isempty(e.uu) ? e.uu .= zero(e.uu) : nothing
    !isempty(e.ux) ? e.ux .= zero(e.ux) : nothing
    return nothing
end

ExpansionTrajectory{T} = Vector{Expansion{T}} where T <: AbstractFloat

function reset!(et::ExpansionTrajectory)
    for e in et
        reset!(e)
    end
end

#######################################################
#              COST FUNCTION INTERFACE                #
#######################################################

#TODO this documentation breaks functionality by create ambiguous methods...

# "$(SIGNATURES) Evaluate the cost at state `x` and control `u`"
# stage_cost(cost::CostFunction, x, u) = 0.0
# "$(SIGNATURES) Evaluate the cost at the terminal state `xN`"
# stage_cost(cost::CostFunction, xN) = 0.0
# "$(SIGNATURES) Evaluate the second order expansion at state `x` and control `u`"
# cost_expansion!(Q::Expansion, cost::CostFunction, x, u) = nothing
# "$(SIGNATURES) Evaluate the second order expansion at the terminal state `xN`"
# cost_expansion!(Q::Expansion, cost::CostFunction, xN) = nothing

"""
$(TYPEDEF)
Cost function of the form
    1/2xₙᵀ Qf xₙ + qfᵀxₙ +  ∫ ( 1/2xᵀQx + 1/2uᵀRu + xᵀHu + qᵀx  rᵀu ) dt from 0 to tf
R must be positive definite, Q and Qf must be positive semidefinite

Constructor use any of the following constructors:
```julia
QuadraticCost(Q, R, H, q, r, c)
QuadraticCost(Q, R; H, q, r, c)
QuadraticCost(Q, q, c)
```
Any optional or omitted values will be set to zero(s).
"""
mutable struct QuadraticCost{T} <: CostFunction
    Q::AbstractMatrix{T}                 # Quadratic stage cost for states (n,n)
    R::AbstractMatrix{T}                 # Quadratic stage cost for controls (m,m)
    H::AbstractMatrix{T}                 # Quadratic Cross-coupling for state and controls (n,m)
    q::AbstractVector{T}                 # Linear term on states (n,)
    r::AbstractVector{T}                 # Linear term on controls (m,)
    c::T                                 # constant term
    function QuadraticCost(Q::AbstractMatrix{T}, R::AbstractMatrix{T}, H::AbstractMatrix{T},
            q::AbstractVector{T}, r::AbstractVector{T}, c::T) where T
        if !isposdef(R)
            @warn "R is not positive definite"
        end
        if !ispossemidef(Q)
            err = ArgumentError("Q must be positive semi-definite")
            throw(err)
        end

        new{T}(Q,R,H,q,r,c)
    end
end

get_sizes(cost::QuadraticCost) = (size(cost.Q,1),size(cost.R,1))

function QuadraticCost(Q,R; H=zeros(size(R,1), size(Q,1)), q=zeros(size(Q,1)),
        r=zeros(size(R,1)), c=0.0)
    QuadraticCost(Q,R,H,q,r,c)
end

function QuadraticCost(Q,q,c)
    QuadraticCost(Q,zeros(0,0),zeros(0,size(Q,1)),q,zeros(0),c)
end

"""
$(SIGNATURES)
Cost function of the form
``(x-x_f)^T Q (x_x_f) + u^T R u``
R must be positive definite, Q must be positive semidefinite
"""
function LQRCost(Q::AbstractArray, R::AbstractArray, xf::AbstractVector)
    H = zeros(size(R,1),size(Q,1))
    q = -Q*xf
    r = zeros(size(R,1))
    c = 0.5*xf'*Q*xf
    return QuadraticCost(Q, R, H, q, r, c)
end

"""
$(SIGNATURES)
Cost function of the form
``(x-x_f)^T Q (x_x_f)``
Q must be positive semidefinite
"""
function LQRCostTerminal(Qf::AbstractArray,xf::AbstractVector)
    qf = -Qf*xf
    cf = 0.5*xf'*Qf*xf
    return QuadraticCost(Qf,zeros(0,0),zeros(0,size(Qf,1)),qf,zeros(0),cf)
end

function stage_cost(cost::QuadraticCost, x::AbstractVector, u::AbstractVector)
    0.5*x'cost.Q*x + 0.5*u'*cost.R*u + cost.q'x + cost.r'u + cost.c + u'*cost.H*x
end

function stage_cost(cost::QuadraticCost, x::AbstractVector, u::AbstractVector, dt)
    (0.5*x'cost.Q*x + 0.5*u'*cost.R*u + cost.q'x + cost.r'u + cost.c + u'*cost.H*x)*dt
end

function stage_cost(cost::QuadraticCost, xN::AbstractVector{T}) where T
    0.5*xN'cost.Q*xN + cost.q'*xN + cost.c
end

function cost_expansion!(Q::Expansion{T}, cost::QuadraticCost,
        x::AbstractVector{T}, u::AbstractVector{T}, dt::T) where T
    Q.x .= cost.Q*x + cost.q + cost.H'*u
    Q.u .= cost.R*u + cost.r + cost.H*x
    Q.xx .= cost.Q
    Q.uu .= cost.R
    Q.ux .= cost.H
    Q*dt
    return nothing
end

function cost_expansion!(S::Expansion{T}, cost::QuadraticCost, xN::AbstractVector{T}) where T
    S.xx .= cost.Q
    S.x .= cost.Q*xN + cost.q
    return nothing
end

function gradient!(grad, cost::QuadraticCost,
        x::AbstractVector, u::AbstractVector,dt)
    grad.x .= (cost.Q*x + cost.q + cost.H'*u)
    grad.u .= (cost.R*u + cost.r + cost.H*x)

    grad .*= dt
    return nothing
end

function gradient!(grad, cost::QuadraticCost, xN::AbstractVector)
    grad .= cost.Q*xN + cost.q
    return nothing
end

function hessian!(hess, cost::QuadraticCost,
        x::AbstractVector, u::AbstractVector, dt)
    hess.xx .= cost.Q
    hess.uu .= cost.R
    hess.ux .= cost.H
    hess.xu .= cost.H'

    hess .*= dt
    return nothing
end

function hessian!(hess, cost::QuadraticCost, xN::AbstractVector)
    hess .= cost.Q
    return nothing
end

function copy(cost::QuadraticCost)
    return QuadraticCost(copy(cost.Q), copy(cost.R), copy(cost.H), copy(cost.q), copy(cost.r), copy(cost.c))
end

"""
$(TYPEDEF)
Cost function of the form
    ℓf(xₙ) + ∫ ℓ(x,u) dt from 0 to tf
"""
struct GenericCost <: CostFunction
    ℓ::Function             # Stage cost
    ℓf::Function            # Terminal cost
    expansion::Function     # 2nd order Taylor Series Expansion of the form,  Q,R,H,q,r = expansion(x,u)
    n::Int                  #                                                     Qf,qf = expansion(xN)
    m::Int
end

"""
$(SIGNATURES)
Create a Generic Cost, specifying the gradient and hessian of the cost function analytically

# Arguments
* hess: multiple-dispatch function of the form,
    Q,R,H = hess(x,u) with sizes (n,n), (m,m), (m,n)
    Qf = hess(xN) with size (n,n)
* grad: multiple-dispatch function of the form,
    q,r = grad(x,u) with sizes (n,), (m,)
    qf = grad(x,u) with size (n,)

"""
function GenericCost(ℓ::Function, ℓf::Function, grad::Function, hess::Function, n::Int, m::Int)
    function expansion(x::Vector{T},u::Vector{T}) where T
        Q,R,H = hess(x,u)
        q,r = grad(x,u)
        return Q,R,H,q,r
    end
    expansion(xN) = hess(xN), grad(xN)
    GenericCost(ℓ,ℓf, expansion, n,m)
end

"""
$(SIGNATURES)
Create a Generic Cost. Gradient and Hessian information will be determined using ForwardDiff

# Arguments
* ℓ: stage cost function of the form J = ℓ(x,u)
* ℓf: terminal cost function of the form J = ℓ(xN)
"""
function GenericCost(ℓ::Function, ℓf::Function, n::Int, m::Int)
    linds = LinearIndices(zeros(n+m,n+m))
    xinds = 1:n
    uinds = n .+(1:m)
    inds = (x=xinds, u=uinds, xx=linds[xinds,xinds], uu=linds[uinds,uinds], ux=linds[uinds,xinds])
    expansion = auto_expansion_function(ℓ,ℓf,n,m)

    GenericCost(ℓ,ℓf, expansion, n,m)
end

function auto_expansion_function(ℓ::Function,ℓf::Function,n::Int,m::Int)
    z = zeros(n+m)
    hess = zeros(n+m,n+m)
    grad = zeros(n+m)
    qf,Qf = zeros(n), zeros(n,n)

    linds = LinearIndices(hess)
    xinds = 1:n
    uinds = n .+(1:m)
    inds = (x=xinds, u=uinds, xx=linds[xinds,xinds], uu=linds[uinds,uinds], ux=linds[uinds,xinds])
    function ℓ_aug(z::Vector{T}) where T
        x = view(z,xinds)
        u = view(z,uinds)
        ℓ(x,u)
    end
    function expansion(x::Vector{T},u::Vector{T}) where T
        z[inds.x] = x
        z[inds.u] = u
        ForwardDiff.hessian!(hess, ℓ_aug, z)
        Q = view(hess,inds.xx)
        R = view(hess,inds.uu)
        H = view(hess,inds.ux)

        ForwardDiff.gradient!(grad, ℓ_aug, z)
        q = view(grad,inds.x)
        r = view(grad,inds.u)
        return Q,R,H,q,r
    end
    function expansion(xN::Vector{T}) where T
        ForwardDiff.gradient!(qf,ℓf,xN)
        ForwardDiff.hessian!(Qf,ℓf,xN)
        return Qf, qf
    end
end

stage_cost(cost::GenericCost, x::Vector{T}, u::Vector{T}) where T = cost.ℓ(x,u)
stage_cost(cost::GenericCost, xN::Vector{T}) where T = cost.ℓf(xN)

function cost_expansion!(Q::Expansion{T}, cost::GenericCost, x::Vector{T},
        u::Vector{T}) where T

    e = cost.expansion(x,u)

    Q.x .= e[4]
    Q.u .= e[5]
    Q.xx .= e[1]
    Q.uu .= e[2]
    Q.ux .= e[3]
    return nothing
end

function cost_expansion!(S::Expansion{T}, cost::GenericCost, xN::Vector{T}) where T
    Qf, qf = cost.expansion(xN)
    S.xx .= Qf
    S.x .= qf
    return nothing
end

copy(cost::GenericCost) = GenericCost(cost.ℓ,cost.ℓ,cost.n,cost.m)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
push!(LOAD_PATH,joinpath(@__DIR__, "../src/"))
using Documenter, StableApproxEPW

# Generating documentation with `Documenter.jl`
makedocs(
    modules = [StableApproxEPW],
    clean = false,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename="Stable Approx with EPW",
    authors="Emile Parolin",
    pages=[
        "Getting started" => "index.md",
        "Examples" => "example.md",
        "Source" => "implementation.md",
    ],
)

deploydocs(
    repo = "github.com/EmileParolin/evanescent-plane-wave-approx.git",
)

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
import Benchmark
Benchmark.run_benchmark(
    num_facilities = 100,
    num_customers = 100,
    num_locations = 100,
    time_limit_sec = 0.001,
    max_iters = 1,
)

max_stars_repo_name: jump-dev/MOIPaperBenchmarks

max_issues_repo_name: jump-dev/MOIPaperBenchmarks

max_forks_repo_name: jump-dev/MOIPaperBenchmarks


text:
using Statistics
using Plots
using BenchmarkTools
using Profile
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using Random
using Distributions
using StatsBase
using FITSIO
using SparseArrays
using MultivariateStats
using FileIO

push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils
#include("Deriv_Utils_New.jl")
#import Deriv_Utils_New
push!(LOAD_PATH, pwd()*"/scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs

##How Different are the dbns of SFD images, SFD+white_noise, SFD+white_noise+smoothing
Nx=128
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf, ) = size(filter_hash["filt_index"])
coeff_mask = falses(2+Nf+Nf^2)
coeff_mask[Nf+3:end] .= true
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true) #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense

#In S20
sfdall = readsfd(Nx)
Num=2000
sfdsub = sfdall[:, :, 1:Num]
oridbn = get_dbn_coeffs(sfdsub, filter_hash, dhc_args, coeff_mask)
noisydbn = zeros(Num, size(oridbn)[2])
smoothdbn1 = zeros(Num, size(oridbn)[2])
noisysub = zeros(Nx, Nx, Num)
smoothsub = zeros(Nx, Nx, Num)
for s=1:Num
  tim = sfdsub[:,:, s]
  noise = reshape(rand(Normal(0.0, std(tim)/5),Nx^2), (Nx, Nx))
  init = tim+noise
  noisysub[:, :, s] .= init
  smoothsub[:, :, s] .= imfilter(init, Kernel.gaussian(1.0))
  noisydbn[s, :] .= DHC_compute_apd(init, filter_hash; dhc_args...)
  smoothdbn1[s, :] .= DHC_compute_apd(imfilter(init, Kernel.gaussian(1.0)), filter_hash; dhc_args...)
end

#sfdvecv = collect(eachslice(sfdsub, dims=3))
#apdsfd = mapslices(apodizer, sfdsub, dims=[1, 2])
#apddbn = get_dbn_coeffs(apdsfd, filter_hash, dhc_args, coeff_mask)
#Skip the below section

##3 dbns
#=
alldbns = vcat(oridbn, noisydbn, apddbn)
M = fit(PCA, alldbns'; maxoutdim=3)
s20_pca = transform(M, alldbns')

p = scatter(s20_pca[1,1:Num],s20_pca[2,1:Num], marker=:circle,linewidth=0,label="Ori")
scatter!(s20_pca[1, (2*Num)+1:end],s20_pca[2, (2*Num)+1:end], marker=:circle,linewidth=0,label="Apd")
#scatter!(s20_pca[1, Num+1:Num*2],s20_pca[2, Num+1:Num*2], marker=:circle,linewidth=0,label="Smoothed Noisy")
scatter!(s20_pca[1, Num+1:Num*2],s20_pca[2, Num+1:Num*2], marker=:circle,linewidth=0,label="Noisy")
#plot!(xlims=(-1e4, 1e3), ylims=(-100, 1e3))
plot!(p,xlabel="PC1",ylabel="PC2")
plot!(title="Reduction computed over all 3 distributions")

#4 dbns
alldbns = vcat(oridbn, apddbn, smoothdbn1, noisydbn)'
M = fit(PCA, alldbns; maxoutdim=3)
s20_pca = transform(M, alldbns)

p = scatter(s20_pca[1,1:Num],s20_pca[2,1:Num], marker=:circle,linewidth=0,label="Ori")
scatter!(s20_pca[1, Num+1:2*Num],s20_pca[2, Num+1:2*Num], marker=:circle,linewidth=0,label="Apd")
scatter!(s20_pca[1, (2*Num)+1:Num*3],s20_pca[2, (2*Num)+1:Num*3], marker=:circle,linewidth=0,label="Smoothed Noisy")
scatter!(s20_pca[1, (3*Num)+1:end],s20_pca[2, (3*Num)+1:end], marker=:circle,linewidth=0,label="Noisy")
#plot!(xlims=(-1e4, 1e3), ylims=(-100, 1e3))
plot!(p,xlabel="PC1",ylabel="PC2")
plot!(p, legend=:bottomleft)
plot!(title="Reduction computed over all 4 distributions")

p=scatter(s20_pca[1, (2*Num)+1:Num*3],s20_pca[2, (2*Num)+1:Num*3], marker=:circle,linewidth=0,label="Smoothed Noisy")
scatter!(s20_pca[1, (3*Num)+1:end],s20_pca[2, (3*Num)+1:end], marker=:circle,linewidth=0,label="Noisy")
#plot!(xlims=(-1e4, 1e3), ylims=(-100, 1e3))
plot!(p,xlabel="PC1",ylabel="PC2")
plot!(p, legend=:bottomleft)
plot!(title="Reduction computed over all 4 distributions: plotting only 2")
=#
##4 dbns
alldbns = vcat(oridbn, smoothdbn1, noisydbn)'

function pca_norm(inputdbn)
  #inputdbn: Rows are features, Cols are observations
  normed = (inputdbn .- mean(inputdbn, dims=2))./std(inputdbn, dims=2)
  M = fit(PCA, normed; maxoutdim=3)
  return mean(inputdbn, dims=2), std(inputdbn, dims=2), M, transform(M, normed)
end

dbn_mean, dbn_std, M, s20_pca = pca_norm(alldbns)

p = scatter(s20_pca[1,1:Num],s20_pca[2,1:Num], marker=:circle,linewidth=0,label="Ori")
scatter!(s20_pca[1, Num+1:2*Num],s20_pca[2, Num+1:2*Num], marker=:circle,linewidth=0,label="Smoothed Noisy")
scatter!(s20_pca[1, (2*Num)+1:Num*3],s20_pca[2, (2*Num)+1:Num*3], marker=:circle,linewidth=0,label="Noisy")
#scatter!(s20_pca[1, (3*Num)+1:end],s20_pca[2, (3*Num)+1:end], marker=:circle,linewidth=0,label="Noisy")
#plot!(xlims=(-1e4, 1e3), ylims=(-100, 1e3))
plot!(p,xlabel="PC1",ylabel="PC2")
plot!(p, legend=:bottomleft)
plot!(title="S20 PCA of apodized images")
plot!(xlims=(-50, 50), ylims=(-50, 50))
savefig(p, "scratch_NM/PCA/sfd_apd_nolog_std20per.png")
#very little variance along PC3
#=
p = scatter(s20_pca[1,1:Num],s20_pca[2,1:Num],s20_pca[3,1:Num], marker=:circle,linewidth=0,label="Ori")
scatter!(s20_pca[1, (2*Num)+1:end],s20_pca[2, (2*Num)+1:end],s20_pca[3, (2*Num)+1:end], marker=:circle,linewidth=0,label="Apd")
scatter!(s20_pca[1, Num+1:Num*2],s20_pca[2, Num+1:Num*2],s20_pca[3, Num+1:Num*2], marker=:circle,linewidth=0,label="Noisy")
plot!(xlims=(-1e4, 1e3), ylims=(-100, 1e3))
plot!(p,xlabel="PC1",ylabel="PC2", zlabel="PC3")
plot!(title="Reduction computed over all 3 distributions")
=#

##S20+Log

ltz = findall(noisysub.<0)
noisysub_bd = copy(noisysub)
noisysub_bd[ltz] .= 1e-6

ltz = findall(smoothsub.<0)
smoothsub_bd = copy(smoothsub)
smoothsub_bd[ltz] .= 1e-6

loridbn = get_dbn_coeffs(log.(sfdsub), filter_hash, dhc_args, coeff_mask)
lnoisydbn = get_dbn_coeffs(log.(noisysub_bd), filter_hash, dhc_args, coeff_mask)
lsmoothdbn1 = get_dbn_coeffs(log.(smoothsub_bd), filter_hash, dhc_args, coeff_mask)

alldbns_log = vcat(loridbn, lsmoothdbn1, lnoisydbn)'
ldbn_mean, ldbn_std, M, ls20_pca = pca_norm(alldbns_log)

p = scatter(ls20_pca[1,1:Num],ls20_pca[2,1:Num], marker=:circle,linewidth=0,label="Ori")
scatter!(ls20_pca[1, Num+1:2*Num],ls20_pca[2, Num+1:2*Num], marker=:circle,linewidth=0,label="Smoothed Noisy")
scatter!(ls20_pca[1, (2*Num)+1:Num*3],ls20_pca[2, (2*Num)+1:Num*3], marker=:circle,linewidth=0,label="Noisy")
#plot!(xlims=(-1e4, 1e3), ylims=(-100, 1e3))
plot!(p,xlabel="PC1",ylabel="PC2")
plot!(p, legend=:bottomleft)
plot!(title="S20 PCA of Apodized Log Images")
plot!(xlims=(-50, 50), ylims=(-50, 50))
savefig(p, "scratch_NM/PCA/sfd_apd_log.png")

##Considering ONLY S20
coeff_mask = falses(2+Nf+Nf^2)
coeff_mask[Nf+3:end] .= true
alldbns20 = alldbns[coeff_mask, :]

dbn_mean, dbn_std, M, s20_pca = pca_norm(alldbns20)
p = scatter(s20_pca[1,1:Num],s20_pca[2,1:Num], marker=:circle,linewidth=0,label="Ori")
scatter!(s20_pca[1, Num+1:2*Num],s20_pca[2, Num+1:2*Num], marker=:circle,linewidth=0,label="Smoothed Noisy")
scatter!(s20_pca[1, (2*Num)+1:Num*3],s20_pca[2, (2*Num)+1:Num*3], marker=:circle,linewidth=0,label="Noisy")
#scatter!(s20_pca[1, (3*Num)+1:end],s20_pca[2, (3*Num)+1:end], marker=:circle,linewidth=0,label="Noisy")
#plot!(xlims=(-1e4, 1e3), ylims=(-100, 1e3))
plot!(p,xlabel="PC1",ylabel="PC2")
plot!(p, legend=:bottomleft)
plot!(title="S20 PCA of apodized images")
plot!(xlims=(-50, 50), ylims=(-50, 50))
savefig(p, "scratch_NM/PCA/sfd_apd_nolog_std20per.png")

#Log
alldbns_logs20 = alldbns_log[coeff_mask, :]
ldbn_mean, ldbn_std, M, ls20_pca = pca_norm(alldbns_logs20)
p = scatter(ls20_pca[1,1:Num],ls20_pca[2,1:Num], marker=:circle,linewidth=0,label="Ori")
scatter!(ls20_pca[1, Num+1:2*Num],ls20_pca[2, Num+1:2*Num], marker=:circle,linewidth=0,label="Smoothed Noisy")
scatter!(ls20_pca[1, (2*Num)+1:Num*3],ls20_pca[2, (2*Num)+1:Num*3], marker=:circle,linewidth=0,label="Noisy")
#plot!(xlims=(-1e4, 1e3), ylims=(-100, 1e3))
plot!(p,xlabel="PC1",ylabel="PC2")
plot!(p, legend=:bottomleft)
plot!(title="S20 PCA of Apodized Log Images")
plot!(xlims=(-50, 50), ylims=(-50, 50))
savefig(p, "scratch_NM/PCA/sfd_apd_log.png")

##Reconstruction
Nx=64
im = readsfd(Nx)
true_img = im[:, :, 1]

noise = rand(Normal(0.0, std(true_img)/5.0), Nx, Nx)
init = imfilter(true_img + noise, Kernel.gaussian(1.0))
#heatmap(init)
#heatmap(true_img)
#PCA
#ltz = findall(init .< 0)
#init_bd = copy(init)
#init_bd[ltz] .= 1e-6

#trues = DHC_compute_apd(log.(true_img), filter_hash; dhc_args...)[coeff_mask]
#inits = DHC_compute_apd(log.(init_bd), filter_hash; dhc_args...)[coeff_mask]

#trues_norm = (trues .- ldbn_mean)./ldbn_std
#inits_norm = (inits .- ldbn_mean)./ldbn_std
#trues_pca = transform(M, trues_norm)
#inits_pca = transform(M, inits_norm)

filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf, ) = size(filter_hash["filt_index"])
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true) #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
coeff_mask = falses(2+Nf+Nf^2)
coeff_mask[Nf+3:end] .= true

white_noise_args = Dict(:loc=>0.0, :sig=>std(true_img)/5.0, :Nsam=>1000, :iso=>false, :norm=>false, :smooth=>true, :smoothval=>1.0) #only if you're using noise based covar
optim_settings = Dict([("iterations", 1000), ("norm", false)])
recon_settings = Dict([("target_type", "ground_truth"), ("covar_type", "white_noise"), ("log", true), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", std(true_img).^(-2)), ("white_noise_args", white_noise_args)]) #Add constraints
#"white_noise_args": args for white_noise sims if GaussianLoss and Covar_type='white_noise'. :loc...

s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings)
recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov
fname_save = "scratch_NM/PCA/log_apd_gt.jld2"
recon_settings["safemode"] = true
recon_settings["fname_save"] = fname_save
recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)

Visualization.plot_diffscales([true_img, init, recon_img, recon_img - true_img], ["GT", "Init", "Reconstruction", "Residual"], fname="scratch_NM/NewWrapper/3-24/kernel1smoothing/sfd_dbntargwhitenoisecov_lam0-1_smoothinit_apd_logspace.png")

println("Mean Abs Res", mean(abs.(init - true_img)), " ", mean(abs.(recon_img - true_img)))

aptrue = apodizer(true_img)
apinit = apodizer(init)
aprecon = apodizer(recon_img)
Visualization.plot_diffscales([apodizer(true_img), apodizer(init), apodizer(recon_img), apodizer(recon_img) - apodizer(true_img)], ["GT", "Init", "Reconstruction", "Residual"], fname="scratch_NM/NewWrapper/3-24/kernel1smoothing/APDsfd_dbntargwhitenoisecov_lam0-1_smoothinit_apd_logspacekern1.png")
Visualization.plot_synth_QA(apodizer(true_img), apodizer(init), apodizer(recon_img), fname="scratch_NM/NewWrapper/3-24/kernel1smoothing/APDsfd_dbntargwhitenoisecov_lam0-1_smoothinit_apd_logspace_kern1_6p.png")
println("Mean Abs Res: Init = ", mean(abs.(apinit - aptrue)), " Recon = ", mean(abs.(aprecon - aptrue)))

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

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
"""
    with_preview(f, on_preview, executor, interval::Real, [throttle]) -> result

Create a `previewer` and call `f(previewer)` whose `result` is returned as-is.

The `previewer` object passed to `f` is a transducer that acts like an identity
transducer (i.e., `Map(identity)`) with respect to the final result. However,
in another "background" task, it also calls the inner reducing function to
combine and complete the *intermediate* result and then passed it to
`on_preview`. The interval of calling `on_preview` can be configured by
`interval` and `throttle`.
"""
with_preview

"""
    with_distributed_preview(f, on_preview, interval, [throttle]) -> result

Equivalent to [`with_preview`](@ref) called with `DistributedEx` executor.
"""
with_distributed_preview

"""
    with_threaded_preview(f, on_preview, interval, [throttle]) -> result

Equivalent to [`with_preview`](@ref) called with `ThreadedEx` executor.
"""
with_threaded_preview

struct BasecaseId
    pid::typeof(myid())
    lid::UInt
end

const PROCESS_LOCAL_ID = Threads.Atomic{UInt}(0)

BasecaseId() = BasecaseId(myid(), Threads.atomic_add!(PROCESS_LOCAL_ID, UInt(1)))

remote_channel(size = 0) = RemoteChannel(() -> Channel(size))

channel_for(f, args...) = f(args...)
channel_for(::Executor, args...) = Channel(args...)
channel_for(::DistributedEx, args...) = remote_channel(args...)

struct Preview{C,S,T<:Throttle} <: Transducer
    channel::C
    start_channel::S
    throttle::T
end

with_distributed_preview(
    f,
    on_preview,
    interval::Real,
    thw::Throttle = IntervalThrottle(interval),
) = with_preview(f, on_preview, remote_channel, interval, thw)

with_threaded_preview(
    f,
    on_preview,
    thw::Throttle,
    interval::Real = IntervalThrottle(interval),
) = with_preview(f, on_preview, Channel, thw, interval)

function with_preview(f, on_preview, make_channel, th::IntervalThrottle)
    @warn(
        "DEPRECATED: Please use `with_preview(f, on_preview, executor, interval)`" *
        " instead of `with_preview(f, on_preview, executor, IntervalThrottle(interval))`",
        maxlog = 2,
    )
    return with_preview(f, on_preview, make_channel, th.interval)
end

function with_preview(
    f,
    on_preview,
    make_channel,
    interval::Real,
    thw::Throttle = IntervalThrottle(interval),
)
    interval >= 0 || error("`interval` must be positive; got: interval = $interval")
    # TODO: better buffer size detection
    channel = channel_for(make_channel, Threads.nthreads() * nprocs())
    start_channel = channel_for(make_channel, 0)
    previewer = Preview(channel, start_channel, thw)
    accref = LockedRef{Any}()
    isdone() = !isopen(channel)
    @sync try
        @async try
            updater_loop(accref, previewer)
        catch err
            @debug "`updater_loop` task" exception = (err, catch_backtrace())
            rethrow()
        finally
            close(channel)
            close(start_channel)
            close(accref)
        end
        @async try
            preview_loop(on_preview, interval, accref, isdone)
        catch err
            @debug "`preview_loop` task" exception = (err, catch_backtrace())
            rethrow()
        finally
            close(channel)
            close(start_channel)
            close(accref)
        end
        @syncthrow f(previewer)
    finally
        close(channel)
        close(start_channel)
        close(accref)
    end
end

# A very hacky solution for using the reducing function in `preview_loop`.
# Ideally, we can do something like this via fold initialization protocol.
function send_start!(start_channel, starter)
    try
        put!(start_channel, starter)
    catch err
        # Errors here are expected here since `start_channel` will be closed
        # after the first element is taken.
        # @debug "`send_start!`" exception = (err, catch_backtrace())
    end
    return
end

function preview_loop(on_preview, interval, accref, isdone)
    acc0 = tryfetch(accref)
    acc0 === nothing && return
    function preview_loop_impl()
        on_preview(something(acc0))
        while true
            # @info "`preview_loop`: Sleeping $interval seconds..."
            sleep(interval)
            isdone() && return
            local acc = tryfetch(accref)
            acc === nothing && return
            on_preview(something(acc))
        end
    end
    Base.invokelatest(preview_loop_impl)
end

function updater_loop(accref, previewer)
    rf, init = try
        take!(previewer.start_channel)
    catch err
        @debug "`updater_loop`: on `take!(start_channel)`" exception =
            (err, catch_backtrace())
        return
    end
    close(previewer.start_channel)
    Base.invokelatest(updater_loop, accref, previewer, rf, init)
end

function updater_loop(accref, previewer, rf, init)
    channel = previewer.channel
    accs = Dict{BasecaseId,Any}()
    while true
        bid, a = try
            take!(channel)
        catch err
            @debug "`updater_loop`: on `take!(channel)`" exception =
                (err, catch_backtrace())
            return
        end
        accs[bid] = a

        # TODO: use monoid cached tree
        acc = start(rf, init)
        for x in values(accs)
            acc = combine(rf, acc, x)
        end
        tryset!(accref, complete(rf, acc)) || break
    end
end

function Transducers.start(rf::R_{Preview}, init)
    iacc = start(inner(rf), init)
    bid = BasecaseId()
    tstate = Throttles.init(xform(rf).throttle)
    return wrap(rf, (bid, tstate, false, init), iacc)
end

Transducers.next(rf::R_{Preview}, acc, input) =
    wrapping(rf, acc) do (bid, tstate, started, init), iacc
        iacc = next(inner(rf), iacc, input)
        go, tstate = Throttles.step(xform(rf).throttle, tstate)
        if go
            if !started
                send_start!(xform(rf).start_channel, (inner(rf), init))
            end
            started = true
            put!(xform(rf).channel, (bid, iacc))
        end
        return (bid, tstate, started, init), iacc
    end
# Note: Postpone `send_start!` until we hit the first `go`. Otherwise, with
# `ThreadedEx` and alike, each base case `Task` will hit the yield point right
# after it is started. It makes the scheduler to (likely) re-use the OS thread
# (which is very bad for performance).

Transducers.complete(rf::R_{Preview}, acc) = complete(inner(rf), last(unwrap(rf, acc)))

function Transducers.combine(rf::R_{Preview}, a, b)
    (bid, tstatea, starteda, init), ira = unwrap(rf, a)
    (_, tstateb, startedb), irb = unwrap(rf, b)
    irc = combine(inner(rf), ira, irb)
    go, tstatea = Throttles.step(xform(rf).throttle, tstatea)
    if !go
        go, _ = Throttles.step(xform(rf).throttle, tstateb)
    end
    started = starteda || startedb
    if go
        if started
            send_start!(xform(rf).start_channel, (inner(rf), init))
        end
        started = true
        put!(xform(rf).channel, (bid, irc))
    end
    return wrap(rf, (bid, tstatea, started, init), irc)
end


text:
# mathematical constants related to statistics

@irrational twoπ   6.2831853071795864769 2 * big(π)
@irrational fourπ  12.566370614359172954 4 * big(π)
@irrational halfπ  1.5707963267948966192 big(π) / 2
@irrational quartπ 0.7853981633974483096 big(π) / 4

@irrational invπ     0.31830988618379067154 inv(big(π))
@irrational twoinvπ  0.63661977236758134308 2 / big(π)
@irrational fourinvπ 1.27323954473516268615 4 / big(π)
@irrational inv2π    0.159154943091895335769 inv(2 * big(π))
@irrational inv4π    0.079577471545947667884 inv(4 * big(π))

@irrational sqrt2     1.4142135623730950488 sqrt(big(2))
@irrational sqrt3     1.7320508075688772935 sqrt(big(3))
@irrational sqrtπ     1.7724538509055160273 sqrt(big(π))
@irrational sqrt2π    2.5066282746310005024 sqrt(2 * big(π))
@irrational sqrt4π    3.5449077018110320546 sqrt(4 * big(π))
@irrational sqrthalfπ 1.2533141373155002512 sqrt(big(π) / 2)

@irrational invsqrt2  0.7071067811865475244 inv(sqrt(big(2)))
@irrational invsqrtπ  0.5641895835477563 inv(sqrt(big(π)))
@irrational invsqrt2π 0.3989422804014326779 inv(sqrt(2 * big(π)))

@irrational loghalf -0.6931471805599453094 log(inv(big(2)))
@irrational logtwo  0.6931471805599453094 log2
@irrational logten  2.302585092994046 log(big(10))
@irrational logπ    1.1447298858494001741 log(big(π))
@irrational log2π   1.8378770664093454836 log(2 * big(π))
@irrational log4π   2.5310242469692907930 log(4 * big(π))