text:
u = [0.4, 0.26]

@testset "PolynomialLimbDark" begin
    ld = PolynomialLimbDark(u)

    @test sprint(show, ld) == "PolynomialLimbDark([-1.0, 0.4, 0.26])"

    @test repr("text/plain", ld) == """
    PolynomialLimbDark
     u_n: [-1.0, 0.4, 0.26]"""

    @test repr("text/html", ld) == """
    PolynomialLimbDark
     un: [-1.0, 0.4, 0.26]"""
end

@testset "QuadLimbDark" begin
    ld = QuadLimbDark(u)

    @test sprint(show, ld) == "QuadLimbDark([-1.0, 0.4, 0.26])"

    @test sprint(show, "text/plain", ld) == """
    QuadLimbDark
     u_n: [-1.0, 0.4, 0.26]"""

    @test sprint(show, "text/html", ld) == """
    QuadLimbDark
     un: [-1.0, 0.4, 0.26]"""
end

@testset "SecondaryLimbDark" begin
    ld = SecondaryLimbDark(u)

    @test sprint(show, ld) == "SecondaryLimbDark(PolynomialLimbDark, PolynomialLimbDark, 1.0)"

    @test sprint(show, "text/plain", ld) == """
    SecondaryLimbDark
     primary: PolynomialLimbDark([-1.0, 0.4, 0.26])
     secondary: PolynomialLimbDark([-1.0, 0.4, 0.26])
     ratio: 1.0"""
end

@testset "IntegratedLimbDark" begin
    ld = IntegratedLimbDark(u)

    @test sprint(show, ld) == "IntegratedLimbDark(PolynomialLimbDark, 21)"

    @test sprint(show, "text/plain", ld) == """
    IntegratedLimbDark
     driver: PolynomialLimbDark([-1.0, 0.4, 0.26])
     N: 21"""
end

@testset "SimpleOrbit" begin
    orbit = SimpleOrbit(duration=1, period=3)

    @test sprint(show, orbit) == "SimpleOrbit(P=3, T=1, t0=0, b=0)"

    @test sprint(show, "text/plain", orbit) == """
    SimpleOrbit
     period: 3
     duration: 1
     t0: 0
     b: 0"""
end

@testset "KeplerianOrbit" begin
    orbit = KeplerianOrbit(
        ρ_star = 1.0,
        R_star = cbrt(3.0 / (4.0 * π)),
        period = 2.0,
        ecc = 0.0,
        t_0 = 0.0,
        incl = π / 2.0,
        Ω = 0.0,
        ω = 0.0,
    )

    orbit_unit = KeplerianOrbit(
        ρ_star = 1.0u"g/cm^3",
        R_star = cbrt(3.0 / (4.0 * π))u"Rsun",
        period = 2.0u"d",
        ecc = 0.0,
        t_0 = 0.0u"d",
        incl = π / 2.0,
        Ω = 0.0,
        ω = 0.0,
    )

    @test repr("text/plain", orbit) === """
    Keplerian Orbit
     P:      2.0000 d
     t₀:     0.0000 d
     tₚ:     -0.5000 d
     t_ref:  -0.5000 d
     τ:      nothing
     a:      6.6802 R⊙
     aₚ:     -6.6802 R⊙
     aₛ:     0.0000 R⊙
     Rₚ:     0.0000 R⊙
     Rₛ:     0.6204 R⊙
     ρₚ:     NaN M⊙/R⊙³
     ρₛ:     1.0000 M⊙/R⊙³
     r:      0.0000
     aRₛ:    10.7685
     b:      0.0000
     ecc:    0.0000
     cos(i): 0.0000
     sin(i): 1.0000
     cos(ω): 1.0000
     sin(ω): 0.0000
     cos(Ω): 1.0000
     sin(Ω): 0.0000
     i:      1.5708 rad
     ω:      0.0000 rad
     Ω:      0.0000 rad
     Mₚ:     0.0000 M⊙
     Mₛ:     1.0000 M⊙"""

    @test repr("text/plain", orbit_unit) === """
    Keplerian Orbit
     P:      2.0000 d
     t₀:     0.0000 d
     tₚ:     -0.5000 d
     t_ref:  -0.5000 d
     τ:      nothing
     a:      3.6958 R⊙
     aₚ:     -3.6958 R⊙
     aₛ:     0.0000 R⊙
     Rₚ:     0.0000 R⊙
     Rₛ:     0.6204 R⊙
     ρₚ:     NaN M⊙ R⊙^-3
     ρₛ:     0.1693 M⊙ R⊙^-3
     r:      0.0000
     aRₛ:    5.9576
     b:      0.0000
     ecc:    0.0000
     cos(i): 0.0000
     sin(i): 1.0000
     cos(ω): 1.0000
     sin(ω): 0.0000
     cos(Ω): 1.0000
     sin(Ω): 0.0000
     i:      1.5708 rad
     ω:      0.0000 rad
     Ω:      0.0000 rad
     Mₚ:     0.0000 M⊙
     Mₛ:     0.1693 M⊙"""
end


text:
import RowEchelon
import LinearAlgebra

#=
Homework 7
Section 3.2
Problems 10
=#

A = [1 3 -1 0 -2 ;
     0 2 -4 -1 -6 ;
     -2 -6 2 3 9 ;
     3 7 -3 8 -7 ;
     3 5 5 2 7]

detA = LinearAlgebra.det(A)

println("Original Matrix:")
println(A)
println("Determinate:")
println(detA)

ext: jl


text:
export DCTOp

mutable struct DCTOp{T} <: AbstractLinearOperator{T}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod! :: Function
  tprod! :: Nothing
  ctprod! :: Function
  nprod :: Int
  ntprod :: Int
  nctprod :: Int
  args5 :: Bool
  use_prod5! :: Bool
  allocated5 :: Bool
  plan
  dcttype::Int
end

"""
  DCTOp(T::Type, shape::Tuple, dcttype=2)

returns a `DCTOp <: AbstractLinearOperator` which performs a DCT on a given input array.

# Arguments:
* `T::Type`       - type of the array to transform
* `shape::Tuple`  - size of the array to transform
* `dcttype`       - type of DCT (currently `2` and `4` are supported)
"""
function DCTOp(T::Type, shape::Tuple, dcttype=2)

  function dct_multiply(res, plan, x, α, β::T2, shape, T::Type, factor) where {T2}
    if β == zero(T2)
      res .= (factor*α) .* vec(plan*reshape(x,shape))
    else
      res .= (factor*α) .* vec(plan*reshape(x,shape)) .+ β .* q
    end
  end

  if dcttype == 2
    plan = plan_dct(zeros(T,shape))
    iplan = plan_idct(zeros(T,shape))

    prod! = (res, x, α, β)  -> dct_multiply(res, plan, x, α, β, shape, T, one(T))
    tprod! = (res, x, α, β)  -> dct_multiply(res, iplan, x, α, β, shape, T, one(T))

  elseif dcttype == 4
    factor = T(sqrt(1.0/(prod(shape)* 2^length(shape)) ))
    plan = FFTW.plan_r2r(zeros(T,shape),FFTW.REDFT11)
    prod! = (res, x, α, β) -> dct_multiply(res, plan, x, α, β, shape, T, factor)
    tprod! = (res, x, α, β) -> dct_multiply(res, plan, x, α, β, shape, T, factor)
  else
    error("DCT type $(dcttype) not supported")
  end

  return DCTOp{T}(prod(shape), prod(shape), false, false,
                      prod!, nothing, tprod!,
                      0, 0, 0, true, true, true,
                      plan, dcttype)
end

function Base.copy(S::DCTOp)
  return DCTOp(eltype(S), size(S.plan), S.dcttype)
end


text:
module JungleHelperAttachTriggerController

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/AttachTriggerController" AttachTriggerController(x::Integer, y::Integer, entityFilter::String="", triggerFilter::String="")

const placements = Ahorn.PlacementDict(
    "Attach Trigger Controller (Jungle Helper)" => Ahorn.EntityPlacement(
        AttachTriggerController,
        "point",
        Dict{String, Any}(),
        function(trigger)
            trigger.data["nodes"] = [
                (Int(trigger.data["x"]) + 16, Int(trigger.data["y"]))
            ]
        end
    )
)

triggerColor = (0, 127, 14, 1) ./ (255, 255, 255, 1)
triggerEdgeColor = (0, 158, 15, 1) ./ (255, 255, 255, 1)

Ahorn.nodeLimits(entity::AttachTriggerController) = 1, 1

function Ahorn.selection(entity::AttachTriggerController)
    triggerX, triggerY = Int(entity.data["x"]), Int(entity.data["y"])
    entityX, entityY = Int.(entity.data["nodes"][1])

    return [Ahorn.Rectangle(triggerX, triggerY, 8, 8), Ahorn.Rectangle(entityX, entityY, 8, 8)]
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::AttachTriggerController, room::Maple.Room)
    Ahorn.Cairo.save(ctx)

    x, y = Int.(entity.data["nodes"][1])

    Ahorn.rectangle(ctx, x, y, 8, 8)
    Ahorn.clip(ctx)

    Ahorn.drawRectangle(ctx, x, y, 8, 8, triggerColor, triggerEdgeColor)
    Ahorn.drawCenteredText(ctx, "E", x, y, 8, 8)

    Ahorn.restore(ctx)

    xT, yT = Ahorn.position(entity)
    Ahorn.drawArrow(ctx, xT + 4, yT + 4, x + 4, y + 4, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::AttachTriggerController, room::Maple.Room)
    Ahorn.Cairo.save(ctx)

    x, y = Ahorn.position(entity)

    Ahorn.rectangle(ctx, x, y, 8, 8)
    Ahorn.clip(ctx)

    Ahorn.drawRectangle(ctx, x, y, 8, 8, triggerColor, triggerEdgeColor)
    Ahorn.drawCenteredText(ctx, "T", x, y, 8, 8)

    Ahorn.restore(ctx)
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using Documenter, Surrogates
makedocs(
    sitename="Surrogates.jl",
    pages = [
    "index.md"
    "Tutorials" => [
        "Basics" => "tutorials.md",
        "Radials" => "radials.md",
        "Kriging" => "kriging.md",
        "Gaussian Process" =>"abstractgps.md",
        "Lobachevsky" => "lobachevsky.md",
        "Linear" => "LinearSurrogate.md",
        "InverseDistance" => "InverseDistance.md",
        "RandomForest" => "randomforest.md",
        "SecondOrderPolynomial" => "secondorderpoly.md",
        "NeuralSurrogate" => "neural.md",
        "Wendland" => "wendland.md",
        "Polynomial Chaos" => "polychaos.md",
        "Variable Fidelity" => "variablefidelity.md",
        "Gradient Enhanced Kriging" => "gek.md"
       ]
    "User guide" => [
        "Samples" => "samples.md",
        "Surrogates" => "surrogate.md",
        "Optimization" => "optimizations.md"
        ]
    "Benchmarks" => [
        "Sphere function" => "sphere_function.md",
        "Lp norm" => "lp.md",
        "Rosenbrock" => "rosenbrock.md",
        "Tensor product" => "tensor_prod.md",
        "Cantilever beam" => "cantilever.md",
        "Water Flow function" => "water_flow.md",
        "Welded beam function" => "welded_beam.md",
        "Branin function" => "BraninFunction.md",
        "Ackley function" => "ackley.md",
        "Gramacy & Lee Function" => "gramacylee.md",
        "Salustowicz Benchmark" => "Salustowicz.md",
        "Multi objective optimization" => "multi_objective_opt.md"
        ]
    "Contributing" => "contributing.md"
    ]
)

deploydocs(
   repo = "github.com/SciML/Surrogates.jl.git",
)


text:
# Model
include_model("biped")
model = free_time_model(model)

# Visualize
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

θ = pi / 10.0
# q1 = initial_configuration(model, θ)
q1, qT = loop_configurations(model, θ)

visualize!(vis, model, [q1])

# Horizon
T = 21

# Time step
tf = 1.0
h = tf / (T-1)

# Bounds

# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= model.uU
_uu[end] = h
_ul = zeros(model.m)
_ul[model.idx_u] .= model.uL
_ul[end] = 0.01 * h
ul, uu = control_bounds(model, T, _ul, _uu)

xl, xu = state_bounds(model, T,
    x1 = [Inf * ones(model.nq); q1],
    xT = [Inf * ones(model.nq); qT])

# Objective
include_objective(["velocity", "nonlinear_stage"])
q_ref = linear_interpolation(q1, qT, T)
x0 = configuration_to_state(q_ref)

obj_penalty = PenaltyObjective(1.0e5, model.m - 1)

Qq = Diagonal(ones(model.nq))
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
QT = cat(0.5 * Qq, 0.5 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0e-1 * ones(model.nu)..., zeros(model.m - model.nu)...])

obj_tracking = quadratic_time_tracking_objective(
    [t < T ? Q : QT for t = 1:T],
    # [Diagonal(zeros(model.n)) for t = 1:T],
    [R for t = 1:T-1],
    [zeros(model.n) for t = 1:T],
    [zeros(model.m) for t = 1:T],
    1.0)
obj_velocity = velocity_objective(
    [Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
    model.nq)

function l_stage_fh(x, u, t)
    if t > 1
        return 1.0 * (kinematics_2(model,
            view(x, 8:14), body = :leg_2, mode = :ee)[2] - 0.025)^2.0
    else
        return 0.0
    end
end
l_terminal_fh(x) = 0.0
obj_fh = nonlinear_stage_objective(l_stage_fh, l_terminal_fh)
obj = MultiObjective([obj_tracking, obj_penalty, obj_velocity, obj_fh])

# Constraints
include_constraints(["contact", "pinned_foot", "free_time", "loop"])
con_contact = contact_constraints(model, T)
con_pinned = pinned_foot_constraint(model, q1, T)
con_free_time = free_time_constraints(T)
# con_loop = loop_constraints(model, collect([(2:7)...,(9:14)...]), 1, T)
con = multiple_constraints([con_contact, con_free_time, con_pinned])#, con_loop])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
q_ref = linear_interpolation(q1, qT, T)
x0 = configuration_to_state(q_ref) # linear interpolation on state
u0 = [1.0e-5 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
include_snopt()

@time z̄ , info = solve(prob, copy(z0),
    nlp = :SNOPT7,
    tol = 1.0e-3, c_tol = 1.0e-3,
    time_limit = 60 * 2)

check_slack(z̄, prob)

x̄, ū = unpack(z̄, prob)
tf, t, hc = get_time(ū)

# Visualize
visualize!(vis, model, state_to_configuration(x̄), Δt = hc[1])


text:
struct IterativeBatchQs <: MultiObjectiveBatchQs
    model::SVDD.OCClassifier
    inf_measure::SequentialPQs
    rep_measure::Function
    div_measure::Function
    normalization::Function
    k::Int
    λ_inf::Float64
    λ_rep::Float64
    λ_div::Float64

    function IterativeBatchQs(model::SVDD.OCClassifier, informativeness::SequentialPQs; representativeness::Symbol, diversity::Symbol,
        k::Int, λ_inf::T1 where T1<:Real=1, λ_rep::T2 where T2 <: Real=1, λ_div::T3 where T3 <: Real=1)::IterativeBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))

        representativeness_measure = get_rep_measure(representativeness)
        diversity_measure = get_incremental_div_measure(diversity)
        return new(model, informativeness, representativeness_measure, diversity_measure, min_max_normalization, k, λ_inf, λ_rep, λ_div)
    end
end

"""
Select best batch with weighted sum of requirements.
Iterative selection: diversity is only computed to previously selected observations.
"""
function select_batch(qs::IterativeBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    #informativeness
    inf_scores_normalized = qs.normalization(qs_score(qs.inf_measure, x, labels)[candidate_indices])

    # representativeness
    rep_scores_normalized = qs.normalization(qs.rep_measure(qs.model, x, labels, candidate_indices))

    batch_samples = Int[]
    div_scores = Float64[]
    for iteration in 1:qs.k
        combined_scores = Vector{Float64}(undef, num_observations)
        if length(batch_samples) == 0
            # first sample for new batch cannot compute diversity to existing samples
            combined_scores = qs.λ_inf * inf_scores_normalized + qs.λ_rep * rep_scores_normalized
        else
            div_scores = qs.div_measure(qs.model, x, batch_samples[end], candidate_indices, div_scores)
            # unormalized div_scores are needed for incremental diversity computation
            # normalize div_scores to combine scores
            combined_scores = qs.λ_inf * inf_scores_normalized + qs.λ_rep * rep_scores_normalized + qs.λ_div * qs.normalization(div_scores)
            batch_sample_indices = [ind for (ind, val) in enumerate(candidate_indices) if val in batch_samples]
            combined_scores[batch_sample_indices] .= -Inf
        end
        best_sample_index = candidate_indices[argmax(combined_scores)]

        push!(batch_samples, best_sample_index)
    end
    return batch_samples
end


text:
#use transition matrix defined in each model file, but general simulation method is here
#specific for fitting figs. 4a and 4b from Ramentol et al., 2020

using DifferentialEquations, OrdinaryDiffEq
using Random, LinearAlgebra, Statistics
using GCMAES, BlackBoxOptim, ForwardDiff
using Glob, CSV, DataFrames, Dates
using Plots

"""
Notes on how performance can be improved:
* Modify the objective function so that `db_pars` is called only once per evoluation of the objective function.
    i.e. `tm(...)` should not have to call `db_pars`
* Remove bounds-checking from the objective function. Originally, bounds-checking was added for
    straightforward use with algorithms that didn't support constraints.
* Pre-allocate a steady-state vector that is modified in-place by `get_ss(...)`
* If using GCMAES for optimization, the initial step-size (coordinate-wise standard deviation) can be tuned.
    Similarly, parameters for differential evolution can be tuned.
    For more details, see the documentation for the respective packages (GCMAES.jl and BlackBoxOptim.jl)
"""

# load model
model_name = "...."
include("[path to model files]/$model_name")
println("Using model $model_name")

# path to data .csv files
data_path = "..."

### Conditions
DB_CONDITION = "c"
# BB for BlackBoxOptim or GCMAES
OPT_OR_NOT = "GCMAES"
FA_OR_WT = "mut"
OPT_MAX_ITER = 800

# holding voltage for activation `va` and deactivation `vd`
if FA_OR_WT == "mut"
    vhold_a = -80
    vhold_d = -80
else
    vhold_a = -10
    vhold_d = 40
end

SAVE_FIGS = false
if SAVE_FIGS == true
    SAVE_NAME = "wt-d1"
end

#open current time course data as dataframe
function reader(f)
    filename = string(basename(f))[1:end-4]
    df = CSV.read(f, DataFrame)

    #rename voltage(mV) in gv/fv df to Voltage
    if Symbol("Voltage(mV)") in names(df)
        rename!(df, Symbol("Voltage(mV)") => :Voltage)
    end

    #HOW TO INDEX DF WITH CONDITIONS:
    # df[df.Voltage .== -100, 2][1]
    # Note [1]-index at the end; this is needed b/c indexing a df returns an array, even if of size 1.

    return df
end

#################################################
"""LOAD DATA FILES"""
#################################################

if FA_OR_WT == "mut"
    ss_path = string(data_path, "mut_ss.csv")
    ia_path = string(data_path, "mut_it_act.csv")
    id_path = string(data_path, "mut_it_de.csv")
    fa_path = string(data_path, "mut_ft_act.csv")
    fd_path = string(data_path, "mut_ft_de.csv")
else
    ss_path = string(data_path, "wt_ss.csv")
    ia_path = string(data_path, "wt_it_act.csv")
    id_path = string(data_path, "wt_it_de.csv")
    fa_path = string(data_path, "wt_ft_act.csv")
    fd_path = string(data_path, "wt_ft_de.csv")
end

files = [glob(x)[1] for x in [ss_path, ia_path, id_path, fa_path, fd_path]]
# println(files)
# asdf

data_files = [reader(x) for x in files]
data_files[2:end] = [x[:, 2:end] for x in data_files[2:end]]
# println(data_files[1])

#Select detailed balance parameter
if DB_CONDITION == true && "four" in model_name
    DB_CONDITION = 'c'
elseif DB_CONDITION == true && model_name == "six_v4.jl"
    DB_CONDITION = true
end

println("Detailed balance using $DB_CONDITION")

# extract voltage, times from data
function extract_VoltageAndTimes(data::Array{DataFrame, 1};
    convert_to_ms=false, convert_to_s=false)
    """
    Assume there are 2 datasets
    data = Array of dataframes
    convert_to_ms = whether to convert time in data from s to ms
    """
    # number of dataframes, i.e. datasets
    Nd = length(data)

    voltages = Array{Array{Float64,1}}(undef, Nd)
    times = Array{Array{Array{Float64,1}}}(undef, Nd)
    trace = Array{Array{Array{Float64,1}}}(undef, Nd)

    for i = 1:Nd
        colnames = names(data[i])

        # times and voltage columns are present as: "-150_t, -150_i, etc.", so we remove the last two characters
        # since the columns are paired, we take every second column for time or voltage
        v = [parse(Float64, String(n)[1:end-2]) for n in colnames[1:2:end]]
        voltages[i] = v

        # number of traces is half the number of columns, since we have time columns for each trace
        Ntraces = Int(length(colnames)/2)
        tᵢ = Array{Array{Float64,1}}(undef, Ntraces)
        yᵢ = Array{Array{Float64,1}}(undef, Ntraces)

        for j = 1:2:length(colnames)

            if convert_to_ms == true
                data[i][:,j] .*= 1000
            elseif convert_to_s == true
                data[i][:,j] .*= 1e-3
            end

            tⱼ = collect( skipmissing( data[i][:,j] ) )
            # tⱼ = (tⱼ .- tⱼ[1]) .* 1000
            tⱼ = sort(abs.(tⱼ .- tⱼ[1]))
            # println(typeof(tⱼ))

            u = Int((j+1)/2)
            tᵢ[u] = tⱼ
            yᵢ[u] = collect(skipmissing(data[i][:,j+1]))
            if length(yᵢ[u]) < 1
                println(data[i][:,j:j+1])
            end
        end

        times[i] = tᵢ
        trace[i] = yᵢ
    end

    # indices of voltages in 1st dataset that are also in the second
    shared_idx = findall(in(voltages[1]), voltages[2])

    return Dict("voltages" => voltages, "times" => times, "shared_idx" => shared_idx, "traces" => trace)
end

# voltage, times, and shared indices for current-time
extracted_i = extract_VoltageAndTimes(data_files[2:3])
extracted_f = extract_VoltageAndTimes(data_files[4:5])

#################################################
"""SIMULATION METHODS"""
#################################################
function moving_average(vs::Vector,n::Int)
    return [sum(@view vs[i:(i+n-1)])/n for i in 1:(length(vs)-(n-1))]
end

function get_ss(p::AbstractVector{}, vhold, db)::Vector{}
    Q = tm(vhold, p, db)
    return Q[:,2:end] \ (-1 .* Q[:,1])
end

function solve_odes(inits, Q, times::Array{Float64,1})
    prob = ODEProblem(rhs!, eltype(Q).(inits), (times[1], times[end]), Q)
    sol = solve(prob, TRBDF2(), saveat = times, reltol=1e-8, abstol=1e-8, maxiters=1e6)
    return sol
end

function one_sim(pars::AbstractVector{T}, data,
                u0::AbstractVector{}, Qhold, db;
                pmax_min=0.7, states="open") where T

    va, vd = data["voltages"]
    ta, td = data["times"]
    shared = data["shared_idx"]

    y_act = Array{Array{T, 1}, 1}(undef, length(va))
    y_de = Array{Array{T, 1}, 1}(undef, length(vd))

    j = 0
    for i in 1:length(va)
        Qa = tm(va[i], pars, db)
        sol = solve_odes(u0, Qa, ta[i])

        if states == "open"
            y_act[i] = vec(sum(sol[oi:end,:], dims=1))
        else
            y_act[i] = vec(sum(sol[ai,:], dims=1))
        end

        if va[i] in vd
            j += 1
            h = shared[j]

            # deactivation at `vhold_d`
            Qd = tm(vhold_d, pars, db)
            sol = solve_odes(sol[end], Qd, td[h])

            if states == "open"
                y_de[h] = vec(sum(sol[oi:end,:], dims=1))
            else
                y_de[h] = vec(sum(sol[ai,:], dims=1))
            end
        end
    end

    return y_act, y_de
end

function one_sim_LossOnly(pars::AbstractVector{}, data, RawData,
                        u0::AbstractVector{}, Qhold, db;
                        pmax_min=0.7, states="open", method="sse")

    va, vd = data["voltages"]
    ta, td = data["times"]
    shared = data["shared_idx"]
    traces = data["traces"]

    E = 0.
    j = 0

    for i in 1:length(va)
        Qa = tm(va[i], pars, db)
        sol = solve_odes(u0, Qa, ta[i])

        if states == "open"
            y_act = vec(sum(sol[oi:end,:], dims=1))
            x = traces[1][i]
        else
            y_act = vec(sum(sol[ai,:], dims=1))
            x = traces[1][i]
        end

        E += sum( (x .- y_act) .^ 2 )/length(ta[i])

        if va[i] in vd
            j += 1
            h = shared[j]

            sol = solve_odes(sol[end], tm(vhold_d, pars, db), td[h])

            if states == "open"
                y_de = vec(sum(sol[oi:end,:], dims=1))
                x = traces[2][h]
            else
                y_de = vec(sum(sol[ai,:], dims=1))
                x = traces[2][h]
            end

            E += sum( (x.- y_de) .^ 2 )/length(td[h])
        end
    end
    return E
end

function sim(pars::AbstractVector{};
            data_i = extracted_i,
            RawData_i = data_files[2:3],
            data_f = extracted_f,
            RawData_f = data_files[4:5],
            db = DB_CONDITION,
            loss_only = false,
            error_method = "sse"
            )
    """
    data_i and data_f are Dicts with keys `voltages, times, and shared_idx`
        voltages = Array{Array{Float64,1}}
        times = Array{Array{Array{Float64,1}}}
        shared_idx = list of indices of voltages in the 2nd dataset that are also in the 1st dataset
        `error_method` = 'sse' for sum of squared errors, 'LLH' for log likelihood
    """
    #initial state vector (steady state at vhold)
    Qhold = tm(vhold_a, pars, db)
    inits = get_ss(pars, vhold_a, db)

    if loss_only == true
        if error_method == "sse"
            loss = one_sim_LossOnly(pars, data_i, RawData_i, inits, Qhold, db;
                                    pmax_min=0.9, states="open", method="sse")
            loss += one_sim_LossOnly(pars, data_f, RawData_f, inits, Qhold, db;
                                    pmax_min=0.85, states="active", method="sse")
            return loss
        elseif error_method == "LLH"
            loss = one_sim_LossOnly(pars, data_i, RawData_i, inits, Qhold, db;
                                    pmax_min=0.9, states="open", method="LLH")
            loss += one_sim_LossOnly(pars, data_f, RawData_f, inits, Qhold, db;
                                    pmax_min=0.9, states="active", method="LLH")
            return loss
        else
            error("`error_method` can only be 'sse' or 'LLH'.")
        end
    else
        sim_ia, sim_id = one_sim(pars, data_i, inits, Qhold, db; states="open")
        sim_fa, sim_fd = one_sim(pars, data_f, inits, Qhold, db; states="active")
        return sim_ia, sim_id, sim_fa, sim_fd
    end
end

#get a gv curve
function sim_ss(pars::AbstractVector{T}; data_ss = data_files[1], db = DB_CONDITION, loss_only=false) where T

    y = zeros(T, (size(data_ss)[1], 2))
    for i in 1:size(data_ss)[1]
        # println(data_ss)
        v = data_ss[i, 1]

        states = get_ss(pars, v, db)

        y[i,1] = sum(states[oi:end])
        y[i,2] = sum([states[a] for a in ai])
    end

    if loss_only == true
        E = 0.
        for i in 1:2
            if normalize
                y[:,i] .*= 1/maximum(y[:,i])
            end
            E += sum((data_ss[:,i+1] .- y[:,i]) .^ 2)/length(y[:,i])
        end
        return E
    else
        if normalize
            for i in 1:2
                y[:,i] .*= 1/maximum(y[:,i])
            end
        end
        return [data_ss[:,1], y]
    end
end

#################################################
# LOSS AND OBJECTIVE FUNCTIONS
#################################################
function check_ranges()
    # check if popt is in ranges; if not, return closest bound
    in_range = copy(popt)

    for (i, b) in enumerate(ln_par_bds)
        if b[1] < popt[i] < b[2]
            continue
        else
            println("$i th parameter is out of bounds, \n   ",
                    popt[i], " outside  $b")
            if popt[i] > b[2]
                in_range[i] = b[2]
            else
                in_range[i] = b[1]
            end
        end
    end

    if in_range == popt
        println("\n All parameters in-range.")
        println("\n Optimized cost = ", objective(popt))
        println("\n Exp params:    ")
        println(exp.(popt))
        return popt
    else
        println("\n Optimized cost = ", objective(in_range))
        println("\n Exp params:    ")
        println(exp.(in_range))
        return in_range
    end
end

function objective(pars)
    # check for negative values in db_pars
    if any(x -> x < 0, db_pars(pars, DB_CONDITION))
        return 1e6
    else
        # check whether parameters are in bounds
        for (i, b) in enumerate(ln_par_bds)
            if b[1] <= pars[i] <= b[2]
                continue
            else
                # println("   outside bound ", i)
                return 1e6
            end
        end

        # check that maximal steady state open probability is > 0.9
        # if sum(get_ss(pars, 200, DB_CONDITION)[oi:end]) < 0.5
            # println("   maximal GV below threshold")
            # return 1e6
        # end
    end

    return sim(pars; loss_only=true)
end

# gradient of objective function
∇objective(pars) = ForwardDiff.gradient(objective, pars)

#################################################
# OPTIMIZATION
#################################################
# remove parameters and bounds defined by detailed balance
p0, ln_par_bds = prepare_BoundsAndParams(p0, ln_par_bds, DB_CONDITION)

# force parameters to closest bound if outside bounds
for (i, b) in enumerate(ln_par_bds)
    if b[1] <= p0[i] <= b[2]
        continue
    elseif p0[i] < b[1]
        p0[i] = b[1]
    else
        p0[i] = b[2]
    end
end

println("   Initial objective: ", objective(p0))
# asdf

if OPT_OR_NOT == "GCMAES"
    popt, fmin, status = GCMAES.minimize((objective, ∇objective), p0, 0.4,
                                        [b[1] for b in ln_par_bds],
                                        [b[2] for b in ln_par_bds],
                                        maxiter = OPT_MAX_ITER
    )
    println("----------// RESULTS OF GCMA_ES //-----------",
            "\n     Status  ", status,
            "\n     Fmin = ", fmin,
            "\n     Bounds = ", ln_par_bds,
            "\n     Full parameters... \n Real-space  ",
            db_pars(popt, DB_CONDITION),
            "\n Full parameters... \n Log-space  ",
            log.(db_pars(popt, DB_CONDITION)))
elseif OPT_OR_NOT == "BB"
    optctrl = bbsetup(objective,
                    SearchRange = ln_par_bds,
                    Method = :adaptive_de_rand_1_bin_radiuslimited,
                    MaxStepsWithoutProgress = 2000,
                    MinDeltaFitnessTolerance = 1e-12, #regular = 1e-12
                    PopulationSize=250,
                    TraceMode=:verbose,
                    TraceInterval=30.0,
                    MaxSteps= OPT_MAX_ITER,
                    MaxTime = 7200 #2h
                    # CallbackFunction = callback,
                    # CallbackInterval = 60.0
        )

    res = bboptimize(optctrl)
    popt = best_candidate(res)

    println("----------// RESULTS OF BBOPTIMIZE //-----------",
        "\n     Bounds... \n ", ln_par_bds,
        "\n     Full parameters... \n Real-space  ",
        db_pars(popt, DB_CONDITION),
        "Full parameters... \n Log-space  ",
        db_pars(popt, DB_CONDITION; fmt="log"),
        )
else
    popt = p0
end

#################################################
# OUTPUT
#################################################
println("   Objective value: ", objective(popt))

# apply detailed balance to optimized parameters (or p0, if no optimization was done)
# output in log format
popt_full = db_pars(popt, DB_CONDITION; fmt="log")

#################################################
# PLOTTING
#################################################
function plotter()
    """
    Uses `popt`, which is assumed to be log-transformed.
    Detail balance will be applied during simulation, so `popt` should not be full length if `DB_CONDITION` is not `'x'`.
    If optimization was not performed, `popt = p0`.
    `p0` should be defined in the corresponding model file.
    """
    #simulate with optimal parameters
    sim_ia, sim_id, sim_fa, sim_fd = sim(popt; loss_only = false)
    volts, y_ss = sim_ss(popt; loss_only = false)

    # steady state
    plt_g = plot(
            xlabel="Voltage (mV)",
            ylabel="Normalized Conductance",
            legend=:outertopright
        );

    plt_f = plot(
        xlabel="Voltage (mV)",
        ylabel="Normalized dFluorescence",
        legend=:outertopright
    );

    labels = ["GV", "FV"]
    plts = [plt_g, plt_f]

    # plot steady state data
    m = Matrix(data_files[1])
    for j in 2:size(m)[2]
        plot!(plts[j-1], m[:, 1], m[:, j],
                markershape = :circle, markercolor = :blue,
                linecolor = :blue, linealpha = 0.3,
                linewidth = 4, label = "Data/" * labels[j-1]
        )

        plot!(plts[j-1], m[:, 1], y_ss[:, j-1],
                linestyle = :dash, linealpha = 0.5,
                linecolor = :red, linewidth = 4,
                markershape = :circle, markercolor = :red,
                label = "Sim/" * labels[j-1]
        )

        display(plts[j-1])
    end

    function plot_TCourse(title)
        if title == "current"
            ylab = "Normalized Current"
            ex = extracted_i
            y_a, y_d = sim_ia, sim_id
            println([length(x) for x in y_d])
        else
            ylab = "Normalized Fluor"
            ex = extracted_f
            y_a, y_d = sim_fa, sim_fd
        end

        plt_a = plot(
            xlabel="Time (ms)",
            ylabel=ylab,
            legend=:outertopright
        )
        plt_d = plot(
            xlabel="Time (ms)",
            ylabel=ylab,
            legend=:outertopright
        )

        va, vd = ex["voltages"]
        ta, td = ex["times"]
        traces = ex["traces"]
        # shared = ex["shared"]

        for (i, v) in enumerate(va)
            t = ta[i]
            y = traces[1][i]

            plot!(plt_a, t, y,
                seriestype = :scatter,
                markeralpha = 0.35,
                label = string(v)
            )

            # println(length(t), "    ", length(y_a[i]))
            if length(y_a[i]) > 1
                plot!(plt_a, t[1:length(y_a[i])], y_a[i],
                        linestyle = :dash,
                        linewidth = 4,
                        label = ""
                )
            end
        end

        for (i, v) in enumerate(vd)
            # v = vd[i]
            t = td[i]
            y = traces[2][i]

            try
                plot!(plt_d, t, y,
                    seriestype = :scatter,
                    markeralpha = 0.35,
                    label = string(v)
                )
            catch e
                println(e)
            end

            # println(length(t), "    ", length(y_d[i]))
            if length(y_d[i]) > 1
                # plot!(plt_d, t[1:length(y_d[i])], y_d[i],
                plot!(plt_d, t, y_d[i],
                        linestyle = :dash,
                        linewidth = 4,
                        label = ""
                )
            end
        end

        display(plt_a)
        display(plt_d)
    end

    try
        plot_TCourse("current")
    catch e
        println(e)
    end
    try
        plot_TCourse("fluorescence")
    catch e
        println(e)
    end
end

plotter()

max_stars_repo_name: haganenoneko/FunnyChannelGatingReview

max_issues_repo_name: haganenoneko/FunnyChannelGatingReview

max_forks_repo_name: haganenoneko/FunnyChannelGatingReview


text:
### Instructions for adding a new version of the Go toolchain
#
# * find the latest stable releases at https://go.dev/dl/, update the `version`
#   variable and the SHA256 hash of the release tarball in the sources, the
#   expected checksum is provided in the download page.
# * To deploy the shard and automatically update your BinaryBuilderBase's
#   `Artifacts.toml`, use the `--deploy` flag to the `build_tarballs.jl` script.
#   You can build & deploy by running:
#
#      julia build_tarballs.jl --debug --verbose --deploy
#

using BinaryBuilder

include("../common.jl")

name = "Go"
version = v"1.17.7"

sources = [
    ArchiveSource(
        "https://go.dev/dl/go$(version).linux-amd64.tar.gz",
        "02b111284bedbfa35a7e5b74a06082d18632eff824fd144312f6063943d49259",
    )
]

# Bash recipe for building across all platforms
script = raw"""
mv ${WORKSPACE}/srcdir/go ${prefix}/
"""

# We only build for host platform: x86_64-linux-musl
platforms = [
    host_platform,
]

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

# The products that we will ensure are always built
products = [
    ExecutableProduct("go", :go, "go/bin"),
    ExecutableProduct("gofmt", :gofmt, "go/bin"),
]

# Build the tarballs, and possibly a `build.jl` as well.
ndARGS, deploy_target = find_deploy_arg(ARGS)
build_info = build_tarballs(ndARGS, name, version, sources, script, platforms, products, dependencies; skip_audit=true)

if deploy_target !== nothing
    upload_and_insert_shards(deploy_target, name, version, build_info)
end

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
const NodeID = Int64

"""
A SequenceGraphNode represents nodes or vertices in a sequence graph.

Conceptually, a node in a sequence graph is a node which posesses a
biological sequence, such nodes may be connected to other nodes through
'Links'.

SequenceGraphNodes are directional. They have a start (marked '+'),
and they have an end (marked '-').

Furthuremore, nodes are canonical or palindromic sequences, else they are reverted.
"""
struct SequenceGraphNode{S <: BioSequence}
    sequence::S
    active::Bool
end

sequence(sn::SequenceGraphNode) = sn.sequence

isactive(sn::SequenceGraphNode) = sn.active

reverse_complement!(sn::SequenceGraphNode) = reverse_complement!(sequence(sn))

length(sn::SequenceGraphNode) = length(sequence(sn))

function iscanonical(seq::BioSequence{DNAAlphabet{2}})
    i = 1
    j = length(seq)
    @inbounds while i < j
        f = seq[i]
        r = complement(seq[j])
        f < r && return true
        r < f && return false
        i += 1
        j -= 1
    end
    return true
end

iscanonical(sn::SequenceGraphNode) = iscanonical(sequence(sn))


text:
function mixFLDO!(;
    qs   :: AbstractArray{Float64, 1},
    zs   :: AbstractArray{Float64, 1},
    hs   :: AbstractArray{Float64, 1},
    q_ML :: Float64,
    FLDO :: Integer,
    FLDO_ratio_top :: Float64,
    FLDO_ratio_bot :: Float64,
)

    Δq = 0.0

    if FLDO != -1

        Δq = q_ML - qs[FLDO]
        qs[FLDO] =  FLDO_ratio_top * q_ML + FLDO_ratio_bot * qs[FLDO]

    end

    return Δq

end

function unmixFLDOKeepDiff!(;
    qs   :: AbstractArray{Float64, 1},
    zs   :: AbstractArray{Float64, 1},
    hs   :: AbstractArray{Float64, 1},
    h_ML :: Float64,
    FLDO :: Integer,
    Nz   :: Integer,
    Δq   :: Float64,
#    verbose = false
)

    int_layer = 0
    integral = 0.0
    new_q_ML = 0.0

    if FLDO == -1

        for k = 1:Nz
            integral += hs[k] * qs[k]
        end
        new_q_ML = integral / h_ML
        qs[1:Nz] .= new_q_ML

    else

#        verbose && println("FLDO = ", FLDO)
#        verbose && println("Δq: ", Δq)
        for k = 1:FLDO
            integral += hs[k] * qs[k]
#            verbose && println(k, "; hs: ", hs[k], "; qs: ", qs[k])
        end
#        verbose && println("integral = ", integral)
#        verbose && println("h_ML = ", h_ML, "; zs[FLDO+1] = ", zs[FLDO+1])

        new_q_ML = (integral - Δq * (h_ML + zs[FLDO+1])) / ( - zs[FLDO+1] )
        qs[1:FLDO] .= new_q_ML
        qs[FLDO]    = new_q_ML - Δq

    end

    return new_q_ML

end

function remixML!(;
    qs   :: AbstractArray{Float64, 1},
    zs   :: AbstractArray{Float64, 1},
    hs   :: AbstractArray{Float64, 1},
    h_ML :: Float64,
    FLDO :: Integer,
    Nz   :: Integer,
)

    # no need to remix if h_ML is shallower than the first layer
    if FLDO == 1
        return qs[1]
    end

    int_layer = 0
    new_q_ML = 0.0

    if FLDO == -1
        int_layer = Nz
    else
        int_layer = FLDO - 1
        new_q_ML += ( h_ML + zs[FLDO] ) * qs[FLDO]
    end

    for k = 1:int_layer
        new_q_ML += hs[k] * qs[k]
    end
    #=
    if go
        println("hs: ", hs)
        println("h_ML: ", h_ML)
        println("FLDO: ", FLDO)
        println("qs[1:FLDO]: ", qs[1:FLDO])
        println("zs[FLDO]: ", zs[FLDO])
        println("int_layer: ", int_layer)

    end

=#

    new_q_ML /= h_ML

    qs[1:int_layer] .= new_q_ML
    return new_q_ML

#    qs[1:int_layer] .= qs[1]
    return qs[1]

end

#=

function OC_getIntegratedTemperature(
    ocn      :: Ocean,
    i        :: Integer,
    j        :: Integer;
    target_z :: Float64,
)

    return getIntegratedQuantity(
        zs       = ocn.cols.zs[i, j],
        qs       = ocn.cols.Ts[i, j],
        q_ML     = ocn.T_ML[i, j],
        h_ML     = ocn.h_ML[i, j],
        Nz       = ocn.Nz[i, j],
        target_z = target_z,
    )
end

function OC_getIntegratedSalinity(
    ocn      :: Ocean,
    i        :: Integer,
    j        :: Integer;
    target_z :: Float64,
)

    return getIntegratedQuantity(
        zs       = ocn.cols.zs[i, j],
        qs       = ocn.cols.Ss[i, j],
        q_ML     = ocn.S_ML[i, j],
        h_ML     = ocn.h_ML[i, j],
        Nz       = ocn.Nz[i, j],
        target_z = target_z,
    )
end

function OC_getIntegratedBuoyancy(
    ocn      :: Ocean,
    i        :: Integer,
    j        :: Integer;
    target_z :: Float64,
)

    return getIntegratedQuantity(
        zs       = ocn.cols.zs[i, j],
        qs       = ocn.cols.bs[i, j],
        q_ML     = ocn.b_ML[i, j],
        h_ML     = ocn.h_ML[i, j],
        Nz       = ocn.Nz[i, j],
        target_z = target_z,
    )
end

function getIntegratedQuantity(;
    zs       :: AbstractArray{Float64,1},
    qs       :: AbstractArray{Float64,1},
    q_ML     :: Float64,
    h_ML     :: Float64,
    Nz       :: Integer,
    target_z :: Float64,
)

    if target_z < zs[Nz+1]
        throw(ErrorException(
            format("target_z ({:f}) cannot be deeper than the minimum of zs ({:f}).", target_z, zs[Nz])
        ))
    end

    # Integrate mixed layer
    if -target_z < h_ML
        return q_ML * ( - target_z )
    end

    sum_q = 0.0
    sum_q += h_ML * q_ML

    # Test if entire ocean column is mixed layer
    FLDO = getFLDO(zs=zs, h_ML=h_ML, Nz=Nz)
    if FLDO == -1
        return sum_q
    end

    # Integrate FLDO
    if target_z > zs[FLDO+1]
        sum_q += qs[FLDO] * ( (-h_ML) - target_z)
        return sum_q
    end

    sum_q += qs[FLDO] * ( (-h_ML) - zs[FLDO+1])

    # Integrate rest layers
    if FLDO < Nz
        for i = FLDO+1 : Nz
            if target_z < zs[i+1]
                sum_q += qs[i] * (zs[i] - zs[i+1])
            else
                sum_q += qs[i] * (zs[i] - target_z)
                return sum_q
            end
        end
    else
        return sum_q
    end

end
=#

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text: abstract type ElectromagneticField end


text:
function yinyang(n::Int=3)
    radii   = (i * n for i in (1, 3, 6))
    ranges  = collect(collect(-r:r) for r in radii)
    squares = collect(collect((x, y) for x in rnge, y in rnge) for rnge in ranges)
    circles = collect(collect((x, y) for (x,y) in sqrpoints if hypot(x, y) ≤ radius)
                      for (sqrpoints, radius) in zip(squares, radii))
    m = Dict((x, y) => ' ' for (x, y) in squares[end])
    for (x, y) in circles[end] m[(x, y)] = x > 0 ? '·' : '*' end
    for (x, y) in circles[end-1]
        m[(x, y + 3n)] = '*'
		m[(x, y - 3n)] = '·'
    end
    for (x, y) in circles[end-2]
        m[(x, y + 3n)] = '·'
		m[(x, y - 3n)] = '*'
    end
    return join((join(m[(x, y)]  for x in reverse(ranges[end])) for y in ranges[end]), '\n')
end

println(yinyang(4))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ---
# title: 1389. Create Target Array in the Given Order
# id: problem1389
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/create-target-array-in-the-given-order/description/>
# hidden: true
# ---
#
# Given two arrays of integers `nums` and `index`. Your task is to create
# _target_ array under the following rules:
#
#   * Initially _target_ array is empty.
#   * From left to right read nums[i] and index[i], insert at index `index[i]` the value `nums[i]` in  _target_ array.
#   * Repeat the previous step until there are no elements to read in `nums` and `index.`
#
# Return the _target_ array.
#
# It is guaranteed that the insertion operations will be valid.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [0,1,2,3,4], index = [0,1,2,2,1]
#     Output: [0,4,1,3,2]
#     Explanation:
#     nums       index     target
#     0            0        [0]
#     1            1        [0,1]
#     2            2        [0,1,2]
#     3            2        [0,1,3,2]
#     4            1        [0,4,1,3,2]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,4,0], index = [0,1,2,3,0]
#     Output: [0,1,2,3,4]
#     Explanation:
#     nums       index     target
#     1            0        [1]
#     2            1        [1,2]
#     3            2        [1,2,3]
#     4            3        [1,2,3,4]
#     0            0        [0,1,2,3,4]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1], index = [0]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length, index.length <= 100`
#   * `nums.length == index.length`
#   * `0 <= nums[i] <= 100`
#   * `0 <= index[i] <= i`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import Lorenzo: smaller_step, extract_date, check_steps,
                dates_directories, create_urls
@testset "scraper" begin
    @testset "steps" begin
        @test smaller_step("y/m/S") == 'S'
        @test smaller_step("H") == 'H'
        @test_throws  ErrorException smaller_step("N")
    end

    @testset "extract" begin
        @test extract_date("this is %Y%d") == "Yd"
        @test extract_date("%m%n/something/%d") == "mnd"
        @test extract_date("There's no pattern") == ""
    end

    @testset "Check" begin
        @test check_steps("md")
        @test check_steps("mnd") == false
    end

    @testset "extractDir" begin
        @test dates_directories("http://my_path/%Y/%m/%d/example%H%M.txt") == "Ymd"
        @test dates_directories("http://my_path/%Y/%m/example%H%M.txt") == "Ym"

    end

    @testset "urls_directories" begin
        url = "http://my_path/%Y/%m/%d/example%H%m.txt"
        dt0 = DateTime(2010,1,1,1,5)
        dt1 = DateTime(2010,1,1,10,5) # 10 hours range
        @test collect(create_urls(url, dt0, dt1)) == ["http://my_path/2010/01/01/"]
        url = "http://my_path/%Y/%m/%d/%H/example%H%m.txt"
        @test collect(create_urls(url, dt0, dt1)) == ["http://my_path/2010/01/01/$(dec(x,2))/" for x in 1:10]
        url = "http://my_path/%Y/%j/example%H%m.txt"
        dt1 = DateTime(2010,2,3,10,5) # 10 hours range
        @test collect(create_urls(url, dt0, dt1)) == ["http://my_path/2010/$(dec(x,3))/" for x in 1:34]

    end
end


text:
# Context version
struct PointwiseLikelihoodContext{A,Ctx} <: AbstractContext
    loglikelihoods::A
    context::Ctx
end

function PointwiseLikelihoodContext(
    likelihoods=Dict{VarName,Vector{Float64}}(),
    context::AbstractContext=LikelihoodContext(),
)
    return PointwiseLikelihoodContext{typeof(likelihoods),typeof(context)}(
        likelihoods, context
    )
end

NodeTrait(::PointwiseLikelihoodContext) = IsParent()
childcontext(context::PointwiseLikelihoodContext) = context.context
function setchildcontext(context::PointwiseLikelihoodContext, child)
    return PointwiseLikelihoodContext(context.loglikelihoods, child)
end

function Base.push!(
    context::PointwiseLikelihoodContext{Dict{VarName,Vector{Float64}}},
    vn::VarName,
    logp::Real,
)
    lookup = context.loglikelihoods
    ℓ = get!(lookup, vn, Float64[])
    return push!(ℓ, logp)
end

function Base.push!(
    context::PointwiseLikelihoodContext{Dict{VarName,Float64}}, vn::VarName, logp::Real
)
    return context.loglikelihoods[vn] = logp
end

function Base.push!(
    context::PointwiseLikelihoodContext{Dict{String,Vector{Float64}}},
    vn::VarName,
    logp::Real,
)
    lookup = context.loglikelihoods
    ℓ = get!(lookup, string(vn), Float64[])
    return push!(ℓ, logp)
end

function Base.push!(
    context::PointwiseLikelihoodContext{Dict{String,Float64}}, vn::VarName, logp::Real
)
    return context.loglikelihoods[string(vn)] = logp
end

function Base.push!(
    context::PointwiseLikelihoodContext{Dict{String,Vector{Float64}}},
    vn::String,
    logp::Real,
)
    lookup = context.loglikelihoods
    ℓ = get!(lookup, vn, Float64[])
    return push!(ℓ, logp)
end

function Base.push!(
    context::PointwiseLikelihoodContext{Dict{String,Float64}}, vn::String, logp::Real
)
    return context.loglikelihoods[vn] = logp
end

function tilde_observe!(context::PointwiseLikelihoodContext, right, left, vi)
    # Defer literal `observe` to child-context.
    return tilde_observe!(context.context, right, left, vi)
end
function tilde_observe!(context::PointwiseLikelihoodContext, right, left, vn, vi)
    # Need the `logp` value, so we cannot defer `acclogp!` to child-context, i.e.
    # we have to intercept the call to `tilde_observe!`.
    logp = tilde_observe(context.context, right, left, vi)
    acclogp!(vi, logp)

    # Track loglikelihood value.
    push!(context, vn, logp)

    return left
end

function dot_tilde_observe!(context::PointwiseLikelihoodContext, right, left, vi)
    # Defer literal `observe` to child-context.
    return dot_tilde_observe!(context.context, right, left, vi)
end
function dot_tilde_observe!(context::PointwiseLikelihoodContext, right, left, vn, vi)
    # Need the `logp` value, so we cannot defer `acclogp!` to child-context, i.e.
    # we have to intercept the call to `dot_tilde_observe!`.

    # We want to treat `.~` as a collection of independent observations,
    # hence we need the `logp` for each of them. Broadcasting the univariate
    # `tilde_obseve` does exactly this.
    logps = _pointwise_tilde_observe(context.context, right, left, vi)
    acclogp!(vi, sum(logps))

    # Need to unwrap the `vn`, i.e. get one `VarName` for each entry in `left`.
    _, _, vns = unwrap_right_left_vns(right, left, vn)
    for (vn, logp) in zip(vns, logps)
        # Track loglikelihood value.
        push!(context, vn, logp)
    end

    return left
end

# FIXME: This is really not a good approach since it needs to stay in sync with
# the `dot_assume` implementations, but as things are _right now_ this is the best we can do.
function _pointwise_tilde_observe(context, right, left, vi)
    return tilde_observe.(Ref(context), right, left, Ref(vi))
end

function _pointwise_tilde_observe(
    context, right::MultivariateDistribution, left::AbstractMatrix, vi
)
    return tilde_observe.(Ref(context), Ref(right), eachcol(left), Ref(vi))
end

"""
    pointwise_loglikelihoods(model::Model, chain::Chains, keytype = String)

Runs `model` on each sample in `chain` returning a `Dict{String, Matrix{Float64}}`
with keys corresponding to symbols of the observations, and values being matrices
of shape `(num_chains, num_samples)`.

`keytype` specifies what the type of the keys used in the returned `Dict` are.
Currently, only `String` and `VarName` are supported.

# Notes
Say `y` is a `Vector` of `n` i.i.d. `Normal(μ, σ)` variables, with `μ` and `σ`
both being `<:Real`. Then the *observe* (i.e. when the left-hand side is an
*observation*) statements can be implemented in three ways:
1. using a `for` loop:
```julia
for i in eachindex(y)
    y[i] ~ Normal(μ, σ)
end
```
2. using `.~`:
```julia
y .~ Normal(μ, σ)
```
3. using `MvNormal`:
```julia
y ~ MvNormal(fill(μ, n), σ^2 * I)
```

In (1) and (2), `y` will be treated as a collection of `n` i.i.d. 1-dimensional variables,
while in (3) `y` will be treated as a _single_ n-dimensional observation.

This is important to keep in mind, in particular if the computation is used
for downstream computations.

# Examples
## From chain
```julia-repl
julia> using DynamicPPL, Turing

julia> @model function demo(xs, y)
           s ~ InverseGamma(2, 3)
           m ~ Normal(0, √s)
           for i in eachindex(xs)
               xs[i] ~ Normal(m, √s)
           end

           y ~ Normal(m, √s)
       end
demo (generic function with 1 method)

julia> model = demo(randn(3), randn());

julia> chain = sample(model, MH(), 10);

julia> pointwise_loglikelihoods(model, chain)
Dict{String,Array{Float64,2}} with 4 entries:
  "xs[3]" => [-1.42862; -2.67573; … ; -1.66251; -1.66251]
  "xs[1]" => [-1.42932; -2.68123; … ; -1.66333; -1.66333]
  "xs[2]" => [-1.6724; -0.861339; … ; -1.62359; -1.62359]
  "y"     => [-1.51265; -0.914129; … ; -1.5499; -1.5499]

julia> pointwise_loglikelihoods(model, chain, String)
Dict{String,Array{Float64,2}} with 4 entries:
  "xs[3]" => [-1.42862; -2.67573; … ; -1.66251; -1.66251]
  "xs[1]" => [-1.42932; -2.68123; … ; -1.66333; -1.66333]
  "xs[2]" => [-1.6724; -0.861339; … ; -1.62359; -1.62359]
  "y"     => [-1.51265; -0.914129; … ; -1.5499; -1.5499]

julia> pointwise_loglikelihoods(model, chain, VarName)
Dict{VarName,Array{Float64,2}} with 4 entries:
  xs[2] => [-1.6724; -0.861339; … ; -1.62359; -1.62359]
  y     => [-1.51265; -0.914129; … ; -1.5499; -1.5499]
  xs[1] => [-1.42932; -2.68123; … ; -1.66333; -1.66333]
  xs[3] => [-1.42862; -2.67573; … ; -1.66251; -1.66251]
```

## Broadcasting
Note that `x .~ Dist()` will treat `x` as a collection of
_independent_ observations rather than as a single observation.

```jldoctest; setup = :(using Distributions)
julia> @model function demo(x)
           x .~ Normal()
       end;

julia> m = demo([1.0, ]);

julia> ℓ = pointwise_loglikelihoods(m, VarInfo(m)); first(ℓ[@varname(x[1])])
-1.4189385332046727

julia> m = demo([1.0; 1.0]);

julia> ℓ = pointwise_loglikelihoods(m, VarInfo(m)); first.((ℓ[@varname(x[1])], ℓ[@varname(x[2])]))
(-1.4189385332046727, -1.4189385332046727)
```

"""
function pointwise_loglikelihoods(model::Model, chain, keytype::Type{T}=String) where {T}
    # Get the data by executing the model once
    vi = VarInfo(model)
    context = PointwiseLikelihoodContext(Dict{T,Vector{Float64}}())

    iters = Iterators.product(1:size(chain, 1), 1:size(chain, 3))
    for (sample_idx, chain_idx) in iters
        # Update the values
        setval!(vi, chain, sample_idx, chain_idx)

        # Execute model
        model(vi, context)
    end

    niters = size(chain, 1)
    nchains = size(chain, 3)
    loglikelihoods = Dict(
        varname => reshape(logliks, niters, nchains) for
        (varname, logliks) in context.loglikelihoods
    )
    return loglikelihoods
end

function pointwise_loglikelihoods(model::Model, varinfo::AbstractVarInfo)
    context = PointwiseLikelihoodContext(Dict{VarName,Vector{Float64}}())
    model(varinfo, context)
    return context.loglikelihoods
end


text:
"""
$(DocStringExtensions.README)
"""
module SBML

using SBML_jll, Libdl

using DocStringExtensions
using IfElse
using SparseArrays
using Symbolics
using Unitful

include("types.jl")
include("structs.jl")
include("version.jl")

include("converters.jl")
include("math.jl")
include("readsbml.jl")
include("symbolics.jl")
include("utils.jl")

sbml(sym::Symbol)::VPtr = dlsym(SBML_jll.libsbml_handle, sym)

export readSBML, readSBMLFromString, stoichiometry_matrix, flux_bounds, flux_objective
export set_level_and_version, libsbml_convert, convert_simplify_math

end # module


text:
include("utils.jl")
include("gbm.jl")
include("disp.jl")
include("sampler.jl")

max_stars_repo_name: oilneck/InverseIsing

max_issues_repo_name: oilneck/InverseIsing

max_forks_repo_name: oilneck/InverseIsing


text:
const N_GRADS_3D_EXPONENT = 8
const N_GRADS_3D = 1 << N_GRADS_3D_EXPONENT

const NORMALIZER_3D = 0.2781926117527186

struct Grad3
    x::Float32
    y::Float32
    z::Float32
end

const grad3 = Grad3[
    Grad3( 2.22474487139,       2.22474487139,      -1.0),
    Grad3( 2.22474487139,       2.22474487139,       1.0),
    Grad3( 3.0862664687972017,  1.1721513422464978,  0.0),
    Grad3( 1.1721513422464978,  3.0862664687972017,  0.0),
    Grad3(-2.22474487139,       2.22474487139,      -1.0),
    Grad3(-2.22474487139,       2.22474487139,       1.0),
    Grad3(-1.1721513422464978,  3.0862664687972017,  0.0),
    Grad3(-3.0862664687972017,  1.1721513422464978,  0.0),
    Grad3(-1.0,                -2.22474487139,      -2.22474487139),
    Grad3( 1.0,                -2.22474487139,      -2.22474487139),
    Grad3( 0.0,                -3.0862664687972017, -1.1721513422464978),
    Grad3( 0.0,                -1.1721513422464978, -3.0862664687972017),
    Grad3(-1.0,                -2.22474487139,       2.22474487139),
    Grad3( 1.0,                -2.22474487139,       2.22474487139),
    Grad3( 0.0,                -1.1721513422464978,  3.0862664687972017),
    Grad3( 0.0,                -3.0862664687972017,  1.1721513422464978),
    #-------------------------------------------------------------------#
    Grad3(-2.22474487139,      -2.22474487139,      -1.0),
    Grad3(-2.22474487139,      -2.22474487139,       1.0),
    Grad3(-3.0862664687972017, -1.1721513422464978,  0.0),
    Grad3(-1.1721513422464978, -3.0862664687972017,  0.0),
    Grad3(-2.22474487139,      -1.0,                -2.22474487139),
    Grad3(-2.22474487139,       1.0,                -2.22474487139),
    Grad3(-1.1721513422464978,  0.0,                -3.0862664687972017),
    Grad3(-3.0862664687972017,  0.0,                -1.1721513422464978),
    Grad3(-2.22474487139,      -1.0,                 2.22474487139),
    Grad3(-2.22474487139,       1.0,                 2.22474487139),
    Grad3(-3.0862664687972017,  0.0,                 1.1721513422464978),
    Grad3(-1.1721513422464978,  0.0,                 3.0862664687972017),
    Grad3(-1.0,                 2.22474487139,      -2.22474487139),
    Grad3( 1.0,                 2.22474487139,      -2.22474487139),
    Grad3( 0.0,                 1.1721513422464978, -3.0862664687972017),
    Grad3( 0.0,                 3.0862664687972017, -1.1721513422464978),
    Grad3(-1.0,                 2.22474487139,       2.22474487139),
    Grad3( 1.0,                 2.22474487139,       2.22474487139),
    Grad3( 0.0,                 3.0862664687972017,  1.1721513422464978),
    Grad3( 0.0,                 1.1721513422464978,  3.0862664687972017),
    Grad3( 2.22474487139,      -2.22474487139,      -1.0),
    Grad3( 2.22474487139,      -2.22474487139,       1.0),
    Grad3( 1.1721513422464978, -3.0862664687972017,  0.0),
    Grad3( 3.0862664687972017, -1.1721513422464978,  0.0),
    Grad3( 2.22474487139,      -1.0,                -2.22474487139),
    Grad3( 2.22474487139,       1.0,                -2.22474487139),
    Grad3( 3.0862664687972017,  0.0,                -1.1721513422464978),
    Grad3( 1.1721513422464978,  0.0,                -3.0862664687972017),
    Grad3( 2.22474487139,      -1.0,                 2.22474487139),
    Grad3( 2.22474487139,       1.0,                 2.22474487139),
    Grad3( 1.1721513422464978,  0.0,                 3.0862664687972017),
    Grad3( 3.0862664687972017,  0.0,                 1.1721513422464978),
]

const GRADIENTS_3D = Vector{Grad3}(undef, N_GRADS_3D)

function grad(a, seed, xrvp, yrvp, zrvp, dx, dy, dz)
    hash = xor(seed, xrvp, yrvp, zrvp) * HASH_MULTIPLIER
    g = GRADIENTS_3D[(xor(hash, hash >> (64 - N_GRADS_3D_EXPONENT)) & (N_GRADS_3D - 1)) + 1]
    Float32(p4(a) * (g.x * dx + g.y * dy + g.z * dz) / NORMALIZER_3D)
end

function fill_3d()
    j = 0
    len = length(grad3)
    for i = 1:N_GRADS_3D
        GRADIENTS_3D[i] = grad3[j += 1]
        j == len && (j = 0)
    end
end

fill_3d()


text:
using Plots
data = readdlm("dsgedata.txt")
plot(data, label = ["y", "c", "n", "r", "w"], show=true)
savefig("dsgedata.png")

max_stars_repo_name: wenddymacro/Econometrics

max_issues_repo_name: wenddymacro/Econometrics

max_forks_repo_name: wenddymacro/Econometrics


text:
# Download Pluto, change CDNs to local file, but don't download files
using Pluto, DelimitedFiles

#-----------------------------------------------------------------------------# Copy Pluto files
@info "Copying Pluto source code..."
pluto_src = abspath(joinpath(pathof(Pluto), "..", ".."))

pluto_dest = joinpath(@__DIR__(), "Pluto")
@info "Copying $pluto_src -> $pluto_dest"
cp(pluto_src, pluto_dest, force=true)

#-----------------------------------------------------------------------------# make edits
@info "Replacing CDN links..."
function _replace(file, old_new::Pair)
    old, new = old_new
    sed = Sys.isapple() ?  ["sed", "-i", ""] : ["sed", "-i"]
    run(`$sed $("s|$old|$new|g") $file`)
end

frontend = joinpath(pluto_dest, "frontend")
assets = joinpath(frontend, "offline_assets")
mkpath(assets)

lines = split(read(`grep -r "jsdelivr" $frontend`, String), '\n')

replaced_assets = []

for line in lines
    match = findfirst(r"((?<=\")https:\/\/[^\"]*)|((?<=\()https:\/\/[^\)]*)", line)
    if isnothing(match)
        @debug "No match found for line: $line"
    else
        url = line[match]
        push!(replaced_assets, url)
        # @info "Downloading: $url"

        # file = touch(joinpath(assets, basename(url)))
        # try
        #     write(touch(file), read(`curl --silent /dev/null $url`))
        # catch
        #     @warn "$url didn't want to download..."
        # end
        pluto_file = line[1:findfirst(':', line) - 1]
        _replace(pluto_file, url => "/offline_assets/$(basename(url))")
    end
end

writedlm(joinpath(@__DIR__(), "assets.csv"), replaced_assets)

#-----------------------------------------------------------------------------# Fix project_relative_path
_replace(joinpath(pluto_dest, "src", "Pluto.jl"), "pathof(Pluto)" => """joinpath(@__DIR__, "..")""")

@info "Done"


text:
using Stheno: FiniteZeroMean, FiniteZeroKernel, LhsFiniteZeroCrossKernel,
    RhsFiniteZeroCrossKernel, FiniteZeroCrossKernel
using Base: OneTo

@testset "zero" begin

    let
        μ1, μ2 = ZeroMean{Float64}(), FiniteZeroMean(OneTo(10))
        μ3, μ4 = BlockMean([μ1, μ1]), BlockMean([μ2, μ2])
        k1, k2 = ZeroKernel{Float64}(), FiniteZeroCrossKernel(OneTo(4), OneTo(5))
        k3, k4 = LhsFiniteZeroCrossKernel(OneTo(3)), RhsFiniteZeroCrossKernel(OneTo(5))
        k5, k6 = ZeroKernel{Float64}(), FiniteZeroKernel(OneTo(7))

        k7, k8 = BlockCrossKernel([k1 k1; k1 k1]), BlockCrossKernel([k2 k2; k2 k2])
        k9, k10 = BlockCrossKernel([k3 k3; k3 k3]), BlockCrossKernel([k4 k4; k4 k4])

        k11 = BlockKernel([k5, k5], [k5 k5; k5 k5])
        k12 = BlockKernel([k6, k6], [k6 k6; k6 k6])

        # Check addition of zero elements.
        for x in [μ1, μ2, μ3, μ4, k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12]
            @test zero(x) == x
            @test x + x == zero(x)
            @test x * x == x
        end

        # Check (element-wise) multiplication of zero elements.
        @test_throws AssertionError μ2 * FiniteZeroMean(randn(11))
        @test_throws AssertionError k2 * FiniteZeroCrossKernel(randn(5), randn(1))
        @test_throws AssertionError k3 * LhsFiniteZeroCrossKernel(randn(4))
        @test_throws AssertionError k4 * RhsFiniteZeroCrossKernel(randn(2))
        @test_throws AssertionError k6 * FiniteZeroKernel(randn(8))

        @test μ1 * ConstantMean(3.0) === μ1
        @test ConstantMean(2.0) * μ1 === μ1

        @test μ2 * FiniteMean(μ1, randn(10)) === μ2
        @test FiniteMean(μ1, randn(10)) * μ2 === μ2

        @test k1 * EQ() === k1
        @test EQ() * k1 === k1

    end

end


text:
using MiCRM
using Test

include("params_test.jl")
include("MTK_test.jl")
include("DiffEq_test.jl")

# @testset "MiCRM.jl" begin

# end


text:
using CoupledDipoles, Revise, Plots, DelimitedFiles

### --- INITIAL CONDITIONS ARE CHANGED TO AGREE WITH QUTIP SAVED DATA ---
function get_initial_condition(problem::LinearOptics{Scalar})
    return zeros(ComplexF64, problem.atoms.N)
end

function get_initial_condition(problem::NonLinearOptics{MeanField})
    β₀ = zeros(ComplexF64, problem.atoms.N)
    z₀ = -ones(problem.atoms.N)
    u₀ = vcat(β₀, z₀)
    return u₀
end

### ------------ FIXED INPUTS ---------------------
r = [ -1.01015   -0.455099     0.543082;
      -2.96935   -0.0781967   -0.00764304;
      -0.714778  -0.248556     0.683987;
      -1.20512    0.376068    -1.85059;
      -1.89765    0.84797      0.248158;
      0.664911   0.00503915  -0.177228]

atoms = Atom(Cube(), Array(transpose(r)), 10);

Δ = -2.0
s = 0.0022
laser = Laser(PlaneWave3D([0,0,1]), s, Δ);

### ------------ CREATE PROBLEM AND EVOLVE ---------------------
simulationScalar =  LinearOptics(Scalar(), atoms, laser)
simulationMeanField = NonLinearOptics(MeanField(), atoms, laser)

u0_scalar    = get_initial_condition(simulationScalar)
u0_meanfield = get_initial_condition(simulationMeanField)
tspan = (0.0,100)

aS = time_evolution(simulationScalar, u0_scalar, tspan);
aMF = time_evolution(simulationMeanField, u0_meanfield, tspan);

### ------------ GET INTENSITIES ---------------------
iS = zeros( length(aS.t) );
for t ∈ eachindex(iS)
    iS[t] = get_intensity_over_an_angle(simulationScalar, aS.u[t], deg2rad(35))
end

iMF = zeros( length(aMF.t) );
for t ∈ eachindex(iMF)
    iMF[t] = get_intensity_over_an_angle(simulationMeanField, aMF.u[t], deg2rad(35))
end

matrixData = readdlm(pwd()*"/test/benchmarks/Field_qutip.txt", '\t', Float64, '\n')
time_Qutip = matrixData[:,1]
iQutip = matrixData[:,2]

### ------------ STEADY STATES -----------------------
ss_S = get_steady_state(simulationScalar)
ss_MF = get_steady_state(simulationMeanField; time_max=1000)

i_ss_S = get_intensity_over_an_angle(simulationScalar, ss_S, deg2rad(35))
i_ss_MF =get_intensity_over_an_angle(simulationMeanField, ss_MF , deg2rad(35))

### ------------ PLOTS ---------------------
plot(aS.t, iS, label="Scalar", lw=5, ylabel="intensity", xlabel="time")
plot!(aMF.t, iMF, label="Mean Field", lw=4, linestyle=:dash, size=(800,600))
plot!( time_Qutip, iQutip, label="Qutip", lw=4, linestyle=:dot  )
plot!(guidefont=17, tickfont=15, legendfontsize=15, legend=:bottomright)

hline!([i_ss_S], label="i_ss_S", lw=5, linestyle=:dash)
hline!([i_ss_MF], label="i_ss_MF", lw=5, linestyle=:dot)


text:
#
function permutedims(A::AbstractArray, perm)
    dest = similar(A, genperm(axes(A), perm))
    permutedims!(dest, A, perm)
end

"""
    permutedims(m::AbstractMatrix)
Permute the dimensions of the matrix `m`, by flipping the elements across the diagonal of
the matrix. Differs from `LinearAlgebra`'s [`transpose`](@ref) in that the
operation is not recursive.
# Examples
```jldoctest; setup = :(using LinearAlgebra)
julia> a = [1 2; 3 4];
julia> b = [5 6; 7 8];
julia> c = [9 10; 11 12];
julia> d = [13 14; 15 16];
julia> X = [[a] [b]; [c] [d]]
2×2 Matrix{Matrix{Int64}}:
 [1 2; 3 4]     [5 6; 7 8]
 [9 10; 11 12]  [13 14; 15 16]
julia> permutedims(X)
2×2 Matrix{Matrix{Int64}}:
 [1 2; 3 4]  [9 10; 11 12]
 [5 6; 7 8]  [13 14; 15 16]
julia> transpose(X)
2×2 transpose(::Matrix{Matrix{Int64}}) with eltype Transpose{Int64, Matrix{Int64}}:
 [1 3; 2 4]  [9 11; 10 12]
 [5 7; 6 8]  [13 15; 14 16]
```
"""
permutedims(A::AbstractMatrix) = permutedims(A, (2,1))

"""
    permutedims(v::AbstractVector)
Reshape vector `v` into a `1 × length(v)` row matrix.
Differs from `LinearAlgebra`'s [`transpose`](@ref) in that
the operation is not recursive.
# Examples
```jldoctest; setup = :(using LinearAlgebra)

 """

permutedims([1, 2, 3, 4])
"""1×4 Matrix{Int64}:
 1  2  3  4"""

 V = [[[1 2; 3 4]]; [[5 6; 7 8]]]

 """2-element Vector{Matrix{Int64}}:
 [1 2; 3 4]
 [5 6; 7 8]
"""

permutedims(V)
"""1×2 Matrix{Matrix{Int64}}:
 [1 2; 3 4]  [5 6; 7 8]"""

transpose(V)

"""1×2 transpose(::Vector{Matrix{Int64}}) with eltype Transpose{Int64, Matrix{Int64}}:
 [1 3; 2 4]  [5 7; 6 8]"""

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
export delim_vector

"""
    Use: delim_vector(x::AbstractVector, probe_length::Int)

Return the indicies of the first element for each homogenous section in `x`.

# Example
```
julia> x = vec([1.0 1.0 1.0 2.0 2.0 2.0 2.0 3.0 3.0 3.0 3.0]);

julia> delim_vector(x, 1)
3-element Array{Int64,1}:
1
4
8
```

`probe_length` sets the intial guess for the section length. Choosing `probe_length` close the the true length of the section provides modest performance benefits, however simply setting the length to 1 is not a huge hit.

```
julia> @btime delim_vector(x,1)
  1.169 μs (3 allocations: 176 bytes)
4-element Array{Int64,1}:
      1
 100001
 100003
 200003

julia> @btime delim_vector(x,100000)
  945.615 ns (3 allocations: 176 bytes)
4-element Array{Int64,1}:
      1
 100001
 100003
 200003
```
"""
function delim_vector(x::AbstractVector, probe_length::Int)

    delims = Array{Int,1}(undef, 1)
    delims[1] = 1
    n = length(x)

    while true
        j = find_next_delim(x, delims[end], probe_length)
        (j < n) ? (push!(delims, j)) : (return delims)
    end

end


text:
# TODO: rewrite in a better and more julian way.
#
# This is just kind of thrown together without any real planning but it seems to
# work. An obviouse design improvement would be to drop the string identifiers
# and use types or something instead. Conssider this as a functional outline for
# a future product.
s = raw"""
@blue:#5555ff

.round-button{
    display:inline-block; margin:7; radius:12; color: @blue;
    hover{color: @blue; }
    circle{ radius:12;}
}

circle{ class:round-button; circle{ image:"Back.png"} }
"""

mutable struct Vars
    name::String
    Vars(name) = new(name)
end
# include("Documents/code/parseScss.jl")
#-==============================================================================
# Tokenize Sml text
#-==============================================================================
function tokens(s)
	toks = []
from,i = 1,1
while i < length(s)
    from = i
        if s[from] == '/'
    # Comment-------------------------
            if s[i+1] == '*'
                i+=1
                while !(s[i] == '/' && s[i-1] == '*')
                    i+=1
                end
				push!(toks, ["Comment", s[from:i]])
    # Comment-------------------------
            elseif s[i+1] == '/'
	              i+=1
                while s[i+1] != '\n'
                    i+=1
                end
				push!(toks, ["Comment", s[from:i]])
    # Divide-------------------------
            else
				push!(toks, ["Math", s[from:i]])
            end
    # Var-------------------------
		elseif s[from] == '@'
            i+=1
            while isletter(s[i+1]) || s[i+1] == '-'
                i+=1
            end
			push!(toks, ["Vars", Vars(s[from+1:i])])
	# Text-------------------------
		elseif s[from] == '"'
			if s[i+1] == '"' &&  s[i+2] == '"'
				i+=3
				while !(s[i] == '"' && s[i-1] == '"' && s[i-2] == '"' && s[i-3] != '\\')
					 i+=1
				end
				push!(toks, ["Text", s[from:i]])
			else
			     i+=1
			     while s[i] != '"' && !(s[i-1] == '\\' && s[i] == '"')
			          i+=1
			     end
				 push!(toks, ["Text", s[from:i]])
		 	end
    # Identifier-------------------------
        elseif isletter(s[from])
            i+=1
            while isletter(s[i+1]) || s[i+1] == '-'
            i+=1
        end
		push!(toks, ["Identifier", s[from:i]])
    # Class-------------------------
        elseif s[from] == '.'
            i+=1
            while isletter(s[i+1]) || s[i+1] == '-'
                i+=1
            end
		push!(toks, ["Class", s[from+1:i]])
	# Number------------------------
		elseif s[from] == '-' || isdigit(s[i]) || s[i] == '.'
	    	while isdigit(s[i+1]) || s[i+1] == '.'
		        i+=1
	        end
			push!(toks, ["Number", parse(Float64, s[from:i])])
			if s[i+1] == '%'
				i+=1
				push!(toks, ["Unit", s[i]])
			end
	# Math------------------------
	        elseif s[from] == '+' || s[i] == '-' || s[i] == '/' || s[i] == '*' || s[i] == '%'
				push!(toks, ["Math", s[i]])
    # Hex-------------------------
        elseif s[from] == '#'
            i+=1
            while isdigit(s[i+1]) || 'a' <= s[i+1] <= 'f'
                i+=1
            end
			push!(toks, ["Hex", s[from:i]])
    # Space-------------------------
        elseif s[from] == ' ' || s[i] == '\t'
            while s[i+1] == ' ' || s[i+1] == '\t'
                i+=1
            end
			push!(toks, ["Space", s[from:i]])
    # Comma-------------------------
        elseif s[from] == ','
			push!(toks, ["Comma", s[i]])
    # Instance-------------------------
        elseif s[from] == ':'
			push!(toks, ["Instance", s[i]])
    # OpenBl-------------------------
        elseif s[from] == '{'
			push!(toks, ["OpenBl", s[i]])
    # CloseBl-------------------------
        elseif s[from] == '}'
			push!(toks, ["CloseBl", s[i]])
    # EndLine-------------------------
        elseif s[from] == '\n' ||  s[i] == ';'
			push!(toks, ["EndLine", s[i]])
    # CatchAll-------------------------
        else
             println("CatchAll" * s[i])
        end
    i+=1
  end
  return toks
end
#-==============================================================================
# "Class" "Identifier" "Vars" "Comment" "Comma" "Instance"  "OpenBl" "CloseBl"
#-==============================================================================
function crunch(toks)
	stack = []
	push!(toks, ["EndLine",';'], ["EndLine",';']);

# Remove whitspace and other nodes
push!(stack, toks[1])
	for t in 2:length(toks)
		if  !((toks[t][1] == "EndLine" && toks[t-1][1] == "EndLine") || toks[t][1] == "Space" || toks[t][1] == "Comma" || toks[t][1] == "Comment")
		   push!(stack, toks[t])
		end
	end

	# find Attributes
	toks = stack
	stack = []
	t=1
	push!(toks, ["telemer",';'],["telemer",';'],["telemer",';'],["telemer",';'],["telemer",';']);
	while t < length(toks)
		if  toks[t][1] == "Identifier" && toks[t+1][1] == "Instance"
			toks[t][1] = "Attribute"
			push!(stack, toks[t])
			t+=2
		else
			push!(stack, toks[t])
			t+=1
		end
	end

	# Make numbers
	toks = stack
	stack = []
	t=1
	while t <= length(toks)
		if  toks[t][1] == "Number" && (toks[t+1][1] == "Identifier" || toks[t+1][1] == "Unit") # Unit
			push!(stack, ["Identifier", [toks[t][2], toks[t+1][2]]])
			t+=2
		else
			push!(stack, toks[t])
			t+=1
		end
	end

	# Make Arrays      {[\n\s]+}    [^{}]+
	# Hex Number Text Var Identifier
	toks = stack
	stack = []
	t=1
	while t < length(toks)
		b, c = toks[t][1], toks[t+1][1]
		if (b == "Hex" || b == "Number" || b == "Text" || b == "Vars" || b == "Identifier") && (c == "Hex" || c == "Number" || c == "Text" || c == "Vars" || c == "Identifier")
			a = [toks[t][2], toks[t+1][2]]
			t+=2
			while t < length(toks) && (toks[t][1] == "Hex" || toks[t][1] == "Number" || toks[t][1] == "Text" || toks[t][1] == "Vars" || toks[t][1] == "Identifier")
				push!(a, toks[t][2])
				t+=1
			end
			push!(stack, ["Identifier", a])
		else
			push!(stack, toks[t])
			t+=1
		end
	end

t = stack
push!(t, ["EndLine",';'], ["EndLine",';'], ["EndLine",';']);

stack = []
i=1
	while i < length(t)-3
		one, two, three, four = t[i],t[i+1],t[i+2],t[i+3]
		if  one[1] == "Vars" && two[1] == "Instance"
			if three[1] == "Identifier" || three[1] == "Hex" || three[1] == "Vars"
				push!(stack, [one[1], one[2], "Identifier", three[2]])
				i+=3
			end
		# Attribute
		elseif  one[1] == "Attribute" #&& two[1] == "Identifier"
			if two[1] == "Vars" || two[1] == "Text" || two[1] == "Hex" || two[1] == "Identifier" || two[1] == "Number"
				push!(stack, [one[1], one[2], two[1], two[2]])
			i+=2
			end
		elseif  one[1] == "CloseBl"
			push!(stack, one)
			i+=1
		elseif  one[1] == "EndLine" || one[1] == "Space"
			i+=1
		elseif  one[1] == "Identifier" && two[1] == "OpenBl"
			one[1] = "Tag"
			push!(stack, one)
			push!(stack, two)
			i+=2
		else
			push!(stack, one)
			i+=1
		end
	end

	return stack
end
#-==============================================================================
# EndLine CloseBl OpenBl Instance Comma Space Hex Math Unit Number Class Identifier Text Vars Comment
# Identifier
# Text Hex Math
# Class Vars Tag
#-==============================================================================

#-==================================OUTLINE=====================================
#                                               # structure, content, styles
# element: ["div", Dict(atributes), []]         #    yes       yes     yes
# variable: ["MyTemplate", Dict(atributes), []] #    yes     default   yes
# variable: ["Mycolor", Value]                  # ---- just store a value ----
# class: ["div", Dict(atributes), []]           #    no        no      yes
#-==============================================================================
#  [ element/variable/class(name: , Dict(atributes), []) ]
function nest(toks)
	this = Dict("styles"=>Dict(), "templates"=>Dict(), "nodes"=>[])
	t = 1

	function element(toks, named=true)
		name = toks[t][2]
		node = Dict()
		node["nodes"] = []
		# Dict( ">"=>"div", "height"=>value(39), "nodes"=>[] )
		t+=1  #["Vars", "template", []] & ["OpenBl", '{']
			while toks[t][1] != "CloseBl"
				if toks[t+1][1] == "OpenBl"
					push!(node["nodes"], element(toks))
				elseif toks[t][1] == "Identifier" || toks[t][1] == "Text"
					println(toks[t])
				elseif toks[t][1] == "Attribute" || toks[t][1] == "Text"
					node[toks[t][2]] = toks[t][4]
				elseif toks[t][1] == "OpenBl"
					;
				else
					println("Error!", toks[t])
				end
				t+=1
			end
			if named==true
				node[">"] = name
			end
				return node #
	end

	while t < length(toks)
		if toks[t][1] == "Vars"
			if toks[t+1][1] == "OpenBl"  #["Vars", "template", []] ["OpenBl", '{']
				push!(this["templates"], toks[t][2] => element(toks, false))
			else
				push!(this["templates"], toks[t][2] => toks[t][4]) # Plain vanilla Variable
			end
		elseif toks[t][1] == "Tag"
			push!(this["nodes"], element(toks))
		elseif toks[t][1] == "Class"
			push!(this["styles"], toks[t][2] => element(toks, false))
			println(this["styles"])
		else
			# println(toks[t])
		end
		t+=1
	end
	return this
end

function readSml(text)
	return nest( crunch(tokens(text)) )
end

function writeSml(text)
	result = " "
	# TODO: write code!
	return result
end

function addAttributes(node, style)  classname, stylesDict

		attrs = collect(keys(style))
		for a in attrs
			if a == "nodes"
				for n in node["nodes"]
					for st in style["nodes"]
						if n[">"] == st[">"]
							addAttributes(n, st)
						end
					end
				end
			elseif !haskey(node, a) && a != ">"
				node[a] = style[a]
			end
		end
end

function compile(nodes)
	for st in nodes
		if haskey(st, "class")
			# WARNING: class could be an array!
			# println(st["class"])
			if haskey(text["styles"], st["class"])
				style = text["styles"][st["class"]]
				addAttributes(st, style)
			end
		end
		if length(st["nodes"]) > 0
			compile(st["nodes"])
			println("has kids!")
		end
	end
end

text = readSml(s)
text["styles"]
compile(text["nodes"])

print(text["nodes"])


text:
using DataFrames, DataValues, Query, CategoricalArrays, Dates, Printf, XLSX
import Base.:(==), Base.:(<)

macro XLSX2DF(var, xlsx)
    quote
        $(esc(var)) = XLSX.openxlsx(string("../data/", $xlsx, ".xlsx")) do xf DataFrame(XLSX.gettable(xf[$xlsx])...) end
    end
end

mutable struct Balance{T <: Number}
    i::AbstractArray{T, 1}
    o::AbstractArray{T, 1}
    s::T
    f::T
    function Balance{T}(v::T) where {T <: Number}
        self = new()
        self.i = Array{T, 1}()
        self.o = Array{T, 1}()
        self.s = v
        self.f = v
        self
    end
end

Balance(v::Float64) = Balance{Float64}(v)
Balance(v::Int64) = Balance{Float64}(convert(Float64, v))

mutable struct DailyBalance{T <: Number}
    q::Balance{T}
    v::Balance{T}
    function DailyBalance{T}(q::T, v::T) where {T <: Number}
        self = new()
        self.q = Balance{T}(q)
        self.v = Balance{T}(v)
        self
    end
end

DailyBalance(q::Float64, v::Float64) = DailyBalance{Float64}(q, v)
DailyBalance(q::Int64, v::Int64) = DailyBalance{Float64}(convert(Float64, q), convert(Float64, v))
DailyBalance(q::Float64, v::Int64) = DailyBalance{Float64}(q, convert(Float64, v))
DailyBalance(q::Int64, v::Float64) = DailyBalance{Float64}(convert(Float64, q), v)

struct Key
    item::AbstractString
    date::Date
end

Key(item::AbstractString, date::DataValue{Any}) = Key(item, get(date))
Key(item::DataValue{Any}, date::DataValue{Any}) = Key(get(item), get(date))
==(k1::Key, k2::Key) = k1.item == k2.item && k1.date == k2.date
<(k1::Key, k2::Key) = k1.item < k2.item || (k1.item == k2.item && k1.date < k2.date)

function qsa(a::AbstractArray{T, 1}) where {T <: Number}
    v = join(string.(sort(a)), " ")
    length(a) > 1 ? string('"', v, '"') : v
end

function prt(i::AbstractString, d::Any, b::DailyBalance)
    println(join((i, Dates.format(d, "dd/mm/yyyy"),
        qsa(b.q.i), qsa(b.v.i), qsa(b.q.o), qsa(b.v.o),
        @sprintf("%.2f", b.q.s), @sprintf("%.2f", b.v.s),
        @sprintf("%.2f", b.q.f), @sprintf("%.2f", b.v.f)), ",")
    )
end

function calc!(bal::DailyBalance, tipo::AbstractString, qtd::Number, vlr::Number)
    if tipo == "Ent"
        push!(bal.q.i, qtd); push!(bal.v.i, vlr)
        bal.q.f += qtd; bal.v.f += vlr
    elseif tipo == "Sai"
        push!(bal.q.o, qtd); push!(bal.v.o, vlr)
        bal.q.f -= qtd; bal.v.f -= vlr
    end
end

calc!(bal::DailyBalance, tipo::DataValue{Any}, qtd::DataValue{Any}, vlr::DataValue{Any}) = calc!(bal, get(tipo), get(qtd), get(vlr))

@XLSX2DF sld "SaldoITEM"
@XLSX2DF mov "MovtoITEM"

data_min = minimum(sld.data_inicio)
balance = @from s in eachrow(sld) begin
    @where s.data_inicio == data_min
    @select Key(s.item, s.data_inicio) => DailyBalance(s.qtd_inicio, s.valor_inicio)
    @collect Dict{Key, DailyBalance}
end

last_dt = data_min

items = CategoricalArrays.levels(mov.item); sort!(items)

for item in items
    global last_dt = data_min
    for i in (mov |> @filter(_.item == item) |> @orderby(_.data_lancamento) |> collect)
        key = Key(i.item, i.data_lancamento)
        if i.data_lancamento != last_dt
            last_key = Key(item, last_dt)
            last_balance = balance[last_key]
            balance[key] = DailyBalance(last_balance.q.f, last_balance.v.f)
        end
        calc!(balance[key], i.tipo_movimento, i.quantidade, i.valor)
        global last_dt = i.data_lancamento
    end
end

for date in sort(CategoricalArrays.levels(mov.data_lancamento))
    for item in items
        if haskey(balance, Key(item, date))
            prt(item, date, balance[Key(item, date)])
        end
    end
end

# eof

max_stars_repo_name: msausu/micro-analytics-cmp

max_issues_repo_name: msausu/micro-analytics-cmp

max_forks_repo_name: msausu/micro-analytics-cmp


text:
module Vehicles

# 1. Export/Imports
export go!

# 2. Interface documentation
# A vehicle (v) must implement the following functions:

# power_on!(v) - turn on the vehicle's engines
# power_off!(v) = turn off the vehicle's engines
# turn!(v, direction) - steer the behicle to the specifed direction
# move!(v, distance) - move the vehicle by the specified distance
# position!(v) - returns the (x, y) position of the vehicle
# engage_wheels!(v) - engage wheels for landing. Optional.

# 3. Generic definitions for the interface
# hard contracts
# ...
#soft contracts
engage_wheels!(args...) = nothing
function power_on! end
function power_off! end
function turn! end
function move! end
function position! end
# traits
has_wheels(vehicle) = error("Not implemented.")
# 4. Game logic
# Returns a travel plan from current position to destination
function travel_path(position, destination)
    return round(pi/6, digits=2), 1000 # just a test
end

# Space travel logic
function go!(vehicle, destination)
    power_on!(vehicle)
    direction, distance = travel_path(position(vehicle), destination)
    turn!(vehicle, direction)
    move!(vehicle, distance)
    power_off!(vehicle)
    nothing
end # function

# Landing functions
function land!(vehicle)
    engage_wheels!(vehicle)
    println("Landing vehicle: ", vehicle)
end
# Landing (using trait)
function land2!(vehicle)
    has_wheels(vehicle) && engage_wheels!(vehicle)
    println("Landing vehicle: ", vehicle)
end

end # module