text:
"""
    calculate_prf_response - dot product of stimImage with
    proposed pRF over many time frames

    ported from matlab code with a julian twist...

    ds 2020-10-13

"""
function calculate_prf_response(stimImUnwrapped, pRFhat, hrf=[])

# idea: consider image as a 1d array!
# "vectorise" stimImage and also the pRF
    # but avoid allocation by not creating variable...?!
    # pRFhatUnwrapped = vec(pRFhat);

    # calculate response with one matrix product
    r = stimImUnwrapped' * vec(pRFhat);

    # if the user asks for the response w/ hrf included, calculate it
    # this will save time, if not!
    if isempty(hrf)
        return r
    else
        # r_hrf = conv(r, hrf)
        r_hrf = ImageFiltering.imfilter(r, reflect(hrf)) # CONV, not CORR!
        r_hrf .-= mean(r_hrf)
        return r_hrf[1:size(stimImUnwrapped, 2)]
    end
end

"""
mgauss - 2d gaussian (x,y, with params p x0,y0, sdxy)
   inputs: p (parameter vector, [x0,y0, sdxy] )
           x,y (meshgrid format)

          NB! params first, to make optimisation easier

ds 2019-03-20 (break out of mrTools/mgl version)
ds 2020-10-13 julia version
"""
function mgauss(p, x, y)

# unpack parameters for readability
    x0 = p[1];
    y0 = p[2];
    sdxy = p[3];

# compute 2d gaussian
    m = exp.(-(((x .- x0).^2) ./ (2 * (sdxy.^2)) .+
           ((y .- y0).^2) ./ (2 * (sdxy.^2))));

# clamp small values to 0
    m[m .< 0.01] .= 0;

# return m
    m
end

"""
fmribHRF - hrf function used in by FMRIB / Oxford crew

Simple port of the Matlab code used elsewhere

    inputs: t, p [vector of params], scaled [default=1]
    outputs: result
    purpose: calculate a haemodynamic response function according
    to the equations given in the Jezzard et al book (page 252).

               H(t) = ( t / d1 ).^a1 .* exp( -(t-d1)./b1 )
                      .- c*( t / d2 ).^a2 .* exp( -(t-d2)./b2 )

   where di = ai*bi (default params a1, a2, b1, b2, c [6 12 0.9 0.9 0.35] )
"""
function fmribHRF(t=0.0:1.5:30.0, p=[6.0; 12.0; 0.9; 0.9; 0.35], scaled=true)

# if no parameters are passed use defaults from worsley chapter
# and a TR of 1.5s - which often used by Schluppeck et al.
#
# FMRIB book (p252)
# default params: a1, a2, b1, b2, c

# calculate the HRF
    t = collect(t);
    result = ( t ./ ( p[1] * p[3] )  ).^p[1] .* exp.(-( t .- p[1] * p[3] ) ./ p[3]) .-
    p[5] * ( t ./ ( p[2] * p[4] )  ).^p[2] .* exp.(-( t .- p[2] * p[4] ) ./ p[4]);

    if scaled
  # normalize to SUM=1
        result = result ./ sum(result);
    end
end

max_stars_repo_name: schluppeck/prf-explainer

max_issues_repo_name: schluppeck/prf-explainer

max_forks_repo_name: schluppeck/prf-explainer


text:
   using ResTime, PDBTools
   dir="/home/viniciusp/Documents/doutorado/ANALYSE/Repository/restime_tests/md-pulchra/urea/0.5/0"

 # dir="/home/viniciusp/Documents/doutorado/ANALYSE/Repository/restime_tests/Ionic_liquid"
  cd(dir)
  solute="protein"
  # solute (protein) selection
  atoms   = PDBTools.readPDB("processed.pdb")
  protein = PDBTools.select(atoms,solute)
  solute  = ResTime.Selection(protein,nmols=1) # next step, create a function to calculate these things for electrolyte solutions w/out proteins

#  sel            = PDBTools.select(atoms,"name C4 and resname UR")
  sel            = PDBTools.select(atoms,"name OW and resname SOL")
  solvent        = ResTime.Selection(sel,natomspermol=1)
  trajectory     = ResTime.Trajectory("processed.xtc",solute,solvent,format="XTC")

  time, prob = correlation_cell(trajectory,8.)
 # dist, domain, evals, sp, nprot, nsvt, nframes, stime = ]
 #autocorr_cell(trajectory)
 # ResTime.domain_eval_cell!(dist,domain,4.)

  # standard method
   trajectory     = ResTime.Trajectory("processed.xtc",solute,solvent,format="XTC")
   time2, prob2 = correlation(trajectory,8.)

  # dist2, domain2, evals2, sp2, nprot2, nsvt2, nframes2, stime2 = autocorr(trajectory)
  # ResTime.domain_eval!(dist2,domain2,4.)

   using Plots
   plot(time, prob, label="CellList", color="red", linestyle=:dot)
   plot!(time2, prob2, label="Standard", color="blue", linewidth=2)
   savefig("teste_water.png")

# AQUI TEM UM PROBLEMA COM AS MATRIZES DOMAIN

  if   prob ≈ prob2
    println("deu tudo certo no cálculo")
  else
    println("as probabilidades não estão iguais")
  end

ext: jl


text:
module ODEXform
using DifferentialEquations
using Cassette
using SemanticModels.Dubstep

Cassette.@context SolverCtx
function Cassette.overdub(ctx::SolverCtx, args...)
    if Cassette.canrecurse(ctx, args...)
        #newctx = Cassette.similarcontext(ctx, metadata = ctx.metadata)
        return Cassette.recurse(ctx, args...)
    else
        return Cassette.fallback(ctx, args...)
    end
end

function Cassette.overdub(ctx::SolverCtx, f::typeof(Base.vect), args...)
    @info "constructing a vector length $(length(args))"
    return Cassette.fallback(ctx, f, args...)
end

# We don't need to overdub basic math. this hopefully makes execution faster.
# if these overloads don't actually make it faster, they can be deleted.
function Cassette.overdub(ctx::SolverCtx, f::typeof(+), args...)
    return Cassette.fallback(ctx, f, args...)
end
function Cassette.overdub(ctx::SolverCtx, f::typeof(-), args...)
    return Cassette.fallback(ctx, f, args...)
end
function Cassette.overdub(ctx::SolverCtx, f::typeof(*), args...)
    return Cassette.fallback(ctx, f, args...)
end
function Cassette.overdub(ctx::SolverCtx, f::typeof(/), args...)
    return Cassette.fallback(ctx, f, args...)
end

end #module

using LinearAlgebra
using Test
using Cassette
using DifferentialEquations
using SemanticModels.Dubstep

"""   sir_ode(du,u,p,t)

computes the du/dt array for the SIR system. parameters p is b,g = beta,gamma.
"""
sir_ode(du,u,p,t) = begin
    S,I,R = u
    b,g = p
    du[1] = -b*S*I
    du[2] = b*S*I-g*I
    du[3] = g*I
end

function Cassette.overdub(ctx::ODEXform.SolverCtx, f::typeof(sir_ode), args...)
    y = Cassette.fallback(ctx, f, args...)
    # add a lagniappe of infection
    extra = args[1][1] * ctx.metadata.factor
    push!(ctx.metadata.extras, extra)
    args[1][1] += extra
    args[1][2] -= extra
    return y
end

function g()
    parms = [0.1,0.05]
    init = [0.99,0.01,0.0]
    tspan = (0.0,200.0)
    sir_prob = Dubstep.construct(ODEProblem,sir_ode,init,tspan,parms)
    return sir_prob
end

function h()
    prob = g()
    return solve(prob, alg=Vern7())
end

#precompile
@time sol1 = h()
#timeit
@time sol1 = h()

"""    perturb(f, factor)

run the function f with a perturbation specified by factor.
"""
function perturb(f, factor)
    t = (factor=factor,extras=Float64[])
    ctx = ODEXform.SolverCtx(metadata = t)
    val = Cassette.recurse(ctx, f)
    return val, t
end

traces = Any[]
solns = Any[]
for f in [0.0, 0.01, 0.05, 0.10]
    val, t = perturb(h, f)
    push!(traces, t)
    push!(solns, val)
end

for (i, s) in enumerate(solns)
    @show s(100)
    @show traces[i].factor
    @show traces[i].extras[5]
    @show sum(traces[i].extras)/length(traces[i].extras)
end

@testset "ODE perturbation" begin

@test norm(sol1(100) .- solns[1](100),2) < 1e-6
@test norm(sol1(100) .- solns[2](100),2) > 1e-6
@test norm(solns[1](100) .- solns[2](100),2) < norm(solns[1](100) .- solns[3](100),2)
@test norm(solns[1](100) .- solns[2](100),2) < norm(solns[1](100) .- solns[4](100),2)

end


text:
using Aurora
using StatsBase
Z = ReplicatedSample(randn(10))
@test collect(Z.Z) ==  Z.sorted_Z[Z.rank_idx]

μs = randn(1000)  # generate true
zs = sqrt(10) .* randn(1000, 10) .+ μs

Zs1 = ReplicatedSample.(copy.(eachrow(zs)))
Zs2 = ReplicatedSample.(zs)

@test mean.(Zs1) == mean.(Zs2)

Aurora.design_matrix(Zs1, 1)

empty_mat = zeros(9, 1000)
empty_vec = zeros(9)
empty_vec2 = zeros(11)

for j in Base.OneTo(10)
    @test reduce(hcat, sort.(deleteat!.(copy.(getfield.(Zs1, :Z)), j))) == Aurora.design_matrix(Zs1, j)
    @test reduce(hcat, sort.(deleteat!.(copy.(eachrow(zs)), j))) == Aurora.design_matrix(Zs1, j)
    Aurora.design_matrix!(empty_mat, Zs1, j)
    @test Aurora.design_matrix(Zs1, j) == empty_mat
    i = sample(1:1000)
    Aurora.design_row!(empty_vec, Zs1[i], j)
    @test sort(deleteat!(copy(Zs1[i].Z), j)) == empty_vec

    Aurora.design_response_row!(empty_vec2, Zs1[i], j)
    @test [1; sort(deleteat!(copy(Zs1[i].Z), j)); Zs1[i].Z[j]] == empty_vec2
    @test Zs1[i].Z[j] == StatsBase.response(Zs1[i], j)
end


text:
# code for Video 5.4:
julia> using RDatasets
julia> iris = dataset("datasets", "iris")
julia> features = convert(Array, iris[:, 1:4]) # measurements of length
# 1) SepalLength, 2) SepalWidth, 3) PetalLength, 4) PetalWidth
julia> labels = convert(Array, iris[:, 5])     # species

### MultivariateStats
julia> Pkg.add("MultivariateStats")
julia> using MultivariateStats

### GLM and GLMNet
julia> Pkg.add("GLM")
julia> Pkg.add("GLMNet")
julia> using GLM, GLMNet, DataFrames

#calculating a Generalized Linear Model:
julia> data = DataFrame();
julia> data[:x] = iris[1:50, :SepalWidth]
50-element DataArrays.DataArray{Float64,1}:
 3.5
 3.0
 3.2
 3.1
 3.6
 3.9
 3.4

julia> data[:y] = iris[1:50, :SepalLength]
50-element DataArrays.DataArray{Float64,1}:
 5.1
 4.9
 4.7
 4.6
 5.0
 5.4
 4.6
 5.0

# LINEAR REGRESSION ---------------------------------------------------------------------------------------------------
# Recall the linreg solution from section 4.3: SepalLength = 2.64466 + 0.690854 * SepalWidth
# Simple linear model:
julia> model = lm(y ~ x, data)
DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.DensePredQR{Float64}},Flo
at64}:

Coefficients:
             Estimate Std.Error t value Pr(>|t|)
(Intercept)     2.639  0.310014 8.51251   <1e-10
x             0.69049 0.0898989 7.68074    <1e-9

# The model's formula is: SepalLength = 2.639 + 0.69049 * SepalWidth

# lm is a shortcut for:
julia> model = fit(LinearModel, y ~ x, data)
DataFrames.DataFrameRegressionModel{GLM.LinearModel{GLM.DensePredQR{Float64}},Flo
at64}:

## GLM general format:
julia> model = glm(formula, data, family, link)

# Family    Canonical Link Function (= optional)
# ------    -----------------------
# Normal    IdentityLink
# Binomial  LogitLink
# Binomial  ProbitLink
# Poisson   LogLink
# Gamma     InverseLink

# example for family = Normal() / link   = IdentityLink()
julia> model = glm(y ~ x, data, Normal(), IdentityLink())
DataFrames.DataFrameRegressionModel{GLM.GeneralizedLinearModel{GLM.GlmResp{Array{
Float64,1},Distributions.Normal,GLM.IdentityLink},GLM.DensePredChol{Float64,Base.
LinAlg.Cholesky{Float64,Array{Float64,2}}}},Float64}:

Coefficients:
             Estimate Std.Error z value Pr(>|z|)
(Intercept)     2.639  0.310014 8.51251   <1e-16
x             0.69049 0.0898989 7.68074   <1e-13
# The model's formula is: SepalLength = 2.639 + 0.0.69049 * SepalWidth

# Estimates for coefficients:
julia> coef(model)
2-element Array{Float64,1}:
 2.639
 0.69049

# Standard errors of coefficients:
julia> stderr(model)
2-element Array{Float64,1}:
 0.310014
 0.0898989

# Covariance of coefficients:
julia> vcov(model)
2x2 Array{Float64,2}:
  0.0961089  -0.0277044
 -0.0277044   0.00808181

# RSS:
julia> deviance(model)
2.731315191455518

julia> GLM.predict(model)
50-element Array{Float64,1}:
 5.05572
 4.71047
 4.84857
 4.77952
 5.12476
 5.33191
 4.98667
 4.98667
 4.64142
 4.77952
 5.19381
 4.98667
 ⋮
 4.71047
 4.98667
 5.05572
 4.22713
 4.84857
 5.05572
 5.26286
 4.71047
 5.26286
 4.84857
 5.19381
 4.91762

# Plot points along with fitted model
julia> using Gadfly
# plot in GLM_model.png
julia> plot(data, layer(x = :x, y = GLM.predict(model), Geom.line(), Theme(default_color = colorant"black")))

### GLMNet is a package for fitting generalised linear models via penalised maximum likelihood.
#
# The shrinkage method is selected via parameter α, where
#
#   α = 0 -> Ridge regression
#   α = 1 -> Lasso regression (the default).
#
# The strength of the penalty is determined by the (output) parameter λ.
#
# More information on GLMNet can be found at http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html.
# Mode information on GLMNet package is at https://github.com/simonster/GLMNet.jl.
#
julia> X = convert(Array, iris[1:50, [:SepalWidth, :PetalLength]])
50x2 Array{Float64,2}:
 3.5  1.4
 3.0  1.4
 3.2  1.3
 3.1  1.5
 3.6  1.4
 3.9  1.7
 3.4  1.4
 3.4  1.5
 2.9  1.4

julia> y = convert(Array, iris[1:50, :SepalLength])
# Fit model for varying values of λ.
#
julia> path = glmnet(X, y)
Least Squares GLMNet Solution Path (58 solutions for 2 predictors in 195 passes):
58x3 DataFrames.DataFrame
| Row | df | pct_dev   | λ          |
|-----|----|-----------|------------|
| 1   | 0  | 0.0       | 0.259109   |
| 2   | 1  | 0.0936139 | 0.236091   |
| 3   | 1  | 0.171334  | 0.215117   |
| 4   | 1  | 0.235858  | 0.196007   |
| 5   | 1  | 0.289428  | 0.178594   |
| 6   | 1  | 0.333902  | 0.162728   |
| 7   | 1  | 0.370825  | 0.148272   |
| 8   | 1  | 0.401479  | 0.1351     |
⋮
| 50  | 2  | 0.570155  | 0.00271447 |
| 51  | 2  | 0.570172  | 0.00247333 |
| 52  | 2  | 0.570187  | 0.0022536  |
| 53  | 2  | 0.570199  | 0.0020534  |
| 54  | 2  | 0.570209  | 0.00187098 |
| 55  | 2  | 0.570217  | 0.00170477 |
| 56  | 2  | 0.570224  | 0.00155332 |
| 57  | 2  | 0.57023   | 0.00141533 |
| 58  | 2  | 0.570234  | 0.00128959 |

# Intercepts for each model.
#
julia> path.a0
58-element Array{Float64,1}:
 5.006
 4.79572
 4.60413
 4.42955
 4.27048
 4.12555
 3.99348
 3.87316
 3.76352
 3.66362
 ⋮
 2.34439
 2.34078
 2.33749
 2.33449
 2.33176
 2.32927
 2.327
 2.32494
 2.32305

# Coefficients for each model.
#
julia> path.betas
2x58 GLMNet.CompressedPredictorMatrix:
 0.0  0.0613412  0.117233  0.16816  0.214562  …  0.663901  0.664214  0.664498
 0.0  0.0        0.0       0.0      0.0          0.275747  0.276429  0.27705

# Now have we models for a variety of λ, how do we know which one is best? Cross-validation!
#
julia> path = glmnetcv(X, y)
Least Squares GLMNet Cross Validation
58 models for 2 predictors in 10 folds
Best λ 0.044 (mean loss 0.065, std 0.018)

# Extract the average loss across cross-validation runs for each model.
#
julia> path.meanloss
58-element Array{Float64,1}:
 0.12247
 0.115895
 0.106013
 0.0974243
 0.0903268
 0.084464
 0.0796237
 0.07563
 0.0728178
 0.0708761
 ⋮
 0.068732
 0.0687994
 0.0688611
 0.0689175
 0.068969
 0.0690161
 0.069059
 0.0690983
 0.0691341

# Which of these gives the smallest loss:
#
julia> indmin(path.meanloss)
20
#
# And what are the corresponding coefficients?
#
julia> path.path.betas[:,indmin(path.meanloss)]
2-element Array{Float64,1}:
 0.567314
 0.0649205

# RIDGE REGRESSION ----------------------------------------------------------------------------------------------------

# α = 1 by default, but we can select α = 0 (or any other value between 0 and 1 for that matter!).
#
julia> path = glmnet(X, y, alpha = 0)
Least Squares GLMNet Solution Path (100 solutions for 2 predictors in 286 passes):
100x3 DataFrames.DataFrame
| Row | df | pct_dev    | λ         |
|-----|----|------------|-----------|
| 1   | 2  | 1.2581e-36 | 259.109   |
| 2   | 2  | 0.00183637 | 236.091   |
| 3   | 2  | 0.00201492 | 215.117   |
| 4   | 2  | 0.0022108  | 196.007   |
| 5   | 2  | 0.00242564 | 178.594   |
| 6   | 2  | 0.0026613  | 162.728   |
| 7   | 2  | 0.00291975 | 148.272   |
| 8   | 2  | 0.0032032  | 135.1     |
⋮
| 92  | 2  | 0.56081    | 0.0545401 |
| 93  | 2  | 0.562226   | 0.0496949 |
| 94  | 2  | 0.563442   | 0.0452801 |
| 95  | 2  | 0.564483   | 0.0412575 |
| 96  | 2  | 0.565374   | 0.0375923 |
| 97  | 2  | 0.566133   | 0.0342527 |
| 98  | 2  | 0.566779   | 0.0312098 |
| 99  | 2  | 0.567327   | 0.0284372 |
| 100 | 2  | 0.567792   | 0.0259109 |

# OTHER OPTIONS -------------------------------------------------------------------------------------------------------

# There are a range of other optional parameters which allow you to fine tune the performance of GLMNet. Amongst these
# are:
#
# - dfmax and pmax (maximum number of predictors);
# - nlambda, lambda_min_ratio and lambda (values of λ);
# - constraints (upper and lower bounds on each predictor);
# - weights (weight vector for samples).

max_stars_repo_path: Section 5/5_4.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 5/5_4.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 5/5_4.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
struct CostVertex
    cost::Float64
    index::Int64
    predecessor::Int64
end
CostVertex(index::Int64) = CostVertex(0.0, index, -1)
Base.isless(a::CostVertex, b::CostVertex) = a.cost < b.cost

mutable struct Pathfinder{T,Q} <: AbstractRoutingAlgorithm
    ruleset::T
    historical_cost::Vector{Float64}
    congestion_cost_factor::Float64
    historical_cost_factor::Float64
    iteration_limit::Int64
    channels_to_route::Vector{ChannelIndex}
    # Scratch-pad to avoid re-allocation.
    discovered::BitVector
    predecessor::Vector{Int64}
    pq::Q

    # Constructor
    function Pathfinder(
        m::Map,
        routing_struct::RoutingStruct;
        iteration_limit = 50,
        channels_to_route = ChannelIndex.(1:length(routing_struct.channels)),
        congestion_cost_factor = 3.0,
        historical_cost_factor = 3.0,
    )
        ruleset = rules(m)
        routing_graph = getgraph(routing_struct)
        num_vertices = nv(routing_graph.graph)

        # Initialize shortest path structures.
        historical_cost = ones(Float64, num_vertices)
        discovered = falses(num_vertices)
        # Initially, all nodes have no predecessors. The initial state of this
        # structure should be cleared to the correct default values.
        predecessor = zeros(Int64, num_vertices)
        pq = BinaryMinHeap{CostVertex}()

        return new{typeof(ruleset),typeof(pq)}(
            ruleset,
            historical_cost,
            congestion_cost_factor,
            historical_cost_factor,
            iteration_limit,
            channels_to_route,
            # scratchpad
            discovered,
            predecessor,
            pq,
        )
    end
end

MapperCore.rules(p::Pathfinder) = p.ruleset

getarchitecture(::Pathfinder{A}) where {A} = A

function channels_to_route(p::Pathfinder, r::RoutingStruct)
    lt = (x, y) -> isless(getchannel(r, x), getchannel(r, y))
    sort!(p.channels_to_route; lt = lt)

    return p.channels_to_route
end

function soft_reset(p::Pathfinder)
    # Zero out the discovered vector.
    p.discovered .= false
    # Empty out the priority queue.
    while !isempty(p.pq)
        pop!(p.pq)
    end
    return nothing
end

function rip_up_routes(p::Pathfinder, rs::RoutingStruct)
    for channel in p.channels_to_route
        clear_route(rs, channel)
    end
end

function linkcost(pf::Pathfinder, rs::RoutingStruct, index::Integer)
    # Get the link info from the routing structure.
    link = getlink(rs, index)
    base_cost = cost(link)
    # Compute the penalty for present congestion
    p = 1 + max(0, pf.congestion_cost_factor * (1 + occupancy(link) - capacity(link)))

    h = pf.historical_cost[index]
    return base_cost * p * h
end

function currentcost(pathfinder, routing_struct, link_idx)
    link = getlink(routing_struct, link_idx)
    base_cost = cost(link)
    p = 1 + max(0, pathfinder.congestion_cost_factor * (occupancy(link) - capacity(link)))

    return base_cost * p
end

function update_historical_congestion(p::Pathfinder, rs::RoutingStruct)
    for i in eachindex(p.historical_cost)
        link = getlink(rs, i)
        overflow = p.historical_cost_factor * (occupancy(link) - capacity(link))
        if overflow > 0
            p.historical_cost[i] += overflow
        end
    end
    return nothing
end

function routecost(
    pathfinder::Pathfinder, routing_struct::RoutingStruct, channel_idx::ChannelIndex
)
    cost = zero(Float64)
    for link_idx in vertices(getroute(routing_struct, channel_idx))
        cost += currentcost(pathfinder, routing_struct, link_idx)
    end
    return cost
end

# Weighted shortest path
function shortest_path(p::Pathfinder, r::RoutingStruct, channel::ChannelIndex)
    A = getarchitecture(p)
    # Reset the runtime structures.
    soft_reset(p)

    # Unpack the certain variables.
    graph = getgraph(r).graph
    pq = p.pq
    discovered = p.discovered
    predecessor = p.predecessor
    # Add all the start nodes to the priority queue.
    start_vectors = start_vertices(r, channel)

    if length(start_vectors) > 1
        throw(ErrorException("""
            Pathfinder cannot run for channels with more than one source component.
            """))
    end
    # Iterate through the list of possible starting points for this channel and
    # add it to the priority queue.
    for node in first(start_vectors)
        push!(pq, CostVertex(node))
    end

    task = getchannel(r, channel)

    # Get the stop nodes for this channel
    stop_vectors = stop_vertices(r, channel)
    stop_mask = trues(length(stop_vectors))

    # Path type that we will store the final paths in.
    paths = SparseDiGraph{Int64}()

    # Keep iterating until all destination nodes have been reached.
    while count(stop_mask) > 0
        # Empty the priority queue
        soft_reset(p)

        # Get the start nodes. If this is the first iteration, get the start
        # nodes from the data structure definition. Otherwise, use the previously
        # used nodes in "paths" as the starting points
        if nv(paths) == 0
            startnodes = first(start_vectors)
            for i in startnodes
                push!(pq, CostVertex(i))
            end
        else
            startnodes = Int64[]
            for j in vertices(paths)
                discovered[j] = true
                # Add all the neighbors of these nodes to the priority queue.
                for u in outneighbors(graph, j)
                    link = getlink(r, u)
                    if !discovered[u] && canuse(rules(p), link, task)
                        new_cost = linkcost(p, r, u)
                        push!(pq, CostVertex(new_cost, u, j))
                    end
                end
                push!(startnodes, j)
            end
        end

        # Get the stop nodes from the undiscovered destinations.
        # Iterate through the vector{vector{int}}, skipping if the destination
        # has already been found.
        stopnodes = Int64[]
        for i in eachindex(stop_vectors)
            stop_mask[i] || continue
            append!(stopnodes, stop_vectors[i])
        end

        # Initialize search variables
        lastnode = 0
        success = false
        while !isempty(pq)
            # Pop a node out of the priority queue
            v = pop!(pq)
            # Check if this node is one of the stop nodes. If so - exit the loop.
            if in(v.index, stopnodes)
                lastnode = v.index
                predecessor[v.index] = v.predecessor
                success = true
                break
            end
            # If the node has not been discovered yet - mark it as discovered and
            # mark its predecessor for a potential backtrace. Add all neighbors
            # of this node to the queue.
            discovered[v.index] && continue
            # Mark node as discovered and all undiscovered neighbors.
            discovered[v.index] = true
            predecessor[v.index] = v.predecessor
            for u in outneighbors(graph, v.index)
                link_type = getlink(r, u)
                if !discovered[u] && canuse(rules(p), link_type, task)
                    # Compute the cost of taking the new vertex and add it
                    # to the queue.
                    new_cost = v.cost + linkcost(p, r, u)
                    push!(pq, CostVertex(new_cost, u, v.index))
                end
            end
        end

        if success == false
            throw(ErrorException("Shortest Path failed on channel index $channel."))
        end
        # Do a back-trace from the last vertex to determine the path that
        # this connection took through the graph.
        add_vertex!(paths, lastnode)

        i = lastnode
        j = predecessor[i]
        # Iterate until we find a start node.
        while j != -1
            add_vertex!(paths, j)
            add_edge!(paths, j, i)
            i, j = j, predecessor[j]
        end

        # Now, find which destination was reached, and mark it as completed
        found_vector = false
        for i in eachindex(stop_vectors)
            stop_mask[i] || continue
            # Check to see if the last discovered node is in this collection
            if in(lastnode, stop_vectors[i])
                stop_mask[i] = false
                found_vector = true
                break
            end
        end

        # Error checking
        if found_vector == false
            throw(ErrorException("""
                Final Node $previous_index is not a valid stop node.
                """))
        end
    end
    # Wrap the path in an edge path and add it to the RoutingStruct.
    setroute(r, paths, channel)
    return nothing
end

function route!(p::Pathfinder, rs::RoutingStruct)
    iterations_per_update = 10
    num_congested_links = 0

    @info "Running Pathfinder Routing Algorithm"
    for i in 1:(p.iteration_limit)
        # Rip up all routes
        for channel in channels_to_route(p, rs)
            clear_route(rs, channel)
        end

        # Run routing
        for channel in channels_to_route(p, rs)
            shortest_path(p, rs, channel)
        end
        update_historical_congestion(p, rs)
        num_congested_links = 0
        # Count the number of congested links.
        for j in p.channels_to_route
            if iscongested(rs, j)
                num_congested_links += 1
            end
        end

        # Report status update.
        if mod(i, iterations_per_update) == 0
            # Report the number of congested links.
            @info """
                On iteration $i of $(p.iteration_limit).
                Number of congested links: $num_congested_links.
                """
            # If congestion is extremely bad, further routings will be super
            # slow. If over 90% of the links are congested, it's very likely
            # that the placement is uunroutable. Just abort early.
            if num_congested_links > 0.9 * length(p.channels_to_route)
                @warn "Warning. Extreme congestion detected - aborting early."
                break
            end
        end
        num_congested_links == 0 && break
    end

    # Perform cleanup if not congested
    if num_congested_links == 0
        cleanup(p, rs)
    end
    return nothing
end

"""
    cleanup(pathfinder)

Post Pathfinder algorithm cleanup routine. Removes historical congestion and
tries to reroute links. Will only accept a new routing if it decreases the
total number of links. Congestion will not be allowed during this phase, so
it can only improve routing quality.

Cleanup process will continue until no improvement it seen.
"""
function cleanup(pathfinder::Pathfinder, routing_struct::RoutingStruct)
    # Set the historical_cost_factor to zero to essentially remove its impact.
    pathfinder.historical_cost_factor = 0.0
    pathfinder.historical_cost .= 1

    # Set the congestion_cost_factor to a large number so congested routes become
    # very sternly discouraged
    pathfinder.congestion_cost_factor = 1.0E6

    # Reroute each link.
    while true
        @debug "Running Cleanup"
        cost_decreased = false
        for channel_idx in pathfinder.channels_to_route
            # Get initial cost of the channel.
            initial_cost = routecost(pathfinder, routing_struct, channel_idx)

            # Save the old routing
            old_routing = getroute(routing_struct, channel_idx)

            # Ripup and reroute
            clear_route(routing_struct, channel_idx)
            shortest_path(pathfinder, routing_struct, channel_idx)

            # Process the final cost.
            final_cost = routecost(pathfinder, routing_struct, channel_idx)

            # @show initial_cost
            # @show final_cost

            if final_cost < initial_cost
                @debug final_cost - initial_cost
                cost_decreased = true
            end
        end

        # Exit when cost stops decreasing.
        cost_decreased || break
    end
    return nothing
end


text:
export transmit, ThreeRayModel, BasicSeabed

Base.@kwdef mutable struct ThreeRayModel <: PropagationModel
  sources::Vector{AcousticSource} = []
  receivers::Vector{AcousticReceiver} = []
  waterdepth = 25.0
  soundspeed = 1500.0
  density = 1023.0
  seabed = BasicSeabed(1600.0, 1600.0, 0.1)
  surfaceloss = db2amp(1.0)
  bottomloss = db2amp(1.0)
  noise = RedGaussianNoise(σ=3000)
  fs = 96000.0
end

struct BasicSeabed
  soundspeed = 1600.0
  density = 1600.0
  absorption = 0.1
end

function transmit(model::ThreeRayModel, x::AbstractVector, spos, rpos)
  rpos2 = copy(rpos)
  rpos2[end] = -rpos2[end]
  rpos3 = copy(rpos)
  rpos3[end] = -2*model.waterdepth - rpos3[end]
  d1 = norm(rpos - spos)
  d2 = norm(rpos2 - spos)
  d3 = norm(rpos3 - spos)
  t1 = round(Int, d1/model.soundspeed*model.fs)
  t2 = round(Int, d2/model.soundspeed*model.fs)
  t3 = round(Int, d3/model.soundspeed*model.fs)
  breflect = reflectioncoef(
    acos((rpos[end]+spos[end])/d2),
    model.seabed.density, model.seabed.soundspeed, model.seabed.absorption,
    model.density, model.soundspeed)
  y = zeros(length(x))
  y[t1:end] .+= x[1:end-t1+1]/d1
  y[t2:end] .+= x[1:end-t2+1]/d2 / model.surfaceloss
  y[t3:end] .+= x[1:end-t3+1]/d3 / model.bottomloss * breflect
  return y
end


text:
# ---
# title: 1035. Uncrossed Lines
# id: problem1035
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/uncrossed-lines/description/>
# hidden: true
# ---
#
# We write the integers of `A` and `B` (in the order they are given) on two
# separate horizontal lines.
#
# Now, we may draw _connecting lines_ : a straight line connecting two numbers
# `A[i]` and `B[j]` such that:
#
#   * `A[i] == B[j]`;
#   * The line we draw does not intersect any other connecting (non-horizontal) line.
#
# Note that a connecting lines cannot intersect even at the endpoints: each
# number can only belong to one connecting line.
#
# Return the maximum number of connecting lines we can draw in this way.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/04/26/142.png)
#
#
#
#     Input: A = [1,4,2], B = [1,2,4]
#     Output: 2
#     Explanation: We can draw 2 uncrossed lines as in the diagram.
#     We cannot draw 3 uncrossed lines, because the line from A[1]=4 to B[2]=4 will intersect the line from A[2]=2 to B[1]=2.
#
#
# **Example 2:**
#
#
#
#     Input: A = [2,5,1,2,5], B = [10,5,2,1,5,2]
#     Output: 3
#
#
# **Example 3:**
#
#
#
#     Input: A = [1,3,7,1,7,5], B = [1,9,2,5,1]
#     Output: 2
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 500`
#   2. `1 <= B.length <= 500`
#   3. `1 <= A[i], B[i] <= 2000`
#
#
## @lc code=start
using LeetCode

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


text:
using Plots

function michalewicz(x; m=10)
    return -sum(sin(v)*sin(i*v^2/π)^(2m) for
        (i,v) in enumerate(x));
end

# plot

heatmap(0:0.1:4, 0:0.1:4, (x, y)->michalewicz([x, y]))


text:
#=
This file defines internal functions for compatibility across
different Julia versions.
=#

using Compat
using Compat: copyto!, axes, argmax
import Compat.String
using Compat.LinearAlgebra
import Compat.LinearAlgebra: norm, checksquare, LAPACKException,
                             SingularException, ×
import Compat.InteractiveUtils.subtypes

using SparseArrays, Printf
@inline _At_mul_B(A, B) = transpose(A) * B
@inline _A_mul_B!(C, A, B) = mul!(C, A, B)
expmat = exp
export eye


text:
# ## Error bar

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Random
Random.seed!(2018)

f(x) = 2 * x + 1
x = 0:0.1:2
n = length(x)
y = f.(x) + randn(n)
plot(x, y,
    xerr=0.1 * rand(n),
    yerr=rand(n))


text:
function createStructures()
    options = alpheusOptions()
    input = alpheusInput()
    preprocessor = alpheusPreprocessor()
    output = alpheusOutput()

    data = alpheusData(options,input,preprocessor,output)

    return data
end

function fillInput!(data::alpheusData)
    input = data.input

    input.J = J
    input.B = B
    input.θ = θ
    input.M = M
    input.Gt = Gt
    input.L = L
    input.n0 = n0
    input.nt = nt
    input.Q = Q
    input.ηu = ηu
    input.ηd = ηd
    input.f = f
    input.Hu = Hu
    input.Hd = Hd
    input.Vmin = Vmin
    input.Vmax = Vmax
    input.N_turbine = N_turbine
    input.Pmax = Pmax

    nothing
end

function fillPreprocessor!(data::alpheusData)
    preprocessor = data.preprocessor
    input = data.input

    preprocessor.I = tan(input.θ);
    preprocessor.Δx = input.L/(input.J-1);
    preprocessor.x = LinRange(0,input.L,input.J);
    preprocessor.E = (input.L.-preprocessor.x)*preprocessor.I;
    preprocessor.K = input.Pmax / input.Q / g /density

    nothing
end

function fillOptions!(data::alpheusData)
    options = data.options

    options.upstream = upstream
    options.downstream = downstream
    options.enumerate = enumerate
    options.solver = solver
    options.parameters = parameters

    nothing
end

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus


text:
# NTRU protocol implementation in Julia
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016
#
# Reference: https://github.com/NTRUOpenSourceProject/ntru-crypto

importall Base
using Krypto

abstract NTRUKey

type NTRUPubKey <: NTRUKey
end

type NTRUPrivKey <: NTRUKey
end


text:
using NewtonsMethod, LinearAlgebra
using Test

f(x)=-x^2; fp(x)=-2x
g(x)=2(x-4)^3; gp(x) = 6(x-4)^2
h(x)=x^3-2x^2-11x+12; hp(x)=3x^2-4x-11
n(x)=x^2+42.0; np(x)=2x
x₀=1.0; xb₀=one(BigFloat); x1=2.35287527; x2=2.35284172
tol=1E-6 # use manual tolerance for known-function tests

for x in [x₀, xb₀]
    @testset "Fns with known roots" begin
        @test norm(newtonroot(f, fp, x)[1]) < tol
        @test norm(newtonroot(g, gp, x)[1] - 4) < tol
        #known basins of attraction for h
        @test norm(newtonroot(h, hp, x1)[1] - 4) < tol
        @test norm(newtonroot(h, hp, x2)[1] + 3) < tol
        #nonconvergence
        @test newtonroot(n, np, x)[1]==nothing
    end

    @testset "Autodiff known fns" begin
        @test norm(newtonroot(f, x)[1]) < tol
        @test norm(newtonroot(g, x)[1] - 4.0) < tol
        @test norm(newtonroot(h, x1)[1] - 4.0) < tol
        @test norm(newtonroot(h, x2)[1] + 3.0) < tol
        @test newtonroot(n, x)[1]==nothing # redundant?
        @test newtonroot(n, x)==newtonroot(n, np, x) # test of f′=Diff(f)
    end
end

@testset "Maxiter, iter, tol" begin
    # expected nonconvergence behavior for complicated function and small maxiter
    @test newtonroot(h, hp, x1, maxiter=5)[1]==nothing
    @test newtonroot(h, x1, maxiter=5)[1]==nothing
    # OK to set tol as exactly 0 when NM yields exact root
    @test newtonroot(h, hp, x1, tol=0.0)[1] == 4.0
    @test newtonroot(h, x1, tol=0.0)[1] == 4.0
    # NM cannot get exact root, so tol=0 returns nothing when x has arbitrary precision
    @test newtonroot(f, fp, one(BigFloat), tol=0.0)[1] == nothing
    @test newtonroot(f, one(BigFloat), tol=0.0)[1] == nothing
    # when tol is nonzero, tol' > tol yields less precise estimate
    @test norm(newtonroot(f, fp, x₀, tol=1E-4)[1] - 0) > norm(newtonroot(f, fp, x₀, tol=1E-5)[1] - 0)
    # tests for iter. nonideal logic, i think.
    @test newtonroot(h, hp, x1, maxiter=0)[3]==1
    @test newtonroot(h, hp, x1, tol=Inf)[3]==1
end


text:
function compute_value_function(inputs)
"""
function for computing the value function
"""

    # Pull in time-invariant model parameters
    states, mparams, exog_states, opts, space, Binv = parameterization(inputs)

    # Copy initial level for later use
    mparams[:initial_smax] = deepcopy(mparams[:smax])
    mparams[:initial_smin] = deepcopy(mparams[:smin])

    # Load terminal value function
    terminal = BSON.load("terminal/terminal_value.bson")
    space_terminal = terminal["space_terminal"]

    # If not loading a previous run start from time-invariant domain
    if ~opts[:load]

        # Initialize domain bounds
        smax = repeat(mparams[:smax]', mparams[:sim_length_adapt], 1)
        smin = repeat(mparams[:smin]', mparams[:sim_length_adapt], 1)

        # Backwards induce the non-adapted problem
        space, c_store, controls = backwards_induce(mparams, exog_states, opts, space, states, Binv, space_terminal)

    end

    # Time-adaptive grid
    if opts[:adaptive_grid]

        if opts[:load]

            smax = inputs[:smax]
            smin = inputs[:smin]

            # Initialize coefficient storage vectors
            c_store = Vector{Vector}(undef, mparams[:time_horizon] + 1)

        end

        # Iterate and adapt the grid
        space, c_store, states, mparams, opts, smax, smin =
            iterate_adaptive_grid(mparams, exog_states, opts, space, states, Binv, space_terminal, smax, smin, c_store)

    end

    if ~opts[:adapted]
        smax = repeat(mparams[:smax]', mparams[:time_horizon], 1)
        smin = repeat(mparams[:smin]', mparams[:time_horizon], 1)
    end

    return space, c_store, states, mparams, opts, smax, smin

end

function backwards_induce(mparams, exog_states, opts, space, states, Binv, space_terminal)
"""
Loops backwards in time and performs VFI
"""

    # Initialize value array
    v = Array{Float64}(undef, space[:num_nodes], 1)

    # Initialize guess for control
    controls = Vector{Vector}(undef, space[:num_nodes])
    for i in eachindex(controls)
        controls[i] = [.1, 1e-4]
    end

    # Initialize coefficient storage vectors
    c_store = Vector{Vector}(undef, mparams[:time_horizon] + 1)

    # Initialize state_variables to pass to the solver
    state_variables = Array{Array}(undef, space[:num_nodes], 1)

    # Initialize current period exogenous variables
    exog_states_t = Dict{Symbol,Any}()

    # Backwards induce the problem
    for t = mparams[:time_horizon]:-1:0

        # Set current year
        mparams[:year] = t

        # Adjust collocation grid and basis functions if increasing number of grid points
        if any(t .== mparams[:increase_level])

            mparams[:approximation_level] = mparams[:approximation_level] .+ mparams[:increase_dims]
            space_new = fundefn(:cheb, mparams[:approximation_level], mparams[:smin], mparams[:smax])
            states = Array{Float64}(undef, prod(mparams[:approximation_level]), length(mparams[:approximation_level]))
            states, bounds = funnode(space_new)
            space_new[:num_nodes] = size(states, 1)
            space[:num_nodes] = size(states, 1)
            space_new[:c_cont] = zeros(space[:num_nodes])
            Binv = funbase(space_new, states)\Matrix{Float64}(I, space_new[:num_nodes], space_new[:num_nodes])
            state_variables = Array{stateVariables}(undef, space[:num_nodes], 1)
            controls = Vector{Vector}(undef, space[:num_nodes])
            [controls[ii] = [.1, .1] for ii in eachindex(controls)]

        end

        # Pull out time-varying variables
        exog_states_t = reduce_dict(exog_states, mparams[:year] + 1)

        # Put in array
        state_variables = [stateVariables(vec(states[node,:]), controls[node], exog_states_t,
            mparams, space, opts, space_terminal, node) for node in eachindex(state_variables)]

        # Iterate on the Bellman
        controls, v = @time iterate_nodes(controls, mparams, opts, state_variables)

        # Store recovered coefficients, controls, and value
        space[:c_cont] = Binv*v

        # If penalty parameter is too small expected continuation value will be Inf
        # and the fitting process will return NaNs for coefficients
        if any(isnan.(space[:c_cont]))
            error("Penalty parameter too small --> breakdown.")
        end

        # Store coefficient vector
        c_store[t+1] = space[:c_cont]

        # Update space if increased approximation level
        if any(mparams[:year] .== mparams[:increase_level])
            space_new[:c_cont] = Binv*v
            space = copy(space_new)
        end

    end

    return space, c_store, controls
end

function backwards_induce(mparams, exog_states, opts, space, states, Binv, space_terminal, smax, smin, adapts)
"""
Loops backwards in time and performs VFI on an adapted grid
"""

    if opts[:learn_e]
        mparams[:approximation_level] = mparams[:initial_approximation_level]
    else
        mparams[:approximation_level] = mparams[:initial_approximation_level][[1 2 5]]
    end

    mparams[:smax] = smax[mparams[:time_horizon] + 1, :]
    mparams[:smin] = smin[mparams[:time_horizon] + 1, :]
    space = fundefn(:cheb, mparams[:approximation_level], mparams[:smin], mparams[:smax])
    space_new = copy(space)
    states = Array{Float64}(undef, prod(mparams[:approximation_level]), length(mparams[:approximation_level]))
    states, bounds = funnode(space)
    space[:num_nodes] = size(states, 1)
    space[:c_cont] = zeros(space[:num_nodes])
    Binv = funbase(space, states)\Matrix{Float64}(I, space[:num_nodes], space[:num_nodes])

    # Initialize value array
    v = zeros(space[:num_nodes], 1)

    # Initialize guess for control
    controls = Vector{Vector}(undef, space[:num_nodes])
    [controls[ii] = [.1, 1e-4] for ii in eachindex(controls)]

    # Coefficient storage vectors
    c_store = Vector{Vector}(undef, mparams[:time_horizon] + 1)

    println("Reapproximating value function on iteration $(mparams[:current_adapt]) with initial approximation level $(mparams[:initial_approximation_level]) and jumps up at $(mparams[:increase_level_adapt]).")

    for t = mparams[:time_horizon]:-1:0

        mparams[:year] = t

        # Check if there was a previous/current increase in approx level
        prev_adapts = mparams[:current_adapt] .- mparams[:increase_adapt_iter]

        # If there was a previous/current adapt
        if any(prev_adapts .>= 0)

            # Get the last approx increase scheme
            last_adapt = mparams[:current_adapt] - minimum(prev_adapts[prev_adapts .>= 0])

            # If at the right year, increase approx level
            if any(adapts .>= mparams[:increase_adapt_iter]) &&
                any(t .== mparams[:increase_level_adapt][findall(mparams[:increase_adapt_iter] .<= adapts)[1]])

                # Current adaptive iteration
                current_iter = findall(mparams[:increase_adapt_iter] .== last_adapt)[1]

                opts[:learn_e] ? idx = collect(1:5) : idx = [1, 2, 5]

                # Increase approximation level given rules for current adaptive iteration and time t
                mparams[:approximation_level] = (mparams[:approximation_level] .+
                mparams[:increase_dims_adapt][current_iter][t .== mparams[:increase_level_adapt][current_iter], idx])

                println("Increased approximation level to  $(mparams[:approximation_level]).")

                # Update approximation space
                space_new = fundefn(:cheb, mparams[:approximation_level], mparams[:smin], mparams[:smax])
                states = Array{Float64}(undef, prod(mparams[:approximation_level]), length(mparams[:approximation_level]))
                states, bounds = funnode(space_new)
                space_new[:num_nodes] = size(states, 1)
                space[:num_nodes] = size(states, 1)
                space_new[:c_cont] = zeros(space[:num_nodes])
                Binv = funbase(space_new, states)\Matrix{Float64}(I, space_new[:num_nodes], space_new[:num_nodes])
                state_variables = Array{stateVariables}(undef, space[:num_nodes], 1)
                controls = Vector{Vector}(undef, space[:num_nodes])
                [controls[ii] = [.1, 1e-4] for ii in eachindex(controls)]

            end

        end

        # Initialize state_variables to pass to the solver
        state_variables = Array{Array}(undef, space[:num_nodes], 1)

        # Initialize current period exogenous variables
        exog_states_t = Dict{Symbol,Any}()

        # Pull out time-varying variables
        exog_states_t = reduce_dict(exog_states, mparams[:year] + 1)

        # Put in array
        state_variables =
            [stateVariables(vec(states[node, :]), controls[node], exog_states_t,
            mparams, space, opts, space_terminal, node) for node in eachindex(state_variables)]

        # Iterate on the Bellman
        controls, v = @time iterate_nodes(controls, mparams, opts, state_variables)

        if t > 0

            if t < mparams[:time_horizon]
                space = copy(space_new)
            end

            mparams[:smax] = smax[t+1, :]
            mparams[:smin] = smin[t+1, :]
            space_new = fundefn(:cheb, mparams[:approximation_level], smin[t,:], smax[t,:])
            space_new[:num_nodes] = size(states, 1)
            states, bounds = funnode(space_new)

        end

        # Store recovered coefficients, controls, and value
        space[:c_cont] = Binv*v
        c_store[t+1] = space[:c_cont]

        if any(t .== mparams[:increase_level_adapt])
            space_new[:c_cont] = Binv*v
            space = copy(space_new)
        end

        # If penalty parameter is too small expected continuation value will be Inf
        # and the fitting process will return NaNs for coefficients
        if any(isnan.(space[:c_cont]))
            error("Penalty parameter too small.")
        end

    end
    return space, c_store, controls
end

function iterate_nodes(controls, mparams, opts, state_variables)
"""
Iterates over the collocation nodes and maximizes the right hand side of the Bellman
"""

    # Initialize control and value vector
    controls = similar(controls)
    v = Array{Float64}(undef, length(controls))

    # Maximization step
    results = pmap(maximize_bellman, state_variables)

    # Pull out maximization results
    for node in eachindex(controls)

        # Consumption and abatement cost
        cons = results[node][1][1]
        cost = results[node][1][2]

        # Output vectors for controls and value
        controls[node] = [cons, cost]
        v[node] = results[node][2]

    end

    # Print current information
    println("----------------------------------------------------")
    println("Year $(mparams[:year]) - Adaptive grid iteration $(mparams[:current_adapt]).")
    if opts[:learn_e]
        println("Domain upper bound: $(mparams[:smax]).")
        println("Domain lower bound: $(mparams[:smin]).")
    end

    return controls, v

end

function iterate_adaptive_grid(mparams, exog_states, opts, space, states, Binv, space_terminal, smax, smin, c_store)
"""
Iterates to adapt the grid
"""

    # Loop over adaptive runs
    for adapts = 1:mparams[:num_adapts]

        # Save current adaptive iteration
        mparams[:current_adapt] = adapts

        # If not loading a previous run, simulate and re-fit the approximation bounds
        if ~opts[:load]

            smax_old = deepcopy(smax)
            smin_old = deepcopy(smin)

            mparams[:sim_length_temp] = copy(mparams[:sim_length])
            mparams[:sim_length] = copy(mparams[:sim_length_adapt])
            mparams[:num_runs] = mparams[:num_runs_adapt]
            opts[:sim] = false
            space_old = deepcopy(space)
            model_parameters_old = deepcopy(mparams)
            println("Adapting grid to simulated path - iteration $(mparams[:current_adapt]).")

            # Get simulated trajectories
            @time consumption, abatement, capital, co2,
                prior_loc, prior_scale, damage, loc_quad, dam_quad =
                simulate_paths_adapt(space, c_store, mparams, opts, smax, smin)

            mparams[:sim_length] = copy(mparams[:sim_length_temp])

            # Start building new approximation bounds
            if opts[:learn_e]
                smax, smin = @time time_adaptive_grid(mparams, capital, co2, prior_loc, prior_scale, damage, loc_quad, dam_quad, states)
            else
                smax, smin = @time time_adaptive_grid(mparams, capital, co2, damage, dam_quad, states)
            end

            # Do spline smoothing
            smax_new = Array{Float64}(undef, mparams[:sim_length_adapt], 5)
            smin_new = Array{Float64}(undef, mparams[:sim_length_adapt], 5)

            for i = [1, 2, 3, 4, 5]

                if i == 1
                    spline_type_max = 1
                    s_type_max = 0.02
                    s_type_min = 2
                    spline_type_min = 3

                elseif i == 5 || i == 4

                    spline_type_max = 2
                    s_type_max = 0.4
                    s_type_min = 3
                    spline_type_min = 3

                else

                    spline_type_max = 2
                    s_type_max = 0.4
                    s_type_min = 0.4
                    spline_type_min = 2

                end

                spline_max = Spline1D(1:mparams[:sim_length_adapt],
                    smax[1:mparams[:sim_length_adapt], i];
                    w = ones(mparams[:sim_length_adapt]), k = spline_type_max, bc = "nearest", s = s_type_max)

                global smax_new[:, i] = evaluate(spline_max, 1:mparams[:sim_length_adapt])

                spline_min = Spline1D(1:mparams[:sim_length_adapt],
                    smin[1:mparams[:sim_length_adapt], i];
                    w = ones(mparams[:sim_length_adapt]), k = spline_type_min, bc = "nearest", s = s_type_min)

                global smin_new[:, i] = evaluate(spline_min, 1:mparams[:sim_length_adapt])

            end

            smax_new = smax_new[:, [1, 2, 3, 4]]
            smin_new = smin_new[:, [1, 2, 3, 4, 5]]

            # Weighted average of old bounds and new bounds
            for tt = 1:mparams[:sim_length_adapt]
                for up_dim = 1:4
                    smax[tt, up_dim] =
                        mparams[:new_bounds_weight_up][up_dim].*smax_new[tt, up_dim] .+
                        (1 - mparams[:new_bounds_weight_up][up_dim]).*smax_old[tt, up_dim]
                end
                for down_dim = 1:5
                    smin[tt,down_dim] =
                        mparams[:new_bounds_weight_down][down_dim].*smin_new[tt,down_dim] .+
                        (1 - mparams[:new_bounds_weight_down][down_dim]).*smin_old[tt, down_dim]
                end
            end

            smin[:,4] = max.(smin[:,4], 0.)

            # Make sure if there are zeros for scale minimum,
            # It doesnt go positive later because of waviness in splines
            for ii = 2:mparams[:sim_length_adapt]
                if smin[ii-1, 4] == 0. && smin[ii, 4] > 0.
                    smin[ii, 4] = 0.
                end
            end

            smax[:,4] = min.(smax[:,4], mparams[:varstart])
            smin[:,4] .= 0.
            smax[:,2] = min.(smax[:,2], mparams[:initial_smax][2])
            smax[mparams[:sim_length_adapt]+1:end, 1] .= smax[mparams[:sim_length_adapt], 1]
            smax[mparams[:sim_length_adapt]+1:end, 4] .= smax[mparams[:sim_length_adapt], 4]
            smin[mparams[:sim_length_adapt]+1:end, 2] .= smin[mparams[:sim_length_adapt], 2]
            smin[mparams[:sim_length_adapt]+1:end, 1] .= mparams[:initial_smin][1]
            smin[mparams[:sim_length_adapt] + 1:end, 2] .= smin[mparams[:sim_length_adapt], 2]
            smax[1,:] = [mparams[:kstart], mparams[:Mstart], mparams[:expstart], mparams[:varstart], mparams[:lossstart]].*1.001
            smin[1,:] = [mparams[:kstart], mparams[:Mstart], mparams[:expstart], mparams[:varstart], mparams[:lossstart]].*0.999
            smax[:,5] .= 1.0

        end

        # Set up for approximation
        opts[:load] = false
        opts[:adapted] = true

        # Backwards induce the adapted problem
        space, c_store, controls =
            backwards_induce(mparams, exog_states, opts, space, states, Binv, space_terminal, smax, smin, adapts)

    end

    return space, c_store, states, mparams, opts, smax, smin

end

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej


text:
using ABCDBeamTrace
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "beamtrace" begin
    f = 100e-3
    L = 1000e-3
    expander_2x = [ThinLens(f), FreeSpace(3f), ThinLens(2f)]
    system = [expander_2x; FreeSpace(L); reverse(expander_2x)]
    beam = Beam(1000e-9, 1.0e-3)
    @test beamtrace(system, beam) isa Vector{Beam}
end

@testset "discretize" begin
    f = 100e-3
    expander_2x = [ThinLens(f), FreeSpace(3f), ThinLens(2f)]
    @test length(discretize(expander_2x,10)) == 12
end

@testset "waistradiusfunc" begin
    w0 = 1e-3
    f = 100e-3
    expander_2x = [ThinLens(f), FreeSpace(3f), ThinLens(2f)]
    beam = Beam(1000e-9, w0)
    @test waistradiusfunc(expander_2x, beam) isa Function
    @test waistradiusfunc(expander_2x, beam)(0) == w0
    @test waistradiusfunc(expander_2x, beam)(3f) ≈ 2w0 rtol=0.01
end


text:
using SFML

f(x) = sin(x)
const thickness = 2
const max_point_distance = 10000^2
const movespeed = 0.5
pixels_per_unit = 20

settings = ContextSettings()
settings.antialiasing_level = 4
window = RenderWindow(VideoMode(800, 600), "Plot", settings, window_defaultstyle)
set_framerate_limit(window, 60)
event = Event()

view = View(Vector2f(0, 0), Vector2f(800, -600))

circle = CircleShape()
set_position(circle, Vector2f(-400, 0))
set_radius(circle, thickness/2)
set_origin(circle, Vector2f(thickness/2, thickness/2))
set_fillcolor(circle, SFML.red)

global points = CircleShape[]

xaxis = RectangleShape()
set_fillcolor(xaxis, Color(128, 128, 128))
set_size(xaxis, Vector2f(get_size(view).x, 2))
set_origin(xaxis, Vector2f(get_size(view).x / 2, 1))

yaxis = copy(xaxis)
rotate(yaxis, 90)

function redraw()
	set_position(xaxis, Vector2f(get_center(view).x, 0))
	set_position(yaxis, Vector2f(0, get_center(view).y))

	empty!(points)
	center = get_center(view)
	i = center.x - get_size(view).x/2
	while i <= center.x + get_size(view).x
		point = copy(circle)
		try
			set_position(point, Vector2f(i, pixels_per_unit * f(i / pixels_per_unit)))
		catch DomainError
			i += 1
			continue
		end
		push!(points, point)
		i += 1
	end
end

redraw()

old_mousepos = Vector2i(0, 0)
while isopen(window)
	# mousepos = pixel2coords(window, get_mousepos(window), view)
	# println("$(mousepos.x / pixels_per_unit) $(mousepos.y / pixels_per_unit)")
	sleep(0)
	while pollevent(window, event)
		if get_type(event) == EventType.CLOSED
			close(window)
		end
		if get_type(event) == EventType.MOUSE_MOVED
			if is_mouse_pressed(0)
				mousepos = get_mousemove(event)
				move(view, Vector2f(old_mousepos.x - mousepos.x, mousepos.y - old_mousepos.y))
				old_mousepos = mousepos
				redraw()
			end
		end
		if get_type(event) == EventType.MOUSE_BUTTON_PRESSED
			old_mousepos = get_mousebutton(event)
		end
		if get_type(event) == EventType.MOUSE_WHEEL_MOVED
			oldcenter = Vector2f(get_center(view).x / pixels_per_unit, get_center(view).y / pixels_per_unit)
			delta = get_mousewheel(event).delta
			pixels_per_unit += delta / 10
			if pixels_per_unit < 0
				pixels_per_unit = 0.01
			end
			# println(pixels_per_unit)
			set_center(view, Vector2f(oldcenter.x*pixels_per_unit, oldcenter.y*pixels_per_unit))
			redraw()
		end
	end

	set_view(window, view)

	clear(window, SFML.white)
	draw(window, xaxis)
	draw(window, yaxis)
	for i = 1:length(points)
		if i > 1
			global l = Line(get_position(points[i]), get_position(points[i - 1]), thickness)
			# if distance_squared(l.p1, l.p2)*pixels_per_unit < max_point_distance
			# set_fillcolor(l, SFML.red)
			draw(window, l)
			# end
		end
		draw(window, points[i])
	end

	display(window)
end


text:
module MyExample

greet() = print("Hello World!")
include("extra_file.jl")
include("new_file.jl")

export my_f, new_f

end

max_stars_repo_name: Mikejmnez/MyExample

max_issues_repo_name: Mikejmnez/MyExample

max_forks_repo_name: Mikejmnez/MyExample


text:
include("bernoulli.jl")
include("beta_uniform.jl")
include("beta.jl")
include("binom.jl")
include("categorical.jl")
include("exponential.jl")
include("gamma.jl")
include("geometric.jl")
include("inv_gamma.jl")
include("laplace.jl")
include("mvnormal.jl")
include("neg_binom.jl")
include("normal.jl")
include("piecewise_uniform.jl")
include("poisson.jl")
include("uniform_continuous.jl")
include("uniform_discrete.jl")


text:
export FiniteStrainTruss
mutable struct FiniteStrainTruss
    eledim::Int64
    mat  # constitutive law
    elnodes::Array{Int64}   # the node indices in this finite element
    props::Dict{String, Any}
    coords::Array{Float64}
    weights::Array{Float64}
    hs::Array{Array{Float64}}
    stress::Array{Float64}
    l0::Float64                  # length in the undeformed shape
    A0::Float64                  # cross section area of the undeformed shape
    rot_mat::Array{Float64}      # 4 by 4 rotation matrix
end

function FiniteStrainTruss(coords::Array{Float64}, elnodes::Array{Int64}, props::Dict{String, Any}, ngp::Int64=2)
    eledim = 1
    weights, hs = get1DElemShapeData( coords, ngp )
    nGauss = length(weights)
    name = props["name"]
    if name=="Elasticity1D"
        mat = [Elasticity1D(props) for i = 1:nGauss]
    elseif name=="Plasticity1D"
        mat = [Plasticity1D(props) for i = 1:nGauss]
    elseif name=="Viscoplasticity1D"
        mat = [Viscoplasticity1D(props) for i = 1:nGauss]
    elseif name=="NeuralNetwork1D"
        mat = [NeuralNetwork1D(props) for i = 1:nGauss]
    elseif name=="PathDependent1D"
        mat = [PathDependent1D(props) for i = 1:nGauss]
    else
        error("Not implemented yet: $name")
    end
    strain = Array{Float64}(undef, length(weights))

    dx, dy = coords[2,1] - coords[1,1], coords[2,2] - coords[1,2]
    l0 = sqrt(dx^2 + dy^2)
    cosα, sinα = dx/l0, dy/l0
    rot_mat = [cosα   0    sinα   0;
               0     cosα    0    sinα;
               -sinα  0    cosα   0;
               0     -sinα   0    cosα]

    FiniteStrainTruss(eledim, mat, elnodes, props, coords, weights, hs, strain, l0, props["A0"], rot_mat)
end

function getStiffAndForce(self::FiniteStrainTruss, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)

    rot_mat = self.rot_mat
    #rotate displacement to the local coordinate
    lu = rot_mat * state
    Dlu = rot_mat * Dstate

    A0, l0 = self.A0, self.l0

    E = (lu[2]-lu[1])/l0 + 0.5*((lu[2]-lu[1])/l0)^2 + 0.5*((lu[4]-lu[3])/l0)^2
    DE = (Dlu[2]-Dlu[1])/l0 + 0.5*((Dlu[2]-Dlu[1])/l0)^2 + 0.5*((Dlu[4]-Dlu[3])/l0)^2
    # compute  ∂E∂u.T, 4 by 1 array
    ∂E∂u =  [-1/l0+(lu[1]-lu[2])/l0;  1/l0+(lu[2]-lu[1])/l0; (lu[3]-lu[4])/l0; (lu[4]-lu[3])/l0]
    ∂∂E∂∂u =  [1/l0  -1/l0     0        0;
              -1/l0   1/l0     0        0;
               0       0      1/l0    -1/l0;
               0       0     -1/l0    1/l0]

    for k = 1:length(self.weights)

        # #@show E, DE
        S, dS_dE = getStress(self.mat[k], E, DE, Δt)

        self.stress[k] = S

        fint += A0 * self.weights[k] * S * ∂E∂u  # 4x1

        stiff += A0 * self.weights[k] * (dS_dE * ∂E∂u * ∂E∂u' + S * ∂∂E∂∂u) # 4x4
    end

    #rotate back to global coordinate

    fint = rot_mat' * fint
    stiff = rot_mat' * stiff * rot_mat

    return fint, stiff
end

function getInternalForce(self::FiniteStrainTruss, state::Array{Float64}, Dstate::Array{Float64}, Δt::Float64)
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)

    rot_mat = self.rot_mat
    #rotate displacement to the local coordinate
    lu = rot_mat * state
    Dlu = rot_mat * Dstate

    A0, l0 = self.A0, self.l0

    E = (lu[2]-lu[1])/l0 + 0.5*((lu[2]-lu[1])/l0)^2 + 0.5*((lu[4]-lu[3])/l0)^2
    DE = (Dlu[2]-Dlu[1])/l0 + 0.5*((Dlu[2]-Dlu[1])/l0)^2 + 0.5*((Dlu[4]-Dlu[3])/l0)^2
    # compute  ∂E∂u.T, 4 by 1 array
    ∂E∂u =  [-1/l0+(lu[1]-lu[2])/l0;  1/l0+(lu[2]-lu[1])/l0; (lu[3]-lu[4])/l0; (lu[4]-lu[3])/l0]

    for k = 1:length(self.weights)

        # #@show E, DE
        S, _ = getStress(self.mat[k], E, DE, Δt)

        self.stress[k] = S

        fint += A0 * self.weights[k] * S * ∂E∂u  # 4x1

    end

    #rotate back to global coordinate

    fint = rot_mat' * fint

    return fint
end

function getMassMatrix(self::FiniteStrainTruss)
    ndofs = dofCount(self)
    nnodes = length(self.elnodes)
    mass = zeros(ndofs,ndofs)

    A0 = self.A0
    for k = 1:length(self.weights)
        rho = self.mat[k].ρ

        mass += [self.hs[k]*self.hs[k]' zeros(nnodes, nnodes)
                 zeros(nnodes, nnodes)  self.hs[k]*self.hs[k]']  * A0 * rho * self.weights[k]
    end
    lumped = sum(mass, dims=2)
    #@show mass
    mass, lumped
end

function getStrain(self::FiniteStrainTruss, state::Array{Float64})
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    fint = zeros(Float64, ndofs)
    stiff = zeros(Float64, ndofs,ndofs)

    rot_mat = self.rot_mat
    #rotate displacement to the local coordinate
    lu = rot_mat * state

    A0, l0 = self.A0, self.l0

    nGauss = length(self.weights)
    E = zeros(nGauss, 1)
    w∂E∂u = zeros(nGauss, 4, 1)

    for k = 1:nGauss
        # compute  ∂E∂u.T, 8 by 3 array
        E[k,1] = (lu[2]-lu[1])/l0 + 0.5*((lu[2]-lu[1])/l0)^2 + 0.5*((lu[4]-lu[3])/l0)^2
        w∂E∂u[k,:,:] = rot_mat' * [-1/l0+(lu[1]-lu[2])/l0;  1/l0+(lu[2]-lu[1])/l0; (lu[3]-lu[4])/l0; (lu[4]-lu[3])/l0] * A0 * self.weights[k]
    end
    # @info "state ", state , "lu", lu , "w∂E∂u[k,:,:] ", w∂E∂u[1,:,:]

    return E, w∂E∂u
end

function getNodes(self::FiniteStrainTruss)
    return self.elnodes
end

function dofCount(self::FiniteStrainTruss)
    return 2length(self.elnodes)
end

function commitHistory(self::FiniteStrainTruss)
    for m in self.mat
        commitHistory(m)
    end
end

function  getForceAndDforceDstress(self::FiniteStrainTruss, state::Array{Float64},
    stress::Array{Float64,2})
    ndofs = dofCount(self);
    nnodes = length(self.elnodes)
    nStrain = 1
    nGauss = length(self.weights)
    fint = zeros(Float64, ndofs)
    dfint_dstress = zeros(Float64,  ndofs, nGauss * nStrain)

    rot_mat = self.rot_mat
    #rotate displacement to the local coordinate
    lu = rot_mat * state

    A0, l0 = self.A0, self.l0

    E = (lu[2]-lu[1])/l0 + 0.5*((lu[2]-lu[1])/l0)^2 + 0.5*((lu[4]-lu[3])/l0)^2
    # compute  ∂E∂u.T, 4 by 1 array
    ∂E∂u =  [-1/l0+(lu[1]-lu[2])/l0;  1/l0+(lu[2]-lu[1])/l0; (lu[3]-lu[4])/l0; (lu[4]-lu[3])/l0]

    for k = 1:length(self.weights)

        # #@show E, DE
        S = stress[k, :]

        self.stress[k] = S[1]

        fint += A0 * self.weights[k]  * ∂E∂u  * S[1]# 4x1

        dfint_dstress[:, (k-1)*nStrain+1:k*nStrain] = A0 * ∂E∂u * self.weights[k]

    end

    #rotate back to global coordinate

    fint = rot_mat' * fint

    dfint_dstress = rot_mat' * dfint_dstress

    return fint , dfint_dstress
end


text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Simulate the collapse of a self-gravitating scalar field
# in Vaidya spacetime
# FIXME: Compute initial data
#--------------------------------------------------------------------

using NLsolve, Printf

function psi(u::T, v::T)::T where {T<:Number}
    t = v + u
    r = v - u
    if r == 0.0
        return cos(k*t)
    else
        return cos(k*t) * (sin(k*r) / (k*r))
    end
end

function computeUboundary(PS::ProductSpace{S1, S2})::NTuple{3, Field{S2}} where {S1, S2}
    a0 = extractUboundary(Field(PS, (u,v)->2), :incoming)
    η0 = extractUboundary(Field(PS, (u,v)->(v-u)), :incoming)
    ϕ0 = extractUboundary(Field(PS, psi), :incoming)
    ηs = computeη((a0, η0, ϕ0))
    # @show norm(lineconstraints((a0, ηs, ϕ0)))
    # @assert norm(lineconstraints((a0, ηs, ϕ0))) < 1e-9
    return (a0, ηs, ϕ0)
end

function computeVboundary(PS::ProductSpace{S1, S2})::NTuple{3, Field{S1}} where {S1, S2}
    a0 = extractVboundary(Field(PS, (u,v)->2), :incoming)
    η0 = extractVboundary(Field(PS, (u,v)->(v-u)), :incoming)
    ϕ0 = extractVboundary(Field(PS, psi), :incoming)
    ηs = computeη((a0, η0, ϕ0))
    # @show norm(lineconstraints((a0, ηs, ϕ0)))
    # @assert norm(lineconstraints((a0, ηs, ϕ0))) < 1e-9
    return (a0, ηs, ϕ0)
end

function initialguess(PS::ProductSpace{S1, S2})::NTuple{3, Field{ProductSpace{S1, S2}}} where {S1, S2}
    a0 = Field(PS, (u,v)->2)
    η0 = Field(PS, (u,v)->v-u)
    ϕ0 = Field(PS, psi)
    return (a0, η0, ϕ0)
end

function computePatch(PS::ProductSpace{S1, S2}, ubnd::NTuple{3, Field{S2}},
                                                vbnd::NTuple{3, Field{S1}})::NTuple{3, Field{ProductSpace{S1, S2}}} where {S1, S2}
    B = incomingboundary(PS)
    I = identity(PS)
    DU, DV = derivative(PS)
    ∂U = combineUVboundary(ubnd, vbnd, :incoming)

    function constr(U::NTuple{3, Field{S}})::NTuple{3, Field{S}} where {S}
        (a, η, ϕ) = U
        C1 = DU*(DU*η) - (2/a)*(DU*a)*(DU*η) + (4pi*η)*(DU*ϕ)^2
        C2 = DV*(DV*η) - (2/a)*(DV*a)*(DV*η) + (4pi*η)*(DV*ϕ)^2
        C  = sqrt(C1*C1 + C2*C2)
        return (C, C, C)
    end

    function offaxis(U::NTuple{3, Field{S}})::NTuple{3, Field{S}} where {S}
        (a, η, ϕ) = U
        F1 = DU*(DV*a) - (1/a)*(DU*a)*(DV*a) + (a/η)*(DU*(DV*η)) + (4pi*a)*(DU*ϕ)*(DV*ϕ)
        F2 = DU*(DV*η) + (1/η)*(DU*η)*(DV*η) + (1/4)*(1/η)*(a^2)
        F3 = DU*(DV*ϕ) + (1/η)*(DU*η)*(DV*ϕ) + (1/η)*(DV*η)*(DU*ϕ)
        return (F1, F2, F3)
    end

    function onaxis(U::NTuple{3, Field{S}})::NTuple{3, Field{S}} where {S}
        (a, η, ϕ) = U
        A1 = (DU*η) * (DV*η) + (1/4) * a^2
        A2 = DU*(DV*η)
        A3 = DU*(DV*ϕ)
        return (A1, A2, A3)
    end

    function boundary(U::NTuple{3, Field{S}})::NTuple{3, Field{S}} where {S}
        return U - ∂U
    end

    function residual(U::NTuple{3, Field{S}})::NTuple{3, Field{S}} where {S}
        F = offaxis(U) + λ * constr(U)
        enforceregularity!(F, onaxis(U))
        enforcebc!(F, boundary(U))
        return F
    end

    function f!(fvec::Array{T,1}, x::Array{T,1}) where {T}
        fvec[:] = reshapeFromTuple(residual(reshapeToTuple(PS, 3, x)))
    end

    s = nlsolve(f!, reshapeFromTuple(initialguess(PS)); method=:trust_region,
                autodiff=:forward, show_trace=debug >= 3, ftol=1e-9, iterations=40)
    k = reshapeToTuple(PS, 3, s.zero)

    if debug >= 2
        println(" [Before solve] norm(constraints) = ", norm(constraints(k)))
        println("                Converged?        = ", converged(s))
        println(" [After solve]  norm(constraints) = ", norm(constraints(k)))
    end

    return k
end

function pconv(range::StepRange, nh::Int=1)
    println("Testing p convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        np = n_[index]
        params = Parameters((np, np), (nh,nh), urange, vrange, 3)
        q = distribute(params, computePatch, computeUboundary, computeVboundary)
        l_[index] = rmse(extract(map(constraints, q), 1))
        @printf("  n = %i, rmse = %e\n", n_[index], l_[index])
    end
    return (n_, l_)
end

function hconv(range::StepRange, np::Int=4)
    println("Testing h convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        nh = n_[index] = 2^n_[index]
        params = Parameters((np, np), (nh, nh), urange, vrange, 3)
        l_[index] = rmse(extract(map(constraints, distribute(params, computePatch, computeUboundary, computeVboundary)), 1))
        l0 = index[1] > 1 ? l_[index[1] - 1] : 1
        @printf("  n = %3i, rmse[2h] / rmse[h] = %e\n", n_[index], l0 / l_[index])
    end
    return (n_, l_)
end

function constraint(U::NTuple{3, Field{S}})::NTuple{1, Field{S}} where {S}
    (a, η, ϕ) = U
    PS = first(U).space
    DU, DV = derivative(PS)
    C1 = DU*(DU*η) - (2/a)*(DU*a)*(DU*η) + (4pi*η)*(DU*ϕ)^2
    C2 = DV*(DV*η) - (2/a)*(DV*a)*(DV*η) + (4pi*η)*(DV*ϕ)^2
    C  = sqrt(C1*C1 + C2*C2)
    return (C,)
end

#-----------------------------------------
# Setup simulation grid parameters,
#-----------------------------------------
npoints  = (4, 4)
npatches = (4, 4)
urange   = (2.0, 4.0)
vrange   = (5.0, 7.0)
params   = Parameters(npoints, npatches, urange, vrange, 3)
λ        = 0.0
debug    = 1
k        = 0.1

S = distribute(params, computePatch, computeUboundary, computeVboundary, debug)
save("output/data/scalarwave_SG/scalarwave_SG_psi", extract(S, 3))
save("output/data/scalarwave_SG/scalarwave_SG_constraints", extract(map(constraint, S), 1))

# np, ep = pconv(2:2:14, 2)
nh, eh = hconv(1:1:5, 4)
# save("output/data/scalarwave_SG/scalarwave_SG_k0p7_off_axis_pconv", np, ep)
# save("output/data/scalarwave_SG/scalarwave_SG_k0p7_off_axis_hconv", nh, eh)


text: include("priorsampling.jl")


text:
# ---
# title: 344. Reverse String
# id: problem344
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Two Pointers, String
# link: <https://leetcode.com/problems/reverse-string/description/>
# hidden: true
# ---
#
# Write a function that reverses a string. The input string is given as an array
# of characters `char[]`.
#
# Do not allocate extra space for another array, you must do this by **modifying
# the input array  [in-place](https://en.wikipedia.org/wiki/In-
# place_algorithm)** with O(1) extra memory.
#
# You may assume all the characters consist of [printable ascii
# characters](https://en.wikipedia.org/wiki/ASCII#Printable_characters).
#
#
#
# **Example 1:**
#
#
#
#     Input: ["h","e","l","l","o"]
#     Output: ["o","l","l","e","h"]
#
#
# **Example 2:**
#
#
#
#     Input: ["H","a","n","n","a","h"]
#     Output: ["h","a","n","n","a","H"]
#
#
#
## @lc code=start
using LeetCode

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


text:
function multiply(a::Number, b::Number)
  return a * b
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
systemid=parse(Int8,ARGS[1])
tend=parse(Float64,ARGS[2])
atol=parse(Float64,ARGS[3])
rtol=parse(Float64,ARGS[4])

using LinearAlgebra: BLAS
BLAS.set_num_threads(1) # Recording CPU time means we should
                        # stick to one thread

# Load data
load_data = :(begin
    using SQLite
    db = SQLite.DB(ENV["KURABENCH_DB"])

    # First read the system info
    results = DBInterface.execute(db,string(
            "select ",
                "name, nodes, edges, coupling_constant " ,
            "from ",
                "systems ",
            "where " ,
                "id=$systemid;"))
    name = ""
    for row in results
        global name = row.name
        global nodes = row.nodes
        global edges = row.edges
        global coupling_constant = row.coupling_constant
    end
    # Now read connectivity
    # Allocate
    sources = Array{Int32}(undef,edges)
    dests = Array{Int32}(undef,edges)
    data = ones(Int32,edges)

    results = DBInterface.execute(db,string(
            "select ",
                "id,source,dest ",
             "from " ,
                "edges " ,
             "where " ,
                "system=$systemid;"))

    for (i,row) in enumerate(results)
        sources[i] = row.source
        dests[i]   = row.dest
    end
    # Store the node omegas
    omegas = Array{Float64}(undef,nodes)
    x0     = Array{Float64}(undef,nodes)
    results = DBInterface.execute(db,string(
            "select ",
                "id,omega,initialCondition ",
            "from " ,
                "nodes " ,
            "where " ,
                "system=$systemid;"))
    for (i,row) in enumerate(results)
        omegas[i] = row.omega
        x0[i] = row.initialCondition
    end
    my_matrix = sparse(sources,dests,data)
end)

prepare_ode = :(begin
    g = SimpleGraph(my_matrix')
    integrator = DP5
    ### Network dynamics vertex and edge functions

	@inline Base.@propagate_inbounds function kuramoto_vertex!(dv, v, edges, p, t)
		dv[1] = p # parameter for vertices are eigenfrequencies omega
		@inbounds for e in edges
	            dv[1] += e[1]
	    end
	    nothing
	end

    @inline Base.@propagate_inbounds function kuramoto_edge!(e,v_s,v_d,p,t)
        e[1] = p * sin(v_s[1] - v_d[1])
        nothing
    end

    ### Constructing the network dynamics

    odevertex = ODEVertex(f! = kuramoto_vertex!, dim = 1)
    staticedge = StaticEdge(f! = kuramoto_edge!, dim = 1, coupling = :antisymmetric)

    # generating random values between -0.5 and 0.5 for the parameter value ω_0 of the vertices
    e_pars = [coupling_constant for i in 1:g.ne]
    parameters = (omegas, e_pars)

    # setting up the  network_dynamics
    kuramoto_network! = network_dynamics(odevertex, staticedge, g)

    ### Simulation
    # constructing random initial conditions for nodes (variable θ)

    prob = ODEProblem(kuramoto_network!, x0, (0.,tend), parameters, abstol=atol, reltol=rtol)

end)

# Precompile eval
eval(:nothing)

using Pkg
Pkg.activate(@__DIR__)
println("Enviroment Activated, starting program...")
using CPUTime

# Times
tstartup = @CPUelapsed begin
	using NetworkDynamics: ODEVertex, StaticEdge, network_dynamics
	using LightGraphs: SimpleGraph
	using OrdinaryDiffEq: ODEProblem, solve, DP5
    using SparseArrays
end
times = 0.:1.:tend
println("time to startup is $tstartup")
tload    = @CPUelapsed(eval(load_data))
println("time to load is $tload")
tjit     = @CPUelapsed(eval(prepare_ode)) +
           @CPUelapsed(solve(prob, integrator(), saveat=times))

println("time to jit is $tjit")
tintegrate = @CPUelapsed(sol = solve(prob, integrator(), saveat=times))

println("time to solve is $tintegrate")

ttotal = tstartup + tload + tjit + tintegrate
# Now write the config to the experiments table and the states to the states table
#

## Config
tidy = cmd -> chomp(String(read(cmd)))

githash = tidy(`git rev-parse HEAD`)
hostname = tidy(`hostname`)
datetime = tidy(`date`)
gitchanges = String(read(`git status -s -uno`))
argsstring = join(ARGS," ")
intstring = string(integrator)
DBInterface.execute(db,string(
    "insert into experiments ",
    "(runtime,datetime,githash,gitchanges,commandline,host,",
    "system,solver,atol,rtol,tstart,tend,",
    "tstartup,tload,tjit,tintegration,ttotal)",
    "values ",
    "(\"networkdynamics.jl\",",
    "\"$datetime\",\"$githash\",\"$gitchanges\",\"$argsstring\",\"$hostname\",",
    "$systemid,\"$intstring\",$atol,$rtol,0,$tend,",
    "$tstartup,$tload,$tjit,$tintegrate,$ttotal);"
                              )
                   )

results = DBInterface.execute(db,"select last_insert_rowid()")
experimentid = 0
for row in results
    # TODO: How to get a single element from julia iterator?
    global experimentid = row[1]
end

for (i,t) in enumerate(sol.t)
    for (k,state) in enumerate(sol[i])
        state = mod(state,2*pi)
        DBInterface.execute(db,string(
            "insert into states",
            "(experimentid,time,idx,node,state)",
            "values",
            "($experimentid,$t,$k,$k,$state);"
                                     )
                           )
    end
end

# TODO: Get the reference trajectory for this request (if it exists)
# and update the expermients table (see the python_complete file for guide)
results =DBInterface.execute(db,string(
            "select experimentid from experiments where",
            " system=\"$systemid\"",
            " and tend=$tend",
            " and atol=1e-12",
            " and rtol=1e-10",
            " and solver is \"radau\";"
                                     )
                   )
refid = -1
for row in results
    global refid=row.experimentid
    break
end
errors = zeros(length(sol[end]))
for (idx,state) in enumerate(last(sol))
    local results = DBInterface.execute(db,string(
            "select state from states",
            " where experimentid=$refid",
            " and time=$tend",
            " and idx=$idx;"
                                           )
                                 )
    for row in results
        global errors[idx] = abs(row.state-mod(state,2*pi))
    end
end

avgErr = sum(errors) / length(errors)

DBInterface.execute(db,string(
        "update experiments",
        " set err_v_ref=$avgErr,",
        " refeid=$refid",
        " where experimentid=$experimentid;"
                             )
                   )

max_stars_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks

max_issues_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks

max_forks_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks


text:
using Vec
using Base.Test

include("test_vecE2.jl")
include("test_vecE3.jl")
include("test_vecSE2.jl")

include("test_coordinate_transforms.jl")
include("test_geomE2.jl")


text:
using Flux
using Tracker
using Statistics
using Distributed
using OptimalPlasticityRatios

const cQuality  = 1 ./[0.1,1, 5]
const sysfluct_ratios =  0.1:0.3:2.5
const disturbance_sizes = [0.01, 0.02, 0.03]
const num_reps = 8

const numInputs = 12
const numOutputs = 10

const batch_size = 1000
const training_size = 1000
# const test_size = 500

const num_hidden_units = 20
const training_cycles = 8000

function make_student(num_hidden_units::Integer)
    return Chain(
  Dense(numInputs, numOutputs))
end

students = [make_student(num_hidden_units) for qual in cQuality, dis in disturbance_sizes, rat in sysfluct_ratios, rep in 1:num_reps]
teachers = [make_student(num_hidden_units) for rep in 1:num_reps]

teacher_weights = [Tracker.data.(params(teach)) for teach in teachers]

for i in 1:num_reps
  [Flux.loadparams!(student, (1. + 0.001*randn(1)[1]).*teacher_weights[i]) for student in students[:,:,:,i]]
end

# [Flux.loadparams!(student, (1. + 0.001*randn(1)[1]).*teacher_weights) for student in students]

tr_data = randn(Float32,numInputs,training_size)

ranger = [1:m for m in size(students)]
li = [OptimalPlasticityRatios.studentTeacherInfo(students[i,j,k,m], teachers[m], Flux.mse)
for (i,j,k,m) in Base.product(ranger...)]

# li = OptimalPlasticityRatios.studentTeacherInfo.(students,(teacher,),(Flux.mse,))
tr = [OptimalPlasticityRatios.trainingDataOnly(tr_data, batch_size, training_cycles) for el in li]
### rat is ||Δc||: ||Δϵ|| ratio. dis is  ||Δϵ||

n_weights = OptimalPlasticityRatios.get_num_weights(li[1])
disturbance_sizes *= sqrt(n_weights)
println(disturbance_sizes)
upd = [OptimalPlasticityRatios.fixedNormFluctuatingNoisyGD(qual,rat*dis,dis) for qual in cQuality, dis in disturbance_sizes, rat in sysfluct_ratios, rep in 1:num_reps]

to_record = NamedTuple{(:loss,)}((training_cycles-100:training_cycles,))

recs = OptimalPlasticityRatios.make_records_template.(li, (to_record,))
# @time out = OptimalPlasticityRatios.training_loop!.(upd, li, (tr,), recs )

@time  outs = map(OptimalPlasticityRatios.training_loop!, upd,li,tr,recs)

stores = first.(outs)
recs = last.(outs)
d = Dict{Symbol, Any}(:recs => recs)

using BSON

bson("linear_call.bson", d)

include("dataParsing.jl")
recLossOnly = [el[:loss] for el in recs]
dataStore = extract_endmean.(recLossOnly)

ps = make_all_plots(sysfluct_ratios, recLossOnly)

# l = @layout [ a b c ;
#               d e f ;
#               g h i ]

# qq = plot(ps..., layout=l)

# sp = deepcopy(ps)
# for i = 1:3
#   for j = 1:3
#     sp[i,j] = ps[j,i]
#   end
# end

pplot(r::recorder) = plot(r.index, r.value)


text:
using SecondQuantization
using Test

using SecondQuantization: flatten_term
using Latexify

@testset "Normal Ordering" begin
    @boson a b c
    @parameter ω
    N = a' * a

    @test isequal(normal_order(a * a'), normal_order(a' * a + 1))
    @test isequal(normal_order(a * a', rewriter = threaded_normal_order_simplifier(100)), normal_order(a' * a + 1))

    to_test = to_test = b / ω + b * c + a' * a * b / ω + a' * a * b * c
    @test isequal(normal_order(a * (ω * c+1) * b * a' / ω), to_test)
    @test isequal(normal_order(a * (ω * c+1) * b * a' / ω, rewriter = threaded_normal_order_simplifier(100)), to_test)

    @test isequal(normal_order((a + a')^2), flatten_term(+, 1 + a^2 + flatten_term(*, 2 * a' * a) + a'^2))
    @test isequal(normal_order((a + a')^2, rewriter = threaded_normal_order_simplifier(100)), flatten_term(+, 1 + a^2 + flatten_term(*, 2 * a' * a) + a'^2))

    @test isequal(normal_order(a / a), 1)

    @test isequal(normal_order(a * N * N * N), a + 7 * a' * a^2 + 6 * a'^2 * a^3 + a'^3 * a^4)
end

@testset "Latexify" begin
    @boson a b
    @parameter ω δ

    @test latexify(Expr(:latexifymerge, "\\hat{a}")) == latexify(a)
    @test latexify(Expr(:latexifymerge, "\\hat{a}^\\dagger")) == latexify(a')
    @test latexify(Expr(:latexifymerge, "\\hat{a}^\\dagger\\hat{b}")) == latexify(a' * b)
    @test latexify(Expr(:latexifymerge, "\\hat{a}^\\dagger\\left(1+\\hat{b}\\right)")) == latexify(a' * (1 + b))
    @test latexify(Expr(:latexifymerge, "\\frac{1}{\\omega}\\hat{a}+\\frac{1}{\\omega}\\hat{a}^\\dagger")) == latexify((a + a') / ω)
    @test latexify(Expr(:latexifymerge, "-1+\\hat{a}^\\dagger\\hat{a}")) == latexify(a' * a - 1)
    @test latexify(Expr(:latexifymerge, "\\sin \\left(\\hat{a}+\\hat{a}^\\dagger\\right)")) == latexify(sin(a + a'))
    @test latexify(Expr(:latexifymerge, "\\cos \\left(\\hat{a}+\\hat{a}^\\dagger\\right)")) == latexify(cos(a + a'))
    @test latexify(Expr(:latexifymerge, "e^{\\hat{a}+\\hat{a}^\\dagger}")) == latexify(exp(a + a'))
    @test latexify(Expr(:latexifymerge, "\\hat{a}-\\hat{a}^\\dagger")) == latexify(normal_order(a - a'))
    @test latexify(Expr(:latexifymerge, "{\\hat{a}}^{2}")) == latexify(a^2)
    @test latexify(Expr(:latexifymerge, "\\left(\\hat{a}+\\hat{a}^\\dagger\\right)^{2}")) == latexify((a + a')^2)
    @test latexify(Expr(:latexifymerge, "\\delta\\omega+\\hat{a}")) == latexify(normal_order(a + δ * ω))
    @test latexify(Expr(:latexifymerge, "\\frac{\\delta}{\\omega}+\\hat{a}")) == latexify(a + δ / ω)
    @test latexify(Expr(:latexifymerge, "{\\delta}^{\\omega}+\\hat{a}")) == latexify(a + δ ^ ω)
    @test latexify(Expr(:latexifymerge, "\\left(\\delta+\\omega\\right)^{\\omega}+\\hat{a}")) == latexify(a + (δ+ω) ^ ω)
end

@testset "Commutator" begin
    @boson a b
    @parameter ω δ

    @test isequal(commutator(a, b), 0)
    @test isequal(commutator(a, a), 0)
    @test isequal(commutator(a, a'), 1)
    @test isequal(commutator(a', a), -1)
    @test isequal(commutator(a, a' * ω), ω)
    @test isequal(commutator(a * ω^2 / (ω + δ), a'), ω^2 / (ω + δ))
    @test isequal(commutator(a * ω^2, a' * δ^2), ω^2 * δ^2)
end


text:
using Shell
using Test

@testset "ls" begin
    if Sys.iswindows()
        @test Shell.run("dir") === nothing
        @test split(Shell.run("cd", capture=true))[end] == pwd()
        Shell.setchomp(false)
        @test endswith(Shell.run("echo love", capture=true), "\n")

        Shell.setshell("powershell")
        Shell.setchomp(true)
        @test Shell.run("ls") === nothing
        @test Shell.run("(Get-Item -Path . -Verbose).FullName", capture=true) == pwd()
        Shell.setchomp(false)
        @test endswith(Shell.run("pwd", capture=true), "\n")
        @test !endswith(Shell.run("(Get-Item -Path . -Verbose).FullName", capture=true, chomp=true), "\n")
    else
        Shell.setshell("bash")
        @test Shell.run("ls") === nothing
        @test Shell.run("pwd", capture=true) == pwd()
        Shell.setchomp(false)
        @test endswith(Shell.run("pwd", capture=true), "\n")
        @test !endswith(Shell.run("pwd", capture=true, chomp=true), "\n")

        file = "temp file 0"
        @test esc`$file` == "'temp file 0'"
        files = ["temp file 1", "temp file 2"]
        filelist = esc`$files.txt`
        @test filelist == "'temp file 1.txt' 'temp file 2.txt'"
        @test Shell.run("touch $filelist") === nothing
        @test Shell.run("rm $filelist") === nothing
    end
end


text:
module Diagrams

using CombinatorialSpaces.ExteriorCalculus
using Catlab
using Catlab.CategoricalAlgebra
using Catlab.Theories
using Catlab.Programs.DiagrammaticPrograms: parse_diagram, parse_gat_expr
using Unicode

export @decapode

macro decapode(cat, body)
  :(parse_exp_diagram($(esc(cat)), $(Meta.quot(body)), free=true))
end

function parse_exp_diagram(cat, body; kw...)
  objs = filter(exp -> !(exp isa LineNumberNode) && (exp.head == :(::)), body.args)
  obj_map = Dict{Symbol, Expr}()
  for obj in objs
    names = obj.args[1]
    type = obj.args[2]
    # Expand if object is a tuple or single object
    if names isa Symbol
      obj_map[names] = type
    else
      for name in names.args
        obj_map[name] = type
      end
    end
  end
  new_body = quote end
  for exp in body.args
    if !(exp isa LineNumberNode) && exp.head == :call && exp.args[1] == :(==)
      exprs = []
      if exp.args[2] isa Symbol
        res, _ = expand_expr!(exprs, exp.args[3], cat, obj_map; make_var = false)
        push!(exprs, :($(exp.args[2]) == $(res)))
      elseif exp.args[3] isa Symbol
        res, _ = expand_expr!(exprs, exp.args[2], cat, obj_map; make_var = false)
        push!(exprs, :($(exp.args[3]) == $(res)))
      else
        res1, _ = expand_expr!(exprs, exp.args[2], cat, obj_map; make_var = true)
        res2, _ = expand_expr!(exprs, exp.args[3], cat, obj_map; make_var = false)
        push!(exprs, :($(res1) == $(res2)))
      end
      append!(new_body.args, exprs)
    else
      push!(new_body.args, exp)
    end
  end
  parse_diagram(cat, new_body; kw...)
end

i2sub = Dict(
  '0'=>'₀', '1'=>'₁', '2'=>'₂', '3'=>'₃', '4'=>'₄', '5'=>'₅',
  '6'=>'₆', '7'=>'₇', '8'=>'₈', '9'=>'₉', '-'=>'₋'
)
i2sup = Dict(
  '0'=>'⁰', '1'=>'¹', '2'=>'²', '3'=>'³', '4'=>'⁴', '5'=>'⁵',
  '6'=>'⁶', '7'=>'⁷', '8'=>'⁸', '9'=>'⁹', '-'=>'⁻'
)

form_sup(type::Expr) = form_sup(type.args[1])
form_sub(type::Expr) = form_sub(type.args[1])

form_sup(type::Symbol) = startswith("$type", "Dual") ?
                                               i2sup[last("$type")] * '̃' :  i2sup[last("$type")]
form_sub(type::Symbol) = startswith("$type", "Dual") ? i2sub[last("$type")] * '̃' : i2sub[last("$type")]

# Currently every projection needs its own index
function make_proj(ind, types)
  op = "proj" * i2sub["$ind"[1]] * "_" * join(form_sup.(types)) * form_sub(types[ind])
  Symbol(Unicode.normalize(op))
end

function parse_type(type)
  if head(type) == :otimes
    Meta.parse("otimes{$(join(["$(head(t)){$(t.args[1])}" for t in type.args], ","))}")
  else
    Meta.parse("$(head(type)){$(type.args[1])}")
  end
end

get_name(expr) = expr isa Symbol ? expr : get_name(expr.args[1])

function expand_expr!(expr_arr, expr, cat, obj_map; make_var = true)
  if expr isa Symbol
    expr, obj_map[expr]
  elseif length(expr.args) > 2
    if expr.args[1] == :+
      # Only process 2 of these at a time (sum first, evaluate second as sum)
      res1 = expand_expr!(expr_arr, expr.args[2], cat, obj_map; make_var = true)

      res2 = length(expr.args) > 3 ?
      expand_expr!(expr_arr, :(+($(expr.args[3:end]...))), cat, obj_map; make_var = true) :
              expand_expr!(expr_arr, expr.args[3], cat, obj_map; make_var = true)
      type = res1[2]

      sum_name = Symbol(Unicode.normalize("sum" * form_sub(type)))
      new_var = gensym()
      nv_type = :(otimes{$(res1[2]), $(res2[2])})
      push!(expr_arr, :($(new_var)::$(nv_type)))
      push!(expr_arr, :($(res1[1]) == $(make_proj(1, [res1[2], res2[2]]))($(new_var))))
      push!(expr_arr, :($(res2[1]) == $(make_proj(2, [res1[2], res2[2]]))($(new_var))))
      if make_var
        name = gensym()
        push!(expr_arr, :($(name)::$(type)))
        push!(expr_arr, :($(name) == $(sum_name)($(new_var))))
        name, type
      else
        :($(sum_name)($(new_var))), type
      end
    else
      res = map(expr.args[2:end]) do ex
        expand_expr!(expr_arr, ex, cat, obj_map; make_var = true)
      end
      new_var = gensym()
      nv_type = Meta.parse("otimes{$(join(last.(res), ","))}")
      push!(expr_arr, :($(new_var)::$(nv_type)))
      append!(expr_arr, map(enumerate(res)) do (i, obj)
                :($(obj[1]) == $(make_proj(i, last.(res)))($(new_var)))
              end)
      type = parse_type(codom(parse_gat_expr(FinCat(cat), expr.args[1])))
      if make_var
        name = gensym()
        push!(expr_arr, :($(name)::$(type)))
        push!(expr_arr, :($(name) == $(expr.args[1])($(new_var))))
        name, type
      else
        :($(expr.args[1])($(new_var))), type
      end
    end
  else
    type = parse_type(codom(parse_gat_expr(FinCat(cat), expr.args[1])))
    arg_var, arg_type = expand_expr!(expr_arr, expr.args[2], cat, obj_map; make_var = false)
    expr.args[2] = arg_var
    cur_var = expr
    if make_var
      name = gensym()
      push!(expr_arr, :($(name)::$(type)))
      push!(expr_arr, :($(name) == $(expr)))
      cur_var = name
    end
    cur_var, type
  end
end
end


text:
# Tutorial on conventional 2D SLAM example
# This tutorial shows how to use some of the commonly used factor types
# This tutorial follows from the ContinuousScalar example from IncrementalInference
using GraffSDK
using GraffSDK.DataHelpers
using ProgressMeter
using UUIDs
using Caesar

# 1a. Create a Configuration
config = loadGraffConfig()
config.sessionId = "HexDemoSample1_"*replace(string(uuid4())[1:6], "-" => "")
println(getGraffConfig())

# 1b. Check the credentials and the service status
printStatus()
# 1c. Check the session queue length
@info "Session backlog (queue length) = $(getSessionBacklog())"

# 2. Confirm that the robot already exists, create if it doesn't.
println(" - Creating or retrieving robot '$(config.robotId)'...")
robot = nothing
if isRobotExisting()
    println(" -- Robot '$(config.robotId)' already exists, retrieving it...")
    robot = getRobot();
else
    # Create a new one programatically - can also do this via the UI.
    println(" -- Robot '$(config.robotId)' doesn't exist, creating it...")
    newRobot = RobotRequest(config.robotId, "My New Bot", "Description of my neat robot", "Active");
    robot = addRobot(newRobot);
end
println(robot)

# 3. Create or retrieve the session.
# Get sessions, if it already exists, add to it.
println(" - Creating or retrieving data session '$(config.sessionId)' for robot...")
session = nothing
if isSessionExisting()
    println(" -- Session '$(config.sessionId)' already exists for robot '$(config.robotId)', retrieving it...")
    session = getSession()
else
    # Create a new one
    println(" -- Session '$(config.sessionId)' doesn't exist for robot '$(config.robotId)', creating it...")
    newSessionRequest = SessionDetailsRequest(config.sessionId, "A test dataset demonstrating data ingestion for a wheeled vehicle driving in a hexagon.", "Pose2")
    session = addSession(newSessionRequest)
end
println(session)

# 4. Drive around in a hexagon
addVariable(:x0, Pose2)
# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor([:x0], IIF.Prior( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2)) )))

println(" - Adding hexagonal driving pattern to session...")
numVariables = 6
@showprogress for i in 1:numVariables
    psym = Symbol("x$(i-1)")
    nsym = Symbol("x$(i)")
    addVariable(nsym, Pose2)

    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    addFactor([psym;nsym], pp )
end

# 5. Now lets add a couple landmarks
# Ref: https://github.com/dehann/RoME.jl/blob/master/examples/Slam2dExample.jl#L35
response = addVariable(:l0, Point2)
p2br2 = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor([:x0; :l0], p2br2)
addFactor([:x6; :l0], p2br2)

# Landmarks generally require more work once they're created, e.g. creating factors,
# so they are not set to ready by default. Once you've completed all the factor links and want to solve,
# call putReady to tell the solver it can use the new nodes. This is added to the end of the processing queue.
putReady(true)

# 5. Checking solver status, getting data
# Let's wait for all nodes to be processed
@time while getSessionBacklog() > 0
    @info "...Session backlog currently has $(getSessionBacklog()) entries, waiting until complete..."
    sleep(2)
end
@info "Okay processing the last node! ...Yeah so we're going to focus on making that faster in the next release (please go write an issue in GraffSDK to spur us along!)"

# Lets see that everything processed successfully
@info "Session dead queue length = $(getSessionDeadQueueLength())"
if getSessionDeadQueueLength() > 0
    @error "This shouldn't happen, please examine the failed messages below to see what went wrong:"
    deadMsgs = getSessionDeadQueueMessages()
    map(d -> println(d["error"]), deadMsgs)
else
    @info "No messages in the dead queue (errors) on the server, looks good!"
end
# You can ask to reprocess them, or delete them with these commands:
# reprocessDeadQueueMessages()
# deleteDeadQueueMessages()

# Now add some data - we need to make this asynchronous as well. ASAP
imgRequest = DataHelpers.readFileIntoDataRequest(joinpath(dirname(pathof(GraffSDK)), "..", "examples", "pexels-photo-1004665.jpeg"), "TestImage", "Pretty neat public domain image", "image/jpeg");
@showprogress for i in 1:numVariables
    println("  - Adding a simple (largish) image data to the pose x$(i)...")
    # Adding image data
    setData("x$i", imgRequest)
end

# Now get a list of all the variables in the graph
@time nodes = GraffSDK.ls()

# Get all data
getSessionDataEntries()

# Get the estimates to see how far along it is (this just retrieves node.properties["MAP_est"])
# but simplifies getting info
getEstimates()
# You can also give it a list, e.g. just the landmarks
getEstimates(nodes=getLandmarks())

# Other things you can do: Get the full node detail with getNode:
# # By NeoID
# node = getNode( nodes.nodes[1].id)
# # By Graff label
# node = getNode( nodes.nodes[1].label)

# Or checking on solver status:
sessionLatest = getSession()
@info "Graff last solved: $(sessionLatest.lastSolvedTimestamp)"
# Manually request a full session solve if you would like to make sure all is good.
requestSessionSolve()
while sessionLatest.lastSolvedTimestamp == getSession().lastSolvedTimestamp
    @info "Waiting on updated full graph solve..."
    sleep(2)
end
@info "Graff last solved: $(getSession().lastSolvedTimestamp)"