text:
# high-level interface

### Synthesis for LPC variants ###

function synthesis(excitation::AbstractVector,
                   param::SpectralParamState{LinearPredictionCoef},
                   hopsize::Integer)
    f = AllPoleDF(param_order(paramdef(param)))
    synthesis!(f, excitation, to_filtcoef(f, param), hopsize)
end

function synthesis(excitation::AbstractVector,
                   param::SpectralParamState{PartialAutoCorrelation},
                   hopsize::Integer)
    f = AllPoleLatticeDF(param_order(paramdef(param)))
    synthesis!(f, excitation, to_filtcoef(f, param), hopsize)
end

function synthesis(excitation::AbstractVector,
                   param::SpectralParamState{LineSpectralPair},
                   hopsize::Integer)
    f = LSPDF(param_order(paramdef(param)))
    synthesis!(f, excitation, to_filtcoef(f, param), hopsize)
end

### Synthesis for mel-generalized cepstrums

function synthesis{T<:MelGeneralizedCepstrum}(excitation::AbstractVector,
                                              param::SpectralParamState{T},
                                              hopsize::Integer)
    def = paramdef(param)
    ns = try
        # InexactError could happen
        convert(Int, -1/glog_gamma(def))
    catch e
        ns = -1/glog_gamma(def)
        error("-1/γ must be able to convert to Int, but got $(ns)")
    end
    f = MGLSADF(param_order(def), allpass_alpha(def), ns)
    synthesis!(f, excitation, to_filtcoef(f, param), hopsize)
end

function synthesis(excitation::AbstractVector,
                   param::SpectralParamState{MelCepstrum},
                   hopsize::Integer)
    def = paramdef(param)
    f = MLSADF(param_order(def), allpass_alpha(def))
    synthesis!(f, excitation, to_filtcoef(f, param), hopsize)
end

function synthesis(excitation::AbstractVector,
                   param::SpectralParamState{LinearCepstrum},
                   hopsize::Integer)
    f = LMADF(param_order(paramdef(param)))
    synthesis!(f, excitation, to_filtcoef(f, param), hopsize)
end


text:
# manipulation of class labels

## simple classification

# classify

function classify(x::RealVector, ord::Ordering)
    n = length(x)
    v = x[1]
    k::Int = 1
    for i = 2:n
        @inbounds xi = x[i]
        if lt(ord, v, xi)
            v = xi
            k = i
        end
    end
    return k
end

classify(x::RealVector) = classify(x, Forward)

function classify!(r::IntegerVector, x::RealMatrix, ord::Ordering)
    m = size(x, 1)
    n = size(x, 2)
    length(r) == n || throw(DimensionMismatch("Mismatched length of r."))
    for j = 1:n
        @inbounds r[j] = classify(view(x,:,j), ord)
    end
    return r
end

classify!(r::IntegerVector, x::RealMatrix) = classify!(r, x, Forward)

# - this one throws a deprecation
classify(x::RealMatrix, ord::Ordering) = classify!(Array{Int}(undef, size(x,2)), x, ord)
classify(x::RealMatrix) = classify(x, Forward)

# classify with score(s)

function classify_withscore(x::RealVector, ord::Ordering)
    n = length(x)
    v = x[1]
    k::Int = 1
    for i = 2:n
        @inbounds xi = x[i]
        if lt(ord, v, xi)
            v = xi
            k = i
        end
    end
    return (k, v)
end

classify_withscore(x::RealVector) = classify_withscore(x, Forward)

function classify_withscores!(r::IntegerVector, s::RealVector, x::RealMatrix, ord::Ordering)
    m = size(x, 1)
    n = size(x, 2)
    length(r) == n || throw(DimensionMismatch("Mismatched length of r."))
    for j = 1:n
        k, v = classify_withscore(view(x,:,j), ord)
        @inbounds r[j] = k
        @inbounds s[j] = v
    end
    return (r, s)
end

classify_withscores!(r::IntegerVector, s::RealVector, x::RealMatrix) =
    classify_withscores!(r, s, x, Forward)

function classify_withscores(x::RealMatrix{T}, ord::Ordering) where T<:Real
    n = size(x, 2)
    r = Array{Int}(undef, n)
    s = Array{T}(undef, n)
    return classify_withscores!(r, s, x, ord)
end

classify_withscores(x::RealMatrix{T}) where {T<:Real} = classify_withscores(x, Forward)

# classify with threshold

classify(x::RealVector, t::Real, ord::Ordering) =
    ((k, v) = classify_withscore(x, ord); ifelse(lt(ord, v, t), 0, k))

classify(x::RealVector, t::Real) = classify(x, t, Forward)

function classify!(r::IntegerVector, x::RealMatrix, t::Real, ord::Ordering)
    m = size(x, 1)
    n = size(x, 2)
    length(r) == n || throw(DimensionMismatch("Mismatched length of r."))
    for j = 1:n
        @inbounds r[j] = classify(view(x,:,j), t, ord)
    end
    return r
end

classify!(r::IntegerVector, x::RealMatrix, t::Real) = classify!(r, x, t, Forward)

classify(x::RealMatrix, t::Real, ord::Ordering) = classify!(Array{Int}(undef, size(x,2)), x, t, ord)
classify(x::RealMatrix, t::Real) = classify(x, t, Forward)

## label map

struct LabelMap{K}
    vs::Vector{K}
    v2i::Dict{K,Int}

    function LabelMap{K}(vs, v2i) where K
        length(vs) == length(v2i) || throw(DimensionMismatch("lengths of vs and v2i mismatch"))
        new(vs,v2i)
    end
end

LabelMap(vs::Vector{K}, v2i::Dict{K,Int}) where {K}= LabelMap{K}(vs, v2i)

length(lmap::LabelMap) = length(lmap.vs)
keys(lmap::LabelMap) = lmap.vs

function show(io::IO, lmap::LabelMap)
    n = length(lmap)
    println(io, "LabelMap (with $n labels):")
    for (i, v) in enumerate(keys(lmap))
        println(io, "[$i] $v")
    end
end

# build a label map (value -> label) from a sequences of values
function labelmap(xs::AbstractArray{T}) where T
    l = 0
    vs = T[]
    v2i = Dict{T, Int}()
    for x in xs
        if !haskey(v2i, x)
            push!(vs, x)
            v2i[x] = (l += 1)
        end
    end
    return LabelMap(vs, v2i)
end

# use a map to encode discrete values into labels
labelencode(lmap::LabelMap{T}, x) where {T} = lmap.v2i[convert(T, x)]
labelencode(lmap::LabelMap{T}, xs::AbstractArray{T}) where {T} =
    reshape(Int[labelencode(lmap, x) for x in xs], size(xs))

# decode the label to the associated discrete value
labeldecode(lmap::LabelMap{T}, y::Int) where {T} = lmap.vs[y]
labeldecode(lmap::LabelMap{T}, ys::AbstractArray{Int}) where {T} =
    reshape(T[labeldecode(lmap, y) for y in ys], size(ys))

## group labels

function groupindices(k::Int, xs::IntegerVector; warning::Bool=true)
    gs = Array{Vector{Int}}(undef, k)
    for i = 1:k
        gs[i] = Int[]
    end

    warned = false
    n = length(xs)
    for i = 1:n
        @inbounds x = xs[i]
        if x < 1 || x > k
            if warning && !warned
                warn("Some labels are out of bound.")
            end
        end
        push!(gs[x], i)
    end
    return gs
end

function groupindices(lmap::LabelMap{T}, xs::AbstractArray{T}) where T
    k = length(lmap)
    gs = Array{Vector{Int}}(undef, k)
    for i = 1:k
        gs[i] = Int[]
    end
    n = length(xs)
    for i = 1:n
        @inbounds l = labelencode(lmap, xs[i])
        push!(gs[l], i)
    end
    return gs
end

ext: jl


text:
using DataFrames
using Unitful
using PhysicalConstants.CODATA2018

@testset "reco" begin
    data_file = "testdata/evtPar_example.h5"

    _, rmin, rmax, evtdf = read_evtpar([data_file])

    # Single xy calculation
    first_evt = first(evtdf)
    x1, y1 = ATools.calculate_xy(first_evt.xr1, first_evt.yr1, first_evt.r1)
    @test isapprox(x1, first_evt.x1)
    @test isapprox(y1, first_evt.y1)

    # With broadcasting
    xv, yv, zv = ATools.radial_correction(evtdf.xr1, evtdf.yr1, evtdf.zr1, evtdf.r1)
    @test all(isapprox.(zv, evtdf.zr1))
    @test all(isapprox.(xv, evtdf.x1))
    @test all(isapprox.(yv, evtdf.y1))

    # Interaction radius.
    const_val, lin_val, bias = 373.80, -1.22, -0.13
    rad_pred = ATools.predict_interaction_radius(ATools.gpol1([const_val, lin_val]),
                                                 rmin, rmax, bias)
    transform!(evtdf, :zstd1 => rad_pred => :r1x, :zstd2 => rad_pred => :r2x)
    r1_test = const_val .+ first(evtdf.zstd1, 5) .* lin_val .+ bias
    r2_test = const_val .+ first(evtdf.zstd2, 5) .* lin_val .+ bias
    @test all(isapprox.(r1_test, first(evtdf.r1x, 5)))
    @test all(isapprox.(r2_test, first(evtdf.r2x, 5)))

    # Units conversion
    units_evtdf = ATools.set_units(evtdf)
    r1_type     = eltype(units_evtdf.r1)
    t1_type     = eltype(units_evtdf.t1)
    @test      r1_type  <: Unitful.Length
    @test      t1_type  <: Unitful.Time
    @test unit(r1_type) == Unitful.mm
    @test unit(t1_type) == Unitful.ns

    # Time of flight
    flight_time_r = ATools.time_of_flight(units_evtdf, :r1, rmax)
    @test      eltype(flight_time_r)  <: Unitful.Time
    @test unit(eltype(flight_time_r)) == Unitful.ps
    full_width_ps = uconvert(Unitful.ps, (rmax - rmin) * Unitful.mm * 1.6 / SpeedOfLightInVacuum)
    @test all(flight_time_r .<= full_width_ps)

    flight_time_disp = ATools.time_of_flight(units_evtdf, [:xs, :ys, :zs], [:x1, :y1, :z1])
    @test      eltype(flight_time_disp)  <: Unitful.Time
    @test unit(eltype(flight_time_disp)) == Unitful.ps
    expected_first6 = [ 988.117 * Unitful.ps, 1318.657 * Unitful.ps,
                        900.190 * Unitful.ps, 1029.788 * Unitful.ps,
                       1097.411 * Unitful.ps, 1210.753 * Unitful.ps]
    @test all(isapprox.(expected_first6, first(flight_time_disp, 6), atol=0.0005 * Unitful.ps))
end

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 8, n}) where n
    io = openInput(puzzle)
    split(read(io, String), "\n", keepempty=false)
end

function parseInput(input)
    m = match(r"^(nop|acc|jmp) (.\d+)$", input)
    argument = parse(Int, m.captures[2])
    return (operation = m.captures[1], argument )
end

struct InfiniteLoopDetected <: Exception
    location::Int
    accumulator::Int
end

function runInstruction(instruction, state)
    (operation, argument) = instruction

    if (operation == "nop")
        return (loc = state.loc + 1, acc = state.acc)
    elseif (operation == "acc")
        return (loc = state.loc + 1, acc = state.acc + argument)
    elseif (operation == "jmp")
        return (loc = state.loc + argument, acc = state.acc)
    end

    throw(ArgumentError("Got an unkown operation: $operation"))
end

function visitInstruction!(visited, state)
    (loc, acc) = state
    if loc in visited
        throw(InfiniteLoopDetected(loc, acc))
    else
        push!(visited, loc)
    end
end

function runInstructions(instructions)
    acc = 0
    loc = 1
    state = (loc = loc, acc = acc)
    visited = [];
    while (state.loc <= length(instructions))
        visitInstruction!(visited, state)

        state = runInstruction(instructions[state.loc], state)
    end
    return state.acc
end

function solve(::Puzzle{2020, 8, 1}, input)
    instructions = map(parseInput, input)
    try
        runInstructions(instructions)
    catch e
        if (typeof(e) == InfiniteLoopDetected)
            return e.accumulator
        end
        throw(e)
    end
end

function switchOperation((operation, argument))
    if (operation == "nop")
        return (operation = "jmp", argument)
    elseif (operation == "jmp")
        return (operation = "nop", argument)
    end

    throw(ArgumentError("Got an unkown operation: $operation"))
end

function repairBootCode!(instructions, errorLocation)
    instruction = instructions[errorLocation]

    if (instruction[1] ∉ ["jmp", "nop"])
        return false
    end

    instructions[errorLocation] = switchOperation(instruction)
    return true
end

function solve(::Puzzle{2020, 8, 2}, input)
    instructions = map(parseInput, input)
    for errorLocation = 1:length(instructions)
        if (repairBootCode!(instructions, errorLocation) == false)
            continue
        end

        try
            return runInstructions(instructions)
        catch e
            if (typeof(e) == InfiniteLoopDetected)
                repairBootCode!(instructions, errorLocation)
                continue
            end
            throw(e)
        end
    end
    return nothing
end

@testset "Day 8" begin
    TEST = split("""
nop +0
acc +1
jmp +4
acc +3
jmp -3
acc -99
acc +1
jmp -4
acc +6
""", '\n', keepempty=false)

    @test(parseInput("nop +0") == (operation = "nop", argument = 0))
    @test(parseInput("acc +4") == (operation = "acc", argument = 4))
    @test(parseInput("jmp -4") == (operation = "jmp", argument = -4))
    @test(solve(Puzzle(8, 1), TEST) == 5)
    @test(solve(Puzzle(8, 1), input(Puzzle(8, 1))) == 1586)

    @test(solve(Puzzle(8, 2), TEST) == 8)
    @test(solve(Puzzle(8, 2), input(Puzzle(8, 2))) == 703)
end


text:
@recipe function f(mdp::VDPTagMDP, s::TagState)
    ratio --> :equal
    xlim --> (-5, 5)
    ylim --> (-5, 5)
    @series begin
        color := :black
        seriestype := :scatter
        label := "target"
        markersize := 0.1
        [s.target[1]], [s.target[2]]
    end
    @series begin
        color --> :blue
        label --> "agent"
        pts = Plots.partialcircle(0, 2*pi, 100, mdp.tag_radius)
        x, y = Plots.unzip(pts)
        x+s.agent[1], y+s.agent[2]
    end
end

@recipe function f(pomdp::VDPTagPOMDP, h::AbstractPOMDPHistory{TagState})
    mdp = pomdp.mdp
    ratio --> :equal
    xlim --> (-5, 5)
    ylim --> (-5, 5)
    @series begin
        label := "path"
        x = [s.agent[1] for s in state_hist(h)[1:end-1]]
        y = [s.agent[2] for s in state_hist(h)[1:end-1]]
        x, y
    end
    @series begin
        if action_hist(h)[end].look
            color := :blue
        else
            color := :red
        end
        s = state_hist(h)[end-1]
        label := "current agent position"
        pts = Plots.partialcircle(0, 2*pi, 100, mdp.tag_radius)
        x, y = Plots.unzip(pts)
        x+s.agent[1], y+s.agent[2]
    end
    @series begin
        label := "belief"
        belief_hist(h)[end]
    end
    @series begin
        seriestype := :scatter
        label := "current target"
        pos = state_hist(h)[end-1].target
        color --> :orange
        [pos[1]], [pos[2]]
    end
end

@recipe function f(pc::ParticleCollection{TagState})
    seriestype := :scatter
    x = [p.target[1] for p in particles(pc)]
    y = [p.target[2] for p in particles(pc)]
    color --> :black
    markersize --> 0.1
    x, y
end


text:
using Random, Statistics, FreqTables
using LinearAlgebra, Optim, GLM
using DataFrames, CSV, HTTP

function ps2()
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 1
#:::::::::::::::::::::::::::::::::::::::::::::::::::
f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
minusf(x) = x[1]^4+10x[1]^3+2x[1]^2+3x[1]+2
startval = rand(1)   # random starting value
result = optimize(minusf, startval, BFGS())

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function ols(beta, X, y)
    ssr = (y.-X*beta)'*(y.-X*beta)
    return ssr
end

beta_hat_ols = optimize(b -> ols(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(beta_hat_ols.minimizer)

bols = inv(X'*X)*X'*y
df.white = df.race.==1
bols_lm = lm(@formula(married ~ age + white + collgrad), df)
println("β_lm estimates:", bols_lm)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 3
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function logit(β, X, y)

    # your turn
    diff_u = X*β
    loglike = sum(y.*diff_u - log.(1 .+ exp.(diff_u)))
    loglike = -loglike

    return loglike
end

β_llk = optimize(β -> logit(β, X,y),
                 rand(size(X,2)),
                 # LBFGS(),
                 Newton(),
                 Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(β_llk.minimizer)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 4
#:::::::::::::::::::::::::::::::::::::::::::::::::::
b_glm = glm(@formula(married ~ age + white + collgrad), df, Binomial(), LogitLink())
println("β_glm estimates:", b_glm)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 5
#:::::::::::::::::::::::::::::::::::::::::::::::::::
freqtable(df, :occupation) # note small number of obs in some occupations
df = dropmissing(df, :occupation)
df[df.occupation.==10,:occupation] .= 9
df[df.occupation.==11,:occupation] .= 9
df[df.occupation.==12,:occupation] .= 9
df[df.occupation.==13,:occupation] .= 9
freqtable(df, :occupation) # problem solved

X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
df.white = df.race.==1
y = df.occupation

function mlogit(β, X, y)

    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)

    for j=1:J
        bigY[:,j] = y.==j
    end

    β_big = [reshape(β,K,J-1) zeros(K)]

    numerator = zeros(N,J)
    denominator = zeros(N)

    for j=1:J
        numerator[:,j] = exp.(X*β_big[:,j])
        denominator .+= numerator[:,j]
    end

    P = numerator./repeat(denominator, 1, J)
    loglike = sum(bigY.*log.(P))

    return -loglike
end

    alpha_zero = zeros(6*size(X,2))
    alpha_rand = rand(6*size(X,2))
    alpha_true = [.1910213,-.0335262,.5963968,.4165052,-.1698368,-.0359784,1.30684,-.430997,.6894727,-.0104578,.5231634,-1.492475,-2.26748,-.0053001,1.391402,-.9849661,-1.398468,-.0142969,-.0176531,-1.495123,.2454891,-.0067267,-.5382892,-3.78975]
    alpha_start = alpha_true.*rand(size(alpha_true))
    println(size(alpha_true))
    alpha_hat_optim = optimize(a -> mlogit(a, X, y), alpha_start, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    alpha_hat_mle = alpha_hat_optim.minimizer
    println(alpha_hat_mle)

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # bonus: how to get standard errors?
    # need to obtain the hessian of the obj fun
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # first, we need to slightly modify our objective function
    function mlogit_for_h(alpha, X, y)

        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

        T = promote_type(eltype(X),eltype(alpha)) # this line is new
        num   = zeros(T,N,J)                      # this line is new
        dem   = zeros(T,N)                        # this line is new
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j])
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end

    # declare that the objective function is twice differentiable
    td = TwiceDifferentiable(b -> mlogit_for_h(b, X, y), alpha_start; autodiff = :forward)
    # run the optimizer
    alpha_hat_optim_ad = optimize(td, alpha_zero, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    alpha_hat_mle_ad = alpha_hat_optim_ad.minimizer
    # evaluate the Hessian at the estimates
    H  = Optim.hessian!(td, alpha_hat_mle_ad)
    # standard errors = sqrt(diag(inv(H))) [usually it's -H but we've already multiplied the obj fun by -1]
    alpha_hat_mle_ad_se = sqrt.(diag(inv(H)))
    println([alpha_hat_mle_ad alpha_hat_mle_ad_se]) # these standard errors match Stata

 return nothing

end
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 6
#:::::::::::::::::::::::::::::::::::::::::::::::::::
ps2()


text:
using Oceananigans, Oceananigans.Units, GLMakie

# Specify domain and mesh / grid
Nx = Ny = Nz = 128
Lx = Ly = 128
Lz = 64

grid = RectilinearGrid(GPU(), size = (Nx, Ny, Nz), halo = (3, 3, 3),
                       x = (0, Lx), y = (0, Ly), z = (-Lz, 0),
                       topology = (Periodic, Periodic, Bounded))

# Specify coriolis, viscosity / turbulence closure, and boundary condition
coriolis = FPlane(f₀=1e-4)
closure = AnisotropicMinimumDissipation()
boundary_conditions = (; b = FieldBoundaryConditions(top=FluxBoundaryCondition(1e-8)))

model = NonhydrostaticModel(; grid, coriolis, closure, boundary_conditions,
                              advection = WENO5(),
                              tracers = :b,
                              buoyancy = BuoyancyTracer())

# Time-stepping
simulation = Simulation(model, Δt=10minutes, stop_time=4hours)

progress(sim) = @info "Iter: $(iteration(sim)), time: $(prettytime(sim))"
simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

wizard = TimeStepWizard(cfl=0.1)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

run!(simulation)

# Visualize the x-velocity field
fig, ax, pl = heatmap(interior(model.velocities.w)[:, 1, :])
display(fig)


text:
module test_SolveOneNonlinearEquation

import ..ModiaMath: solveOneNonlinearEquation
using Test
@eval using Printf

@testset "ModiaMath.NonlinearEquations: test solveOneNonlinearEquation" begin
    # Test from Modelica.Math.Nonlinear.Examples.solveNonlinearEquations1
    y          = zeros(3)
    y_analytic = zeros(3)

    fun1(u)       = u^2 - 1
    y[1]          = solveOneNonlinearEquation(fun1, -0.5, 10.0)
    y_analytic[1] = 1.0

    fun2(u;w=1.0) = 3*u - sin(w*u) - 1
    y[2]          = solveOneNonlinearEquation(u->fun2(u; w=3.0), 0.0, 5.0)
    y_analytic[2] = 0.6448544035840080891877538

    fun3(u; p=zeros(2), m=0.0) = p[1] + log(p[2]*u) - m*u
    y[3]          = solveOneNonlinearEquation(u->fun3(u; p=[5.0,1.0], m=1.0), 1.0, 8.0)
    y_analytic[3] = 6.9368474072202187221643182

    for i in 1:length(y)
        println("\n... Results of Solve_SingleNonlinearEquations:")
        println("fun", string(i), ":")
        @printf("   analytical zero     = %20.16e\n", y_analytic[i])
        @printf("   numerical zero      = %20.16e\n", y[i])
        @printf("   absolute difference = %20.16e\n", y_analytic[i] - y[i])
        @test isapprox(y[i], y_analytic[i],  atol=1e-12)
    end

    # Test error case
    errorOccured = false
    try
       yy = solveOneNonlinearEquation(u->fun2(u; w=3.0), 0.0, 0.5)
    catch
       errorOccured = true
    end
    @test errorOccured

    # Test tolerance
    yy = solveOneNonlinearEquation(fun1, -0.5, 10.0; tolerance=1e-4)
    @test isapprox(yy, 1.0,  atol=1e-4)

    # Tests from https://github.com/JuliaMath/Roots.jl/test/test_find_zero.jl
    fns = [(x -> x^2 - exp(x) - 3x + 2,   0.257530285439860, -1, 0.5),
           (x -> x^3 + 4x^2 -10,          1.365230013414097, 0.5, 2.0),
           (x -> exp(x)*sin(x) - 2x - 5, -2.523245230732555, -3.0, -1.0),
           (x -> log(x) + sqrt(x) - 5,    8.309432694231571,  7, 9),
           (x -> sqrt(x) - 1/x - 3,       9.633595562832696,  8, 10),
           (x -> exp(-x) - cos(x), 1.292695719373398, 1.0, 2.0),
           (x -> cos(x)^2 - x/5, 1.085982678007472, 0.5, 2.0),
           (x -> x^10 - x^3 - x - 1, -0.674177935277052, -1.0, 0.0),
           (x -> sin(x) - x + 1, 1.934563210752024, 1.0, 20.0),
           (x -> exp(-x^2 + x + 2) - cos(x) + x^3 + 1, -1.0899423334391694, -2, 0),
           (x -> asin(x^2 -1) - x/2 + 1, 0.594810968398369, 0.3, 0.9),
           (x -> tanh(x) - tan(x), 7.068582745628732, 5.5, 7.1)
    ]

    for (i, (f, u_exact, u_min, u_max)) in enumerate(fns)
        u = solveOneNonlinearEquation(f, u_min, u_max)
        @test isapprox(u, u_exact,  atol=1e-12)
    end

    # Some more tests from https://github.com/JuliaMath/Roots.jl/test/test_find_zero.jl
    # based on http://ir.igsnrr.ac.cn/bitstream/311030/8840/1/%E4%BE%AF%E9%BA%9F%E7%A7%91(SCI)2.pdf
    # which have multiple roots
    multiplicity_tests = [
         (x -> (x - sqrt(5))^4 / ((x-1)^2 + 2),    2.236067977499790, 0, 3)
         (x -> (sin(x)^2 - 2x + 1)^5,              0.71483582544138924, 0, 1)
         (x -> (8x*exp(-x^2) -2x - 3)^8,           -1.7903531791589544, -2, -1)
         (x -> (2x*cos(x) + x^2 - 3)^10/(x^2 + 1),  2.9806452794385368, 0, 3)
         (x -> (exp(-x^2 + x + 3) - x + 2)^9,       2.4905398276083051, 0, 4)

         (x -> (exp(-x) + 2sin(x))^4,              3.1627488709263654, 2, 4)
         (x -> (log(x^2 + 3x + 5) - 2x + 7)^8,     5.4690123359101421, 3, 6)
         (x -> (sqrt(x^2 + 2x + 5) - 2sin(x) - x^2 + 3)^5,   2.3319676558839640, 0, 4)
         (x -> (x-2)^4 / ( (x-1)^2 + 1),           2.0000000000000000, 0, 4)
         (x -> abs(x - 2.5)^(15/4) * exp(x),       2.5, 0, 4)

         (x -> (sqrt(x) - 1/x - 1)^7,              2.147899035704787, 0, 4)
         (x -> (log(x) + sqrt(x) - 5)^3,           8.309432694231572, 0, 10)
         (x -> (sin(x)*cos(x) - x^3 + 1)^9,        1.117078770687451, 0, 2)
         (x -> ((x-3)*exp(x))^5,                   3.0000000000000000, 0, 5)
         (x -> (log(x) + sqrt(x^4 + 1) -2)^7,      1.222813963628973, 1, 3)
         ]

    for (i, (f, u_exact, u_min, u_max)) in enumerate(fns)
        u = solveOneNonlinearEquation(f, u_min, u_max)
        @test isapprox(u, u_exact,  atol=1e-12)
    end
end

end


text:
abstract type AbstractThunk <: AbstractDifferential end

Base.Broadcast.broadcastable(x::AbstractThunk) = broadcastable(extern(x))

@inline function Base.iterate(x::AbstractThunk)
    externed = extern(x)
    element, state = iterate(externed)
    return element, (externed, state)
end

@inline function Base.iterate(::AbstractThunk, (externed, state))
    element, new_state = iterate(externed, state)
    return element, (externed, new_state)
end

#####
##### `Thunk`
#####

"""
    Thunk(()->v)
A thunk is a deferred computation.
It wraps a zero argument closure that when invoked returns a differential.
`@thunk(v)` is a macro that expands into `Thunk(()->v)`.

Calling a thunk, calls the wrapped closure.
`extern`ing thunks applies recursively, it also externs the differial that the closure returns.
If you do not want that, then simply call the thunk

```
julia> t = @thunk(@thunk(3))
Thunk(var"##7#9"())

julia> extern(t)
3

julia> t()
Thunk(var"##8#10"())

julia> t()()
3
```

### When to `@thunk`?
When writing `rrule`s (and to a lesser exent `frule`s), it is important to `@thunk`
appropriately.
Propagation rule's that return multiple derivatives are not able to do all the computing themselves.
 By `@thunk`ing the work required for each, they then compute only what is needed.

#### So why not thunk everything?
`@thunk` creates a closure over the expression, which (effectively) creates a `struct`
with a field for each variable used in the expression, and call overloaded.

Do not use `@thunk` if this would be equal or more work than actually evaluating the expression itself. Examples being:
- The expression wrapping something in a `struct`, such as `Adjoint(x)` or `Diagonal(x)`
- The expression being a constant
- The expression being itself a `thunk`
- The expression being from another `rrule` or `frule` (it would be `@thunk`ed if required by the defining rule already)
"""
struct Thunk{F} <: AbstractThunk
    f::F
end

"""
    @thunk expr

Define a [`Thunk`](@ref) wrapping the `expr`, to lazily defer its evaluation.
"""
macro thunk(body)
    # Basically `:(Thunk(() -> $(esc(body))))` but use the location where it is defined.
    # so we get useful stack traces if it errors.
    func = Expr(:->, Expr(:tuple), Expr(:block, __source__, body))
    return :(Thunk($(esc(func))))
end

"""
    unthunk(x)

On `AbstractThunk`s this removes 1 layer of thunking.
On any other type, it is the identity operation.

In contrast to `extern` this is nonrecursive.
"""
@inline unthunk(x) = x

@inline extern(x::AbstractThunk) = extern(unthunk(x))

# have to define this here after `@thunk` and `Thunk` is defined
Base.conj(x::AbstractThunk) = @thunk(conj(unthunk(x)))

(x::Thunk)() = x.f()
@inline unthunk(x::Thunk) = x()

Base.show(io::IO, x::Thunk) = println(io, "Thunk($(repr(x.f)))")

"""
    InplaceableThunk(val::Thunk, add!::Function)

A wrapper for a `Thunk`, that allows it to define an inplace `add!` function,
which is used internally in `accumulate!(Δ, ::InplaceableThunk)`.

`add!` should be defined such that: `ithunk.add!(Δ) = Δ .+= ithunk.val`
but it should do this more efficently than simply doing this directly.
(Otherwise one can just use a normal `Thunk`).

Most operations on an `InplaceableThunk` treat it just like a normal `Thunk`;
and destroy its inplacability.
"""
struct InplaceableThunk{T<:Thunk, F} <: AbstractThunk
    val::T
    add!::F
end

unthunk(x::InplaceableThunk) = unthunk(x.val)
(x::InplaceableThunk)() = unthunk(x)

function Base.show(io::IO, x::InplaceableThunk)
    println(io, "InplaceableThunk($(repr(x.val)), $(repr(x.add!)))")
end

# The real reason we have this:
accumulate!(Δ, ∂::InplaceableThunk) = ∂.add!(Δ)
store!(Δ, ∂::InplaceableThunk) = ∂.add!((Δ.*=false))  # zero it, then add to it.


text:
# Load packages and helper functions
using MPI, Statistics
include("../parallelhelpers.jl")

# Call MPI.init() and
MPI.Init()
comm = MPI.COMM_WORLD
my_rank = MPI.Comm_rank(comm)
n_procs = MPI.Comm_size(comm)
root = 0

# Load the file names
dataLoc = "../../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Distribute filenames by rank
myfnames = fnames[my_rank+1:n_procs:length(fnames)]

# Calculate the word counts
mycounts =  getcounts.(myfnames)

# Gather all counts onto Process 0 and merge
if my_rank > 0
    # send to 0
    msg = MPI.serialize(mycounts)
    print("$my_rank: Sending mycounts $my_rank -> 0\n")
    sreq = MPI.Send(length(msg), 0, my_rank+n_procs, comm) # First send the length of the message
    sreq = MPI.Send(msg, 0, my_rank+(n_procs*2), comm) # Then send the message
else
    allcounts = Array{Array{UInt8,1},1}(undef,n_procs-1)
    for i = 1:n_procs-1
        # receive
        global allcounts
        msglen = Array{Int64,1}(undef,1)
        println("$my_rank: Waiting to receive length from $i")
        rreq = MPI.Recv!(msglen, i,  i+n_procs, comm) # First receive the length of the message
        allcounts[i] = Array{UInt8}(undef,msglen[1])
        println("$my_rank: Waiting to receive msg from $i")
        rreq = MPI.Recv!(allcounts[i], i,  i+(n_procs*2), comm) # Teh receive the message
        println("$my_rank: Received from $i")
    end
    allcounts = vcat(MPI.deserialize.(allcounts)...)
    overallcounts = merge(+,mycounts...,allcounts...)
end

# Now send the overall counts to all the other Processes
if my_rank == 0
    msg = MPI.serialize(overallcounts)
    for i in 1:n_procs-1
        global msg
        print("$my_rank: Sending overallcounts $my_rank -> $i\n")
        sreq = MPI.Send(length(msg), i, i+(n_procs*3), comm) # First send the length of the message
        sreq = MPI.Send(msg, i, i+(n_procs*4), comm) # Then send the message
    end
else
    msglen = Array{Int64,1}(undef,1)
    println("$my_rank: Waiting to receive length from 0")
    rreq = MPI.Recv!(msglen, 0,  my_rank+(n_procs*3), comm)
    overallcounts = Array{UInt8}(undef,msglen[1])
    println("$my_rank: Waiting to receive msg from 0")
    rreq = MPI.Recv!(overallcounts, 0,  my_rank+(n_procs*4), comm)
    overallcounts = MPI.deserialize(overallcounts)
end

# Calculate the normalized counts
normcount = getnormcount.(mycounts,Ref(overallcounts))

# Print out the top 5 words for each file and their counts on the leader process
for i=1:length(normcount)
    top5 = sort(collect(normcount[i]), by=x->x[2], rev=true)[1:min(5,length(normcount[i])),:]
    println(top5)
end

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
# ---------------------------------------- Gravity ----------------------------------------
# Define abstract type and interface.
abstract type Gravity end

# TODO: make an interface which receives position. Right now, for flat Earth and constant
# gravity does not make sense. However, will be usefull in the future.
"""
    get_gravity_accel(g::Gravity)

Get gravity acceleration value (m/s²).
"""
function get_gravity_accel(g::Gravity) end

"""
    get_gravity_horizon(g::Gravity)

Get gravity vector in local horizon coordinates (North, East, Down) (m/s²).
"""
function get_gravity_horizon(g::Gravity) end

"""
    get_gravity_horizon(g::Gravity, θ, ϕ)

Get gravity vector in body axis (m/s²).
"""
function get_gravity_body(g::Gravity, θ, ϕ)
    ax, ay, az = horizon2body(get_gravity_horizon(g)..., 0.0, θ, ϕ)
    return @SVector [ax, ay, az]
end

# ------------------------------ Local Horizon Down Gravity -------------------------------
"""
    LHDownGravity{T<:Number} <: Gravity

Local Horizon Down gravity.

# Fields
- `gD::T`: down gravity acceleration value (m/s²).
"""
struct LHDownGravity{T<:Number} <: Gravity
    """Gravity acceleration (m/s²)"""
    gD::T
end

"""
    LHDownGravity([gD])

If no value of acceleration is given, default acceleration is: `gD=$(gD) m/s²`.
"""
LHDownGravity(gD=gD) = LHDownGravity(gD)

get_gravity_accel(g::LHDownGravity) = g.gD
get_gravity_horizon(g::LHDownGravity) = @SVector [0, 0, g.gD]


text:
#=
Copyright 2021 BlackRock, Inc.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#
const ϵ = 1e-12

function ≈(a::Real, b::Real)
    return (abs(a - b) ≤ ϵ) || (a == b)  # Second case covers infinite values.
end
≉(a::Real, b::Real) = !(a ≈ b)

function ≲(a::Real, b::Real)
    return a ≤ b + ϵ
end
≳(a::Real, b::Real) = b ≲ a
⋦(a::Real, b::Real) = !((a ≳ b) || isnan(a) || isnan(b))
⋧(a::Real, b::Real) = b ⋦ a

"""
    solve_quad(a::Real, b::Real, c::Real)

Solve a quadratic function `a x^2 + b x + c` using the quadratic formula.
"""
function solve_quad(a::Real, b::Real, c::Real)
    if a == 0
        if b == 0.0
            return NaN, NaN
        else
            return -c / b, NaN
        end
    else
        b2m4ac = b^2 - 4 * a * c
        if !(b2m4ac ≳ 0.0)
            return NaN, NaN
        end
        sqrt_b2m4ac = sqrt(pos(b2m4ac))
        if b > 0
            x1 = (-b - sqrt_b2m4ac) / (2 * a)
            x2 = 2 * c / (-b - sqrt_b2m4ac)
        else
            x1 = (-b + sqrt_b2m4ac) / (2 * a)
            x2 = 2 * c / (-b + sqrt_b2m4ac)
        end
        return x1, x2
    end
end

pos(a::Real) = max(a, 0)
neg(a::Real) = pos(-a)
clip(a::Real, lb::Real, ub::Real) = min(max(a, lb), ub)
clip(a::Real, dom::Interval) = clip(a, dom.lb, dom.ub)


text:
##### Beginning of file

module Hosts # Begin submodule Snapshots.Hosts

__precompile__(true)

import ..Types
import ..Utils

import ..delayederror

include(joinpath("BitbucketSrcHost", "BitbucketSrcHost.jl"))
include(joinpath("BitbucketDstHost", "BitbucketDstHost.jl"))

include(joinpath("GitHubSrcHost", "GitHubSrcHost.jl"))
include(joinpath("GitHubDstHost", "GitHubDstHost.jl"))

include(joinpath("GitLabSrcHost", "GitLabSrcHost.jl"))
include(joinpath("GitLabDstHost", "GitLabDstHost.jl"))

end # End submodule Snapshots.Hosts

##### End of file


text:
# Abstract Rational-Quadratic Kernel =======================================================
abstract type AbstractRationalQuadraticKernel{T<:AbstractFloat} <: MercerKernel{T} end

@inline basefunction(::AbstractRationalQuadraticKernel) = SquaredEuclidean()

# Rational-Quadratic Kernel ================================================================
@doc raw"""
    RationalQuadraticKernel([α [,β]])

The rational quadratic kernel is an isotropic Mercer kernel given by the formula:

```
    κ(x,y) = (1 + α‖x-y‖²)⁻ᵝ   α > 0, β > 0
```
where `α` is a scaling parameter and `β` is a shape parameter. The rational quadratic
kernel is a special of the more general gamma-rational-quadratic kernel (see
[`GammaRationalQuadraticKernel`](@ref)) with `γ = 1`.

# Examples

```jldoctest; setup = :(using MLKernels)
julia> RationalQuadraticKernel()
RationalQuadraticKernel{Float64}(1.0,1.0)

julia> RationalQuadraticKernel(2.0f0)
RationalQuadraticKernel{Float32}(2.0,1.0)

julia> RationalQuadraticKernel(2.0f0, 2.0)
RationalQuadraticKernel{Float64}(2.0,2.0)
```
"""
struct RationalQuadraticKernel{T<:AbstractFloat} <: AbstractRationalQuadraticKernel{T}
    α::T
    β::T
    function RationalQuadraticKernel{T}(
            α::Real=T(1),
            β::Real=T(1)
        ) where {T<:AbstractFloat}
        @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 0")
        @check_args(RationalQuadraticKernel, β, β > zero(T), "β > 0")
        return new{T}(α, β)
    end
end
function RationalQuadraticKernel(
        α::T₁=1.0,
        β::T₂=T₁(1)
    ) where {T₁<:Real, T₂<:Real}
    RationalQuadraticKernel{promote_float(T₁, T₂)}(α, β)
end

@inline function kappa(κ::RationalQuadraticKernel{T}, d²::T) where {T}
    return (one(T) + κ.α*d²)^(-κ.β)
end

function convert(
        ::Type{K},
        κ::RationalQuadraticKernel
    ) where {K>:RationalQuadraticKernel{T}} where T
    return RationalQuadraticKernel{T}(κ.α, κ.β)
end

# Gamma Rational-Quadratic Kernel ==========================================================
@doc raw"""
    GammaRationalKernel([α [,β [,γ]]])

The gamma-rational-quadratic kernel is a generalization of the rational-quadratic kernel
with an additional shape parameter `γ`:

```
    κ(x,y) = (1 + α‖x-y‖²ᵞ)⁻ᵝ   α > 0, β > 0, γ ∈ (0,1]
```
where ``\alpha`` is a scaling parameter, ``\beta`` is a shape parameter and ``\gamma`` is a
shape parameter of the Euclidean distance.

# Examples

```jldoctest; setup = :(using MLKernels)
julia> GammaRationalQuadraticKernel()
GammaRationalQuadraticKernel{Float64}(1.0,1.0,1.0)

julia> GammaRationalQuadraticKernel(2.0f0)
GammaRationalQuadraticKernel{Float32}(2.0,1.0,1.0)

julia> GammaRationalQuadraticKernel(2.0f0, 2.0f0)
GammaRationalQuadraticKernel{Float32}(2.0,2.0,1.0)

julia> GammaRationalQuadraticKernel(2.0f0, 2.0f0, 0.5f0)
GammaRationalQuadraticKernel{Float32}(2.0,2.0,0.5)
```
"""
struct GammaRationalQuadraticKernel{T<:AbstractFloat} <: AbstractRationalQuadraticKernel{T}
    α::T
    β::T
    γ::T
    function GammaRationalQuadraticKernel{T}(
            α::Real=T(1),
            β::Real=T(1),
            γ::Real=T(1)
        ) where {T<:AbstractFloat}
        @check_args(GammaRationalQuadraticKernel, α, α > zero(T), "α > 0")
        @check_args(GammaRationalQuadraticKernel, β, β > zero(T), "β > 0")
        @check_args(GammaRationalQuadraticKernel, γ, one(T) >= γ > zero(T), "γ ∈ (0,1]")
        return new{T}(α, β, γ)
    end
end
function GammaRationalQuadraticKernel(
        α::T₁ = 1.0,
        β::T₂ = T₁(1),
        γ::T₃ = one(promote_float(T₁, T₂))
    ) where {T₁<:Real, T₂<:Real, T₃<:Real}
    GammaRationalQuadraticKernel{promote_float(T₁,T₂,T₃)}(α, β, γ)
end

@inline function kappa(κ::GammaRationalQuadraticKernel{T}, d²::T) where {T}
    return (one(T) + κ.α*(d²^κ.γ))^(-κ.β)
end

function convert(
        ::Type{K},
        κ::GammaRationalQuadraticKernel
    ) where {K>:GammaRationalQuadraticKernel{T}} where T
    return GammaRationalQuadraticKernel{T}(κ.α, κ.β, κ.γ)
end


text:
"""
This module provides a collection of multivariate polynomial bases.

The exported names are:

$(EXPORTS)
"""
module Polynomials

using DocStringExtensions
using LinearAlgebra: mul!
using InplaceArrays.Helpers
using InplaceArrays.Inference
using InplaceArrays.Arrays
using InplaceArrays.TensorValues
using InplaceArrays.Fields

import InplaceArrays.Fields: evaluate_field!
import InplaceArrays.Fields: field_cache
import InplaceArrays.Fields: evaluate_gradient!
import InplaceArrays.Fields: gradient_cache
import InplaceArrays.Fields: evaluate_hessian!
import InplaceArrays.Fields: hessian_cache

export MonomialBasis
export QGradMonomialBasis
export QCurlGradMonomialBasis
export change_basis

include("MonomialBases.jl")

include("QGradMonomialBases.jl")

include("QCurlGradMonomialBases.jl")

include("ChangeBasis.jl")

end # module


text:
 function solveBGP(p,eqInit)

  eq = EqObj()   # initialize equilibrium object

  # gausslegendre for integration
  eq.qMinAll = 1.0e-8
  eq.qMaxAll = 10.0
  eq.node,eq.weight = gausslegendre(1000)
  eq.node = ((eq.qMaxAll - eq.qMinAll)/2)*eq.node .+ (eq.qMaxAll + eq.qMinAll)/2    # adjust the limits

  @inline function objFnc(eqnd,x)
     eqfunc!(eqnd,x,eq,p)
  end

  res = nlsolve(objFnc, eqInit, method = :trust_region,inplace = true,show_trace=false,iterations=50)

  if !res.f_converged
       #print("👎")
       #println("Couldn't solve!!!")
       eq = EqObj()  # return empty EqObj, type stability
  else
       #print("👍")
       eqfunc!(similar(res.zero),res.zero,eq,p)   # evaluate at minimizer
  end

  return eq,res

end

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

function eqfunc!(eqnd,eqv,eq,p)
  # reject out of bounds

  if (any(eqv.<=0.0))
    eqnd = -1000.0*ones(size(eqv)[1])

  else

    # load in guesses for solution
    eq.ws      = eqv[1]
    eq.cactiv  = eqv[2:3]
    eq.eyq     = eqv[4:5]
    eq.qbar    = eqv[6]

    # flow of free products
    eq.cactivtot = sum(eq.cactiv)
    eq.cinac     = 1.0 - eq.cactivtot

    # find updated eqvars
    innovation!(eq,p)
    qualityDist!(eq,p)
    qbarActive!(eq,p)
    calcey!(eq,p)
    labordem!(eq,p)

    # equilibrium differences

    eqnd[1]    = eq.skilled_lab - p.Lˢ           # skilled labor market clearing
    eqnd[2:3] .= eq.cactiv     .- eqv[2:3]       # product shares
    eqnd[4:5] .= eq.eyq        .- eqv[4:5]       # innovated product value
    eqnd[6]    = eq.qbarAct     - (p.ε/(p.ε-1.0))^(p.ε - 1)

  end
end

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

function innovation!(eq,p)

  # innovation rates
  eq.x    = p.𝚯.* (((1.0 - p.γ)/eq.ws)*eq.eyq).^((1.0 - p.γ)/p.γ)
  eq.xout = p.θᴱ.*(((1.0 - p.γᴱ)/eq.ws).*sum(p.𝛂.*eq.eyq)).^((1.0 - p.γᴱ)/p.γᴱ)

  # option value of R&D
  eq.optval = p.γ*p.𝚯*((1.0 - p.γ)/eq.ws)^((1.0 - p.γ)/p.γ).*eq.eyq.^(1.0/p.γ)

  # aggregate creative destruction rates
  eq.τs = eq.cactiv.*eq.x .+ p.𝛂*eq.xout
  eq.τ  = sum(eq.τs)

  # minimum qhat
  eq.qmin = (max.(0,eq.ws*p.ϕ .- eq.optval)/p.Π).^(1.0/(p.ε-1.0))

end

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

function qualityDist!(eq,p)

    # -1 - Parameters for delayed differential equation
    alphaDelay = 1.0/(1.0 + (p.λ-1.0)*(1.0 - p.ω))
    betaDelay  = (p.λ - 1.0)*p.ω*eq.qbar*alphaDelay

    # 0 - Growth rate consistent with guessed qbar
    eq.g = eq.τ*(p.λ - 1.0)

    # 1 - Overall dist.function
    function funcFAll!(du,u,h,p,t)
        du[1] = (eq.τ/(eq.g*t))*u[1] - (eq.τ/(eq.g*t))*h(p,t-betaDelay)[1]
    end

    tspan = (1.0e-09,10.0)
    lags = [betaDelay]
    h(p,t)= [0.0]
    probDelay  = DDEProblem(funcFAll!,[1.0e-16],h,tspan, constant_lags = lags)
    eq.solFAll = solve(probDelay,MethodOfSteps(Tsit5()),reltol=1e-8,abstol=1e-8) # to-do: we can directly get output at the nodes by using ``saveat``
    # option and pass it to integration functions. One thing to note: once saveat is used, between node results are linearly interpolated.

    eq.FAllend = eq.solFAll.u[end][1,1]  # needed for normalization

    eq.qAllDist  = vcat_nosplat(eq.solFAll(eq.node,Val{1}),1).*(1/eq.FAllend) # saving for later use

    # 2 -  Gross dist. and active product line dist.
    # different than matlab verison, we convert the problem to IVP by simply reversing time span
    eq.PhiHG = eq.τs[2]/(eq.τ + p.ν)
    eq.PhiLG = 1.0 - eq.PhiHG

    @inline function getFAllVal(t)
      delayT 	= min(max(0.0000001,t - betaDelay),10.0)
      fAll = eq.solFAll(t)[1,1]/eq.FAllend
      fAllDelay = eq.solFAll(delayT)[1,1]/eq.FAllend
      return fAll::Float64,fAllDelay::Float64
    end

    function funcFRest!(du,u,pp,t)
      fAll,fAllDelay = getFAllVal(t)

      du[1] = ((eq.τ + p.ν)/(eq.g*t))*u[1] - (eq.τs[2]/(eq.g*t))*fAllDelay
      du[2] = ((eq.τ + p.ν + p.ψ)/(eq.g*t))*u[2] - (p.ψ/(eq.g.*t))*u[1]
      du[3] = ((eq.τ +       p.ψ)/(eq.g*t))*u[3] - (p.ψ/(eq.g*t))*(fAll - u[1]) - (p.ν/(eq.g*t))*u[2]
    end

    boundHRest = (p.ψ/(eq.τ + p.ν + p.ψ))*eq.PhiHG
    boundLRest = (p.ψ/(eq.τ + p.ψ))*eq.PhiLG + (p.ν/(eq.τ + p.ψ))*boundHRest

    u0 = [eq.PhiHG,boundHRest,boundLRest]
    probFRest   = ODEProblem(funcFRest!,u0,(tspan[2],tspan[1]))
    eq.solFRest = solve(probFRest, Tsit5(),reltol=1e-8,abstol=1e-8)

    # 3 - Update PhiL and PhiH
    FRestcut    = eq.solFRest(max.(1.0e-09,eq.qmin))
    FAllcut     = eq.solFAll(max.(1.0e-09,eq.qmin))./eq.FAllend

    eq.PhiHExo   = eq.solFRest(10.0)[2]
    eq.PhiLExo   = eq.solFRest(10.0)[3]

    eq.PhiHObso  = FRestcut[1,2] - FRestcut[2,2]
    eq.PhiLObso  = (FAllcut[1,1] - FRestcut[1,1]) - FRestcut[3,1]

    eq.cactiv[1] = (eq.PhiLG - eq.PhiLExo) - eq.PhiLObso
    eq.cactiv[2] = (eq.PhiHG - eq.PhiHExo) - eq.PhiHObso
    eq.cactivtot = sum(eq.cactiv)
    eq.cinac     = 1.0 - eq.cactivtot

end

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

function qbarActive!(eq,p)
  outL,outH = funQbarAct(eq,p)
  eq.qbarAct  = ((eq.qMaxAll - eq.qMinAll)/2)*(dot(outL,eq.weight) + dot(outH,eq.weight))
end
########

function funQbarAct(eq,p)
        outRest  = eq.solFRest(eq.node,Val{1})

        outL = (eq.qAllDist .- vcat_nosplat(outRest,1) .- vcat_nosplat(outRest,3)).*(eq.node.^(p.ε-1)).*(eq.node.>=eq.qmin[1])
        outH = (vcat_nosplat(outRest,1)                .- vcat_nosplat(outRest,2)).*(eq.node.^(p.ε-1)).*(eq.node.>=eq.qmin[2])
        return outL, outH
end

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

function calcey!(eq,p)

  eq.r = p.ρ + p.σ*eq.g

  outL, outH = eyfunc(eq,p)

  eq.eyq = [((eq.qMaxAll - eq.qMinAll)/2)*dot(outL,eq.weight), ((eq.qMaxAll - eq.qMinAll)/2)*dot(outH,eq.weight)]

end
###################################################
function eyfunc(eq,p)

  qPlus = eq.node .+ (p.λ - 1).*(p.ω*eq.qbar .+ (1 - p.ω).*eq.node)

  eyL = zfunc(qPlus,0.0,1,eq,p).*eq.qAllDist     # LOW TYPE
  eyH = (zfunc(qPlus,p.ν,2,eq,p) .+ zfunc(qPlus,0.0,1,eq,p) .- zfunc(qPlus,p.ν,1,eq,p)).*eq.qAllDist      # HIGH TYPE

  return eyL, eyH

end

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

function zfunc(q, x, s, eq, p)

  psit = eq.r + eq.τ + p.ψ + x

  kappa1 = psit + (p.ε - 1.0)*eq.g
  kappa2 = psit

  coeff1 = p.Π
  coeff2 = eq.optval[s] - eq.ws*p.ϕ

  # this min means that when q <= qmin, qrat = 1.0, so the value contribution is zero
  qrat = min.(1.0,eq.qmin[s]./q)

  zout = coeff1/kappa1.*q.^(p.ε-1.0).*(1.0 .- qrat.^(kappa1/eq.g)) .+ coeff2/kappa2.*(1.0 .- qrat.^(kappa2/eq.g))

end

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

function labordem!(eq,p)

  # managerial fixed cost
  eq.cfix = eq.cactivtot*p.ϕⁿˢ

  # production labor
  eq.wu = (p.ε - 1.0)/p.ε   # unskilled wage over Q.

  # incumbent R&D labor
  eq.cx = (eq.x./(p.𝚯ⁿˢ.^p.γ)).^(1.0/(1.0-p.γ))
  eq.crnd = sum(eq.cactiv.*eq.cx)

  # entrant R&D labor
  eq.cout = (eq.xout/(p.θᴱⁿˢ^p.γᴱ))^(1.0/(1.0-p.γᴱ))

  # total skilled labor
  eq.skilled_lab = eq.crnd + eq.cfix + eq.cout
end

#########################################################
# Utility Functions

vcat_nosplat(y,l) =  eltype(y[l])[el[l] for el in y] # this is for converting weird output of diff solver