text:
#=
demo_dHDP_Gaussian1DGaussian1D

A demo for dynamic Hierarchical Dirichlet Process mixture models with Gaussian1DGaussian1D
Bayesian components.

23/06/2016
Adham Beyki, odinay@gmail.com
=#

using BNP
srand(123)

## --- synthesizing the data --- ##
#=
true_gg = 2.0
true_aa = 0.5
true_aw = 2
true_bw = 3
n_groups = 1000
n_group_j   = 100 * ones(Int, n_groups)
KK_truncation = 150

beta_tilde = rand(Distributions.Beta(1, true_gg), KK_truncation)
beta_tilde[KK_truncation] = 1.0
my_beta = BNP.stick_breaking(beta_tilde)

pi_tilde = zeros(Float64, n_groups, KK_truncation)
my_pi = zeros(Float64, n_groups, KK_truncation)
for jj = 1:n_groups
    for kk = 1:KK_truncation-1
        gamma1 = true_aa * my_beta[kk]
        gamma2 = true_aa * (1 - sum(my_beta[1:kk]))
        pi_tilde[jj, kk] = rand(Distributions.Beta(gamma1, gamma2))
    end
    pi_tilde[jj, KK_truncation] = 1.0
    my_pi[jj, :] = BNP.stick_breaking(pi_tilde[jj, :][:])
end

w_tilde = rand(Distributions.Beta(true_aw, true_bw), n_groups-1)
insert!(w_tilde, 1, 1)

true_rr = Array(Vector{Int}, n_groups)
true_zz = Array(Vector{Int}, n_groups)

logww_prime = [0.0]
for jj = 1:n_groups
    if jj == 1
        logww = [0.0]
    else
        logww = zeros(Float64, jj)
        logww[1:jj-1] = logww_prime + log(1 - w_tilde[jj])
        logww[jj] = log(w_tilde[jj])
        logww_prime = logww
    end
    pp_rr = exp(logww)

    true_rr[jj] = zeros(Int, n_group_j[jj])
    true_zz[jj] = zeros(Int, n_group_j[jj])
    for ii = 1:n_group_j[jj]
        true_rr[jj][ii] = sample(pp_rr)
        kk = sample(my_pi[true_rr[jj][ii], :][:])
        true_zz[jj][ii] = kk
    end
end

zz_flat = Int[]
for jj = 1:n_groups
    append!(zz_flat, true_zz[jj])
end

zz_unq = sort(unique(zz_flat))
true_KK = length(zz_unq)
for jj = 1:n_groups
    for kk = 1:true_KK
        idx = find(x -> x==zz_unq[kk], true_zz[jj])
        true_zz[jj][idx] = kk
    end
end

vv = 0.001          # fixed variance
ss = 2
true_atoms = [Gaussian1D(ss*kk, vv) for kk = 1:true_KK]

xx      = Array(Vector{Float64}, n_groups)
true_nn = zeros(Int, n_groups, true_KK)

for jj = 1:n_groups

    xx[jj] = zeros(Float64, n_group_j[jj])
    for ii = 1:n_group_j[jj]
        kk = true_zz[jj][ii]
        xx[jj][ii] = sample(true_atoms[kk])
        true_nn[jj, kk] += 1
    end
end
=#

srand(123)

true_gg     = 2.0
true_aa     = 0.5
n_groups    = 100
n_group_j   = 100 * ones(Int, n_groups)
join_tables = true

true_tji, true_njt, true_kjt, true_nn, true_mm, true_zz, true_KK = BNP.gen_CRF_data(n_group_j, true_gg, true_aa, join_tables)

vv = 0.001          # fixed variance
ss = 2
true_atoms = [Gaussian1D(ss*kk, vv) for kk = 1:true_KK]

true_nn = zeros(Int, n_groups, true_KK)
xx = Array(Vector{Float64}, n_groups)
for jj = 1:n_groups
    xx[jj] = zeros(Float64, n_group_j[jj])
    for ii = 1:n_group_j[jj]
        kk = true_zz[jj][ii]
        xx[jj][ii] = sample(true_atoms[kk])
        true_nn[jj, kk] += 1
    end
end

## ------- inference -------- ##
# constructing the Bayesian component of LDA model
m0 = mean(mean(xx))
v0 = 10.0
q0 = Gaussian1DGaussian1D(m0, v0, vv)

# constructing the HDP model
dhdp_KK_init = 1
dhdp_gg = 2.0
dhdp_g1 = 0.1
dhdp_g2 = 0.1
dhdp_aa = 0.5
dhdp_a1 = 0.1
dhdp_a2 = 0.1
dhdp_aw = 2.0
dhdp_bw = 3.0
dhdp = dHDP(q0, dhdp_KK_init, dhdp_gg, dhdp_g1, dhdp_g2, dhdp_aa, dhdp_a1, dhdp_a2, dhdp_aw, dhdp_bw)

# sampling
zz = Array(Vector{Int}, n_groups)
for jj = 1:n_groups
    zz[jj] = ones(Int, n_group_j[jj])
    zz[jj] = rand(1:dhdp.KK, n_group_j[jj])
end

n_burnins   = 500
n_lags      = 1
n_samples   = 500
sample_hyperparam = true
n_internals = 10
store_every = 1000
results_path = ""
filename    = "demo_dHDP_Gaussian1DGaussian1D_"

KK_list, KK_dict, w_tilde = truncated_gibbs_sampler(dhdp, xx, zz,
    n_burnins, n_lags, n_samples, sample_hyperparam, n_internals, store_every, results_path, filename)

KK_hist = hist(KK_list, 0.5:maximum(KK_list)+0.5)[2]
candidate_KK = indmax(KK_hist)
posterior_components, nn = posterior(dhdp, xx, KK_dict[candidate_KK], candidate_KK)


text:
@doc """
    ≺( X, Y )

Return `true` if vector `X` dominates `Y`, i.e. if all components
of `X` not bigger than those of `Y` and if there is at least one
component of `X` that is strictly smaller than the corresponding
`Y` component.
"""
function ≺( X :: Vector{R} where R<:Real, Y :: Vector{R} where R<:Real )
    if length(X) != length(Y)
        error("Cannot compare vectors of unequal lengths.")
    end

    for i = eachindex(X)
        if Y[i] < X[i]
            return false
        end
    end
    return any( X .< Y )
end
#=
@doc """
    ≺( box1, box2 )

Return `true` if all result vectors (indexed by the respective fields
`db_indices`) of `box1` dominate all result vectors in `box2`,
else return `false`.
"""
function ≺( box1 :: Box, box2 :: Box )
    for X ∈ results( box1; include_nothing = false )
        for Y ∈ results( box2; include_nothing = false )
            if !≺(X,Y)
                return false
            end
        end
    end
    return true
end

@doc """
    ≺( X, box2 )

Return `true` if `X` dominates all result vectors of `box2`.
"""
function ≺( X :: Vector{R} where R<:Real, box2 :: Box )
    for Y ∈ results( box2; include_nothing = false )
        if !≺(X,Y)
            return false
        end
    end
    return true
end

"Return `true` if there is no result vector of `box2` that dominates `X`."
function point_not_dominated_by_box(  X :: Vector{R} where R<:Real, box2 :: Box )
    for Y ∈ results( box2; include_nothing = false )
        if ≺(Y, X)
            return false
        end
    end
    return true
end
=#

ext: jl


text:
module TestWiringDiagramLayouts

using Test

using Catlab.Theories, Catlab.WiringDiagrams
using Catlab: Graphics
using Catlab.Graphics.WiringDiagramLayouts: LeftToRight, position

layout_diagram(f; kw...) = Graphics.layout_diagram(f; outer_ports_layout=:fixed, kw...)
box_layouts(d::WiringDiagram) = map(box -> box.value, boxes(d))
box_values(d::WiringDiagram) = map(layout -> layout.value, box_layouts(d))

# Generators
A, B = Ob(FreeSymmetricMonoidalCategory, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

d = layout_diagram(f, orientation=LeftToRight)
@test box_values(d) == [f]
fv = first(box_ids(d))
((f_in_x, f_in_y),) = map(position, input_ports(box(d,fv)))
((f_out_x, f_out_y),) = map(position, output_ports(box(d,fv)))
@test f_in_x < f_out_x
@test f_in_y == f_out_y

# Composition
d = layout_diagram(compose(f,g), orientation=LeftToRight)
@test box_values(d) == [f,g]
(fx, fy), (gx, gy) = map(position, box_layouts(d))
@test fx < gx
@test fy == gy

# Products
d = layout_diagram(otimes(f,g), orientation=LeftToRight)
@test box_values(d) == [f,g]
(fx, fy), (gx, gy) = map(position, box_layouts(d))
@test fx == gx
@test fy < gy

# Identities and braidings
@test box_values(layout_diagram(id(A))) == []
@test box_values(layout_diagram(braid(A,B))) == []

# Diagonals and codiagonals
A, B, C = Ob(FreeBiproductCategory, :A, :B, :C)

for expr in (mcopy(A), delete(A), mmerge(A), create(A))
  @test map(l -> l.shape, box_layouts(layout_diagram(expr))) == [:junction]
end

for expr in (mcopy(A⊗B), delete(A⊗B), mmerge(A⊗B), create(A⊗B))
  d = layout_diagram(expr)
  @test map(l -> l.shape, box_layouts(d)) == repeat([:junction], 2)
end

for expr in (mcopy(A⊗B⊗C), delete(A⊗B⊗C), mmerge(A⊗B⊗C), create(A⊗B⊗C))
  d = layout_diagram(expr)
  @test map(l -> l.shape, box_layouts(d)) == repeat([:junction], 3)
end

end


text:
using Lifting
using Test
cd(@__DIR__)

SampleExercise = Exercise(;
    name = "Swiss Bar Bench",
    modality = ["Block", "Semi-Supinated"],
    equipment = "Swiss Bar",
    size = "1 inch",
    muscles = ["Pecs", "Triceps", "Front Delts"],
    trainingMax = 100,
    roundBase = 2.5,
    roundMode = floor,
)
SampleExercise2 = Exercise(;
    name = "Squat",
    modality = "Pause",
    trainingMax = 130,
    roundBase = 2.5,
    roundMode = floor,
)
SampleScheme = SetScheme(;
    type = ["Long Rest", "Longer Rest", "Longest Rest"],
    sets = [1, 2, 1],
    reps = [12, 14, 10],
    intensity = [7, 8.5, 9.5],
    roundMode = [floor, floor, ceil],
    rpeMode = true,
)
SampleScheme2 = SetScheme(;
    type = ["Long Rest", "Longer Rest", "Longest Rest"],
    sets = [1, 2, 1],
    reps = [12, 14, 10],
    intensity = [7, 8.5, 9.5],
    roundMode = [floor, floor, ceil],
    rpeMode = true,
)
SampleScheme3 =
    SetScheme(; type = "Long Rest", sets = 1, reps = 12, intensity = 7, rpeMode = true)
SampleProgression = Progression(;
    type = LinearProgression(),
    name = "Progression Name",
    setScheme = SampleScheme,
)
SampleProgression2 = Progression(;
    type = LinearProgression(),
    name = "Progression Name",
    sessions = 1,
    period = 2,
    setScheme = [deepcopy(SampleScheme), SampleScheme2],
)
SampleProgression3 = Progression(;
    type = LinearProgression(),
    name = "Progression Name",
    setScheme = SampleScheme3,
)
struct SampleProgramme <: AbstractProgramme end
exerProg = Dict()
push!(
    exerProg,
    "SampleExercise" => (exercise = SampleExercise, progression = SampleProgression),
)
push!(
    exerProg,
    "SampleExercise2" => (exercise = SampleExercise2, progression = SampleProgression2),
)
push!(
    exerProg,
    "SampleExercise3" => (exercise = SampleExercise2, progression = SampleProgression3),
)

import Lifting: makeDays
function makeDays(::SampleProgramme, exerProg::Dict)
    week = [[], [], [], [], [], [], []]
    SampleExercise = exerProg["SampleExercise"].exercise
    SampleProgression = exerProg["SampleExercise"].progression

    SampleExercise2 = exerProg["SampleExercise2"].exercise
    SampleProgression2 = exerProg["SampleExercise2"].progression

    SampleExercise3 = exerProg["SampleExercise3"].exercise
    SampleProgression3 = exerProg["SampleExercise3"].progression

    push!(week[1], SampleExercise, SampleProgression)   # Day 1
    push!(week[1], SampleExercise, SampleProgression)   # Day 1
    insert!(week[1], 1, SampleExercise2, SampleProgression2)   # Day 1
    insert!(week[1], 3, SampleExercise3, SampleProgression3)   # Day 1
    push!(week[2], "Rest")  # Day 2
    push!(week[3], SampleExercise, SampleProgression)   # Day 3
    push!(week[4], "Rest")  # Day 4
    push!(week[5], SampleExercise3, SampleProgression3)   # Day 5
    push!(week[6], "Rest")  # Day 6
    push!(week[7], "Rest")  # Day 7
    return week
end

week = makeDays(SampleProgramme(), exerProg)
sampleProgramme = Programme(SampleProgramme(), "sampleProgramme", exerProg, week)

@test println(sampleProgramme[1]) === println(sampleProgramme, 1)
@test println(sampleProgramme) === println(sampleProgramme, 1:7)
@test println(sampleProgramme[1:3]) === println(sampleProgramme, 1:7)
@test println(sampleProgramme[1, 3, 5]) === println(sampleProgramme, 1:7)

write("sampleProgramme.csv", sampleProgramme; log = false)
rm("sampleProgramme.csv")
@test true

write("sampleProgramme.csv", sampleProgramme, 3; log = false)
rm("sampleProgramme.csv")
@test true

write("sampleProgramme.csv", sampleProgramme, 1:2:6; log = false)
rm("sampleProgramme.csv")
@test true

write("sampleProgramme.csv", sampleProgramme; log = true)
keyArr, date, day1, Δdays, reps, wght, rpe = loadLogFile(sampleProgramme)
@test length(keyArr) ==
      length(date) ==
      length(day1) ==
      length(Δdays) ==
      length(reps) ==
      length(wght) ==
      length(rpe)
rm("Log_sampleProgramme.csv")

week = makeDays(SampleProgramme(), exerProg)
sampleProgramme = Programme(SampleProgramme(), "sampleProgrammeP", exerProg, week)
keyArr, date, day1, Δdays, reps, wght, rpe = loadLogFile(sampleProgramme)
tm, change = calcTrainingMaxLogs(sampleProgramme, keyArr, reps, wght)
@test length(tm) == 2

figs =
    plotData(sampleProgramme, keyArr, Δdays, tm; xlabel = "Days", ylabel = "Weight", lw = 3)
for (i, fig) in enumerate(figs)
    scatterData!(
        fig,
        sampleProgramme,
        keyArr[i],
        Δdays,
        tm;
        shape = :circle,
        markersize = 5,
        label = "",
    )
end
@test figs[1].n == 2
figs = plotData(
    sampleProgramme,
    keyArr,
    Δdays,
    tm;
    xlabel = "Days",
    ylabel = "Weight",
    lw = 3,
    label = "",
)
for (i, fig) in enumerate(figs)
    scatterData!(
        fig,
        sampleProgramme,
        keyArr[i],
        Δdays,
        tm;
        shape = :circle,
        markersize = 5,
    )
end
@test figs[1].n == 2

figs = scatterData(
    sampleProgramme,
    keyArr,
    Δdays,
    tm;
    xlabel = "Days",
    ylabel = "Weight",
    lw = 3,
)
for (i, fig) in enumerate(figs)
    plotData!(
        fig,
        sampleProgramme,
        keyArr[i],
        Δdays,
        tm;
        shape = :circle,
        markersize = 5,
        label = "",
    )
end
@test figs[1].n == 2
figs = scatterData(
    sampleProgramme,
    keyArr,
    Δdays,
    tm;
    xlabel = "Days",
    ylabel = "Weight",
    lw = 3,
    label = "",
)
for (i, fig) in enumerate(figs)
    plotData!(fig, sampleProgramme, keyArr[i], Δdays, tm; shape = :circle, markersize = 5)
end
@test figs[1].n == 2

tm, change = updateMaxes(sampleProgramme, keyArr, reps)
@test tm == [56.32850633965048, 77.22038362825882]
@test change == [-42.5, -52.5]

old1 = sampleProgramme.exerProg["SampleExercise"].exercise.trainingMax
old2 = sampleProgramme.exerProg["SampleExercise2"].exercise.trainingMax

updateMaxes!(sampleProgramme, keyArr, reps)
new1 = sampleProgramme.exerProg["SampleExercise"].exercise.trainingMax
new2 = sampleProgramme.exerProg["SampleExercise2"].exercise.trainingMax
@test new1 == old1 - 42.5
@test new2 == old2 - 52.5

tm, change = adjustMaxes("SampleExercise", exerProg, 13, weight = 100)
@test tm == 110.76583082526619
@test change == 52.5
adjustMaxes!("SampleExercise", exerProg, 13, weight = 100)
@test sampleProgramme.exerProg["SampleExercise"].exercise.trainingMax == new1 + 52.5
@test sampleProgramme.exerProg["SampleExercise2"].exercise.trainingMax == new2

old = sampleProgramme.exerProg["SampleExercise"].exercise.trainingMax
tm, change = adjustMaxes("SampleExercise", exerProg, 10, weight = 92.5)
new = tm
@test tm == 92.5
@test change == 2.5
adjustMaxes!("SampleExercise", exerProg, 10, weight = 92.5)
@test sampleProgramme.exerProg["SampleExercise"].exercise.trainingMax == old + change

old = sampleProgramme.exerProg["SampleExercise"].exercise.trainingMax
tm, change = adjustMaxes("SampleExercise", exerProg, 1, weight = 300)
new = tm
@test tm == 213.98592973630116
@test change == 100.0
adjustMaxes!("SampleExercise", exerProg, 1, weight = 300)
@test sampleProgramme.exerProg["SampleExercise"].exercise.trainingMax == old + change

max_stars_repo_name: dcelisgarza/Lifting


text:
# supertype for all logit model results

using StatsBase

abstract type LogitModel <: RegressionModel end

StatsBase.coefnames(r::LogitModel) = r.coefnames
StatsBase.coef(r::LogitModel) = r.coefs
StatsBase.islinear(r::LogitModel) = false
StatsBase.loglikelihood(r::LogitModel) = r.final_ll
StatsBase.vcov(r::LogitModel) = r.vcov
StatsBase.stderror(r::LogitModel) = sqrt.(diag(StatsBase.vcov(r)))


text:
function buildmodel!(op_model::PowerOperationModel, sys::PowerSystems.PowerSystem; args...)

    #TODO: Add check model spec vs data functions before trying to build

    netinjection = instantiate_network(op_model.transmission, sys)

    for category in op_model.generation
        constructdevice!(op_model.model, netinjection, category.device, category.formulation, op_model.transmission, sys; args...)
    end

    if op_model.demand != nothing
        for category in op_model.demand
            constructdevice!(op_model.model, netinjection, category.device, category.formulation, op_model.transmission, sys; args...)
        end
    end

    #=
    for category in op_model.storage
        op_model.model = constructdevice!(category.device, network_model, op_model.model, devices_netinjection, sys, category.constraints)
    end
    =#
    if op_model.services != nothing
        service_providers = Array{NamedTuple{(:device, :formulation),Tuple{DataType,DataType}}}([])
        [push!(service_providers,x) for x in vcat(op_model.generation,op_model.demand,op_model.storage) if x != nothing]
        for service in op_model.services
            op_model.model = constructservice!(op_model.model, service.service, service.formulation, service_providers, sys; args...)
        end
    end

    constructnetwork!(op_model.model, op_model.branches, netinjection, op_model.transmission, sys; args..., PTDF = op_model.ptdf)

    @objective(op_model.model, Min, op_model.model.obj_dict[:objective_function])

   return op_model

end

function build_sim_ts(ts_dict::Dict{String,Any}, steps, periods, resolution, date_from, lookahead_periods, lookahead_resolution ; args...)
    # exmaple of time series assembly
    # TODO: once we refactor PowerSystems, we can improve this process

    steps_dict = Dict()

    function _subset_ts(ts_dict,start,finish)
        return ts_dict[(ts_dict.DateTime .>= start) .& (ts_dict.DateTime .< finish),:]
    end

    for step in 1:steps
        step_stamp = date_from + ((resolution * periods) + (lookahead_periods * lookahead_resolution)) * (step-1)
        step_end = step_stamp + (resolution * periods) + (lookahead_periods * lookahead_resolution)
        steps_dict[step_stamp] = deepcopy(ts_dict)
        steps_dict[step_stamp]["load"] = _subset_ts(steps_dict[step_stamp]["load"],step_stamp,step_end)
        for cat in keys(steps_dict[step_stamp]["gen"])
            steps_dict[step_stamp]["gen"][cat] = _subset_ts(steps_dict[step_stamp]["gen"][cat],step_stamp,step_end)
        end
    end

    return steps_dict
end

function buildsimulation!(sys::PowerSystems.PowerSystem, op_model::PowerOperationModel, ts_dict::Dict{String,Any}; args...)

    name = :name in keys(args) ? args[:name] : "my_simulation"

    model = op_model

    resolution = :resolution in keys(args) ? args[:resolution] : PowerSystems.getresolution(sys.loads[1].scalingfactor)

    date_from = :date_from in keys(args) ? args[:date_from] : minimum(timestamp(sys.loads[1].scalingfactor))

    date_to = :date_to in keys(args) ? args[:date_to] : maximum(timestamp(sys.loads[1].scalingfactor))

    periods = :periods in keys(args) ? args[:periods] : (resolution < Hour(1) ? 1 : 24)

    steps = :steps in keys(args) ? args[:steps] : Int64(floor((length(sys.loads[1].scalingfactor)-1)/periods))

    if steps != (length(sys.loads[1].scalingfactor)-1)/periods
        @warn "Time series length and simulation definiton inconsistent, simulation may be truncated, simulating $steps steps."
    end

    lookahead_periods = :lookahead_periods in keys(args) ? args[:lookahead_periods] : 0

    lookahead_resolution = :lookahead_resolution in keys(args) ? args[:lookahead_resolution] : resolution

    @info "Simulation defined for $steps steps with $periods * $resolution periods per step (plus $lookahead_periods * $lookahead_resolution lookahead periods), from $date_from to $date_to"

    dynamic_analysis = false;

    timeseries = build_sim_ts(ts_dict, steps, periods, resolution, date_from, lookahead_periods, lookahead_resolution ; args...)

    PowerSimulationsModel(name,model, steps, periods, resolution, date_from, date_to,
            lookahead_periods, lookahead_resolution, dynamic_analysis, timeseries)

end

function buildsimulation!(sys::PowerSystems.PowerSystem, op_model::PowerOperationModel; args...)

    ts_dict = Dict{String,Any}()

    ts_dict["load"] = DataFrame(Dict([(l.name,values(l.scalingfactor)) for l in sys.loads]))
    ts_dict["load"][:DateTime] = TimeSeries.timestamp(sys.loads[1].scalingfactor)

    ts_dict["gen"] = Dict()

    if !isa(sys.generators.renewable,Nothing)
        # TODO: do a better job of classifying generators in the timeseries dict and reflect that here. For now, i'm just using PV as a placeholder
        ts_dict["gen"]["PV"] = DataFrame(Dict([(g.name,values(g.scalingfactor)) for g in sys.generators.renewable]))
        ts_dict["gen"]["PV"][:DateTime] = timestamp(sys.generators.renewable[1].scalingfactor)
        ts_dict["gen"]["WIND"] = DataFrame(Dict([(g.name,values(g.scalingfactor)) for g in sys.generators.renewable]))
        ts_dict["gen"]["WIND"][:DateTime] = timestamp(sys.generators.renewable[1].scalingfactor)
    end

    if !isa(sys.generators.hydro,Nothing)
        ts_dict["gen"]["Hydro"] = DataFrame(Dict([(g.name,values(g.scalingfactor)) for g in sys.generators.hydro]))
        ts_dict["gen"]["Hydro"][:DateTime] = timestamp(sys.generators.hydro[1].scalingfactor)
    end

    buildsimulation!(sys, op_model, ts_dict; args...)
end


text:
"""
Slider crank with ball
q = (θ1,θ2,x3,x4), u = (θ1dot,θ2dot,x3dot,x4dot),
g(q) = x4 - (x3 + w)
W = del g / del q = [0; 0; -1; 1]
"""

using LinearAlgebra
using NLsolve
using PyPlot
using TimeStepping

const l1 = 0.5
const l2 = 0.75
const m1 = 0.25
const m2 = 0.4
const m3 = 0.5
const m4 = 1.0
const J1 = 0.05
const J2 = 0.07
const g = 9.81
final_time = 2.0
q0 = [π-π/4; -π/4; 0.0; 1.0]
u0 = [0.0; 0.0; 0.0; -1.0]

function loop_fixed(q::AbstractArray)
    return [
        l1*cos(q0[1]) + l2*cos(q[1]) - q[2],
        l1*sin(q0[1]) + l2*sin(q[1])
    ]
end

sol = nlsolve(q->loop_fixed(q[1:2]), [-π/4, 1.0], autodiff=:forward)
q0[2:3] = sol.zero

function loop(q::AbstractArray)
    return [
        l1*cos(q[1]) + l2*cos(q[2]) - q[3],
        l1*sin(q[1]) + l2*sin(q[2])
    ]
end

function jac(q::AbstractArray)
    jac = zeros(Float32, 2, 4)
    jac[1,1] = -l1*sin(q[1])
    jac[1,2] = -l2*sin(q[2])
    jac[1,3] = -1
    jac[2,1] = l1*cos(q[1])
    jac[2,2] = l2*cos(q[2])
    return jac
end

function jacdot(q::AbstractArray, u::AbstractArray)
    jacdot = zeros(Float32, 2, 4)
    jacdot[1,1] = -l1*cos(q[1])*u[1]
    jacdot[1,2] = -l2*cos(q[2])*u[2]
    jacdot[2,1] = -l1*sin(q[1])*u[1]
    jacdot[2,2] = -l2*sin(q[2])*u[2]
    return jacdot
end

function gap(q::AbstractArray, u::AbstractArray)
    W = zeros(Float32, 4,1)
    W[3,1] = -1.0
    W[4,1] = 1.0
    return [q[4]-q[3]-0.1], W
end

function dynamics(q::AbstractArray, u::AbstractArray)
    M = zeros(Float32,4, 4)
    M[1,1] = 1/4*m1*l1*l1 + m2*l1*l1 + J1
    M[1,2] = M[2,1] = 1/2*m2*l1*l2*cos(q[1]-q[2])
    M[2,2] = 1/4*m2*l2*l2 + J2
    M[3,3] = m3
    M[4,4] = m4

    h = zeros(Float32, 4)
    h[1] = -m2*sin(q[1]-q[2])u[1]*u[2] + 1/2*m2*sin(q[1]-q[2])*u[2]*u[2] +
        (1/2*m1+m2)*g*l1*cos(q[1])
    # h[1] += 1.0*(q[1] - q0[1]) + 1.0*u[1]  # simple controller to keep the slider-crank fixed.
    h[2] = m2*sin(q[1]-q[2])u[1]*u[2] - 1/2*m2*sin(q[1]-q[2])*u[1]*u[1] + m2*g*l1*cos(q[2])
    h[3] = 0.0
    h[4] = 0.0
    return M, -h
end

# m = Moreau(gap, dynamics, loop, jac, jacdot, q0, u0)
sc_ball = Integrator(gap, dynamics, loop, jac, jacdot, convert.(Float32, q0),
    convert.(Float32, u0); Δt=1e-3, ε=0.5)
integrate(sc_ball, convert(Float32, final_time))

fig = figure(1, clear=true)
fig.add_subplot(1,1,1)
fig.axes[1].clear()
fig.axes[1].plot(sc_ball.t, getindex.(sc_ball.q,4), "k", linewidth=2, label="ball position")
fig.axes[1].plot(sc_ball.t, getindex.(sc_ball.q,3).+0.1, "b", linewidth=2, label="slider (bd) position")

# fig.axes[1].set_title("Ball motion", fontsize=18)
fig.axes[1].set_xlabel(L"$t$", fontsize=18)
fig.axes[1].set_ylabel(L"$x_3(t)$, $x_4(t)$", fontsize=18)
fig.axes[1].legend(loc="best", fontsize="x-large")

fig.axes[1].autoscale_view()
fig.tight_layout()

fig = figure(2, clear=true)
fig.add_subplot(1,1,1)
fig.axes[1].clear()
fig.axes[1].plot(sc_ball.t, getindex.(sc_ball.u,4), "k", linewidth=2, label="ball velocity")
fig.axes[1].plot(sc_ball.t, getindex.(sc_ball.u,3), "b", linewidth=2, label="slider velocity")

# fig.axes[1].set_title("Ball motion", fontsize=18)
fig.axes[1].set_xlabel(L"$t$", fontsize=18)
fig.axes[1].set_ylabel(L"$\dot{x}_3(t)$, $\dot{x}_4(t)$", fontsize=18)
fig.axes[1].legend(loc="best", fontsize="x-large")

fig.axes[1].autoscale_view()
fig.tight_layout()


text:
using Distributions
using PyPlot
using RDatasets
using MAT

include("./cluster_eval.jl")

srand(0) #seed

@everywhere function kpp_init(X::Array{Float64,2}, K::Int64)

    n, d = size(X)
    mu = zeros(K,d)
    dist = Inf * ones(n)
    Lambda = 1

    idx = round(Int64, rand(Uniform(1,n)))
    mu[1,:] = X[idx,:]
    for i = 2:K
        D = X - repmat(mu[i-1,:]',n,1)
        dist = min(dist, sum(D.*D,2))
        idx_cond = rand(Uniform(0,1)) .< cumsum(dist/sum(dist))
        idx_mu = find(idx_cond)[1]
        mu[i,:] = X[idx_mu,:]
        Lambda = maximum(dist)
    end
    Lambda
end

@everywhere function dpmeans_fit(X::Array{Float64, 2})

    #init params
    K = 1
    K_init = 9
    d = size(X,2)
    z = mod(randperm(size(X,1)),K)+1
    mu = rand(MvNormal(zeros(d), 10*eye(d)))
    sigma = 1
    nk = zeros(K)
    pik = ones(K)/K

    mu = mean(X,1)
    Lambda = kpp_init(X, K_init)
    println("lambda: ", Lambda)

    max_iter = 100
    obj = zeros(max_iter)
    em_time = zeros(max_iter)

    obj_tol = 1e-3
    n, d = size(X)

    for iter = 1:max_iter
        dist = zeros(n,K)

        #assignment step
        for kk = 1:K
            Xm = X - repmat(mu[kk,:]',n,1)
            dist[:,kk] = sum(Xm.*Xm,2)
        end

        #update labels
        dmin = minimum(dist,2)
        z = mapslices(indmin, dist, 2)
        #find(dist->dist==0, dist)
        idx = find(dmin->dmin > Lambda, dmin)
        if (size(idx)[1] > 0)
            K = K + 1
            z[idx] = K
            mu = cat(1, mu, mean(X[idx,:],1))
            Xm = X - repmat(mu[K,:]',n,1)
            dist = cat(2, dist, sum(Xm.*Xm,2))
        end

        #update step
        nk = zeros(K)
        for kk = 1:K
            nk[kk] = counts(z, kk)[1]
            idx2 = find(z->z==kk, z)
            if (length(idx2)>0)
                mu[kk,:] = mean(X[idx2,:],1)
            end
        end
        pik = nk / sum(nk)

        #compute objective
        for kk = 1:K
            idx = find(z->z==kk, z)
            obj[iter] = obj[iter] + sum(dist[idx,kk])
        end
        obj[iter] = obj[iter] + Lambda * K

        #check convergence
        if (iter > 1 && abs(obj[iter] - obj[iter-1]) < obj_tol * obj[iter])
            println("converged in ", iter, " iterations.")
            break
        end
        em_time[iter] = 0

    end

    z, obj, em_time
end

#generate data (X is NxD)
load_data = "synthetic"  #"synthetic", "iris", "wine", "soybean", "pima", "car"

if load_data == "synthetic"
    println("synthetic dataset")
    file = matread("./data/clusters.mat")
    X = transpose(file["data"])
    z_true = round(Int64, file["labels"])
elseif load_data == "iris"
    println("iris dataset")
    iris = dataset("datasets","iris")
    X = convert(Array, iris[:,1:4])
    species = convert(Array, iris[:,5])
    z_true = zeros(length(species))
    cnt = 1
    for sp in unique(species)
        idx = find(species->species == sp, species)
        z_true[idx] = cnt
        cnt = cnt + 1
    end
    z_true = round(Int64, z_true)
elseif load_data == "wine"
    println("wine dataset")
    file = matread("./data/wine.mat")
    X = 1.0 .* file["wine"]
    z_true = round(Int64, file["label"])
    z_true = squeeze(z_true, 2)
elseif load_data == "soybean"
    println("soybean dataset")
    file = matread("./data/soybean.mat")
    X = file["data"]
    z_true = round(Int64, file["label"])
    z_true = squeeze(z_true, 2)
elseif load_data == "pima"
    println("pima dataset")
    file = matread("./data/pima.mat")
    X = file["data"]
    z_true = round(Int64, file["label"])
    z_true = squeeze(z_true, 2)
elseif load_data == "car"
    println("car dataset")
    file = matread("./data/car.mat")
    X = file["data"]
    z_true = round(Int64, file["label"])
    z_true = squeeze(z_true, 2)
else
    println("choose a dataset")
    quit()
end

#DP means
tic()
labels, obj, em_time = dpmeans_fit(X)
ftime = toc()

#evaluate clustering
labels_ord = zeros(length(labels))
cnt = 1
for kk in unique(labels)
    cidx = find(labels->labels == kk, labels)
    labels_ord[cidx] = cnt
    cnt = cnt + 1
end
labels_ord = round(Int64, labels_ord)
println("K true: ", length(unique(z_true)))
println("K dpmeans: ", length(unique(labels_ord)))
cluster_eval(labels_ord, z_true)

#generate plots
figure()
scatter(X[:,1], X[:,2], c=z_true)
title("ground truth clustering")

figure()
scatter(X[:,1], X[:,2], c=labels_ord)
title("dp-means clustering")


text:
function createmesh!(data::Dict{Symbol, Any}, params::Parameters)
  el = 0
  seg = data[:segments]
  sgmedium = data[:media][data[:medium]] / 231.     # in [in^3]
  for i in 1:length(seg)
    mat = data[:materials][seg[i].material]
    radius = seg[i].od / 2.0
    dm2 = seg[i].od^2 - seg[i].id^2
    dm4 = seg[i].od^4 - seg[i].id^4
    dp2 = seg[i].od^2 + seg[i].id^2
    if typeof(seg[i]) in [Bit, Stabilizer, Bentsub, Rig]
      if (el+1) < data[:noofelements] && radius > params.nodes[el+1].stringradius
        params.nodes[el+1].fc = seg[i].fc
        params.nodes[el+1].stringradius = radius
      end
    else
      nels = floor(Int64, 1.0+data[:ratio]*data[:segments][i].length/data[:segments][i].od)
      dl = seg[i].length
      elength = 12.0 * dl / nels

      # Cross sectional area: csa = pi * d^2 / 4
      csa = pi * dm2 / 4.0
      sgstring = mat[:sg]
      eweight = elength * csa * (sgstring - sgmedium)

      # See http://en.wikipedia.org/wiki/Area_moment_of_inertia
      # Area moment of inertia: Pi * r^4 / 4 or Pi * d^4 / 64

      # See http://en.wikipedia.org/wiki/Polar_moment_of_inertia
      # Polar moment of inertia: Pi * r^4 / 2 or Pi * d^4 / 32

      ea = csa * mat[:ez]
      er = pi * mat[:er] * dm4 / 64
      gj = pi * mat[:g] * dm4 / 32

      # g is in [m / sec^2], 32.2 [ft/sec^2] or 386.4 [in/sec^2]
      # In US-units the assumption is that sgstring & sgmedium are forces [lbf],
      # not [lbm].

      #println([nels, elength, sgstring, sgmedium, eweight, ea, eiz, eir, gj])

      emass = elength * csa * sgstring
      mmass = elength * pi * seg[i].id^2 * sgmedium / 4.0

      # See Wikipedia, list of moments of inertia:
      #       http://en.wikipedia.org/wiki/List_of_moments_of_inertia
      # Iz = elementaxialmoment = elementmass * (od^2 + id^2) / 8
      # Ix = Iy = (1 / 12) * emass * ( 3 * (od^2 + id^2) / 4 + elength^2)
      #         = emass * elength^2 / 12 + elementaxialmoment / 2

      einertialmass = emass + mmass
      eaxialmoment = emass * dp2 / 8.0
      eradialmoment = emass * elength^2 / 12.0 + eaxialmoment / 2.0

      for j in 1:nels
      	el += 1

      	# params.nodes[1] is usually the bit or a 1st support for deflection testing
      	# params.nodes[1] = 0.0, add length of element.

        params.nodes[el+1].z = params.nodes[el].z + elength
        if typeof(seg[i]) in [Pipe, Collar]
          params.elements[el].ea = ea
          params.elements[el].er = er
          params.elements[el].gj = gj
        else
          # Adjust these values for Stabilizers. :od is not the whole story.
          params.elements[el].ea = 0.58 * ea
          params.elements[el].er = 0.68 * er
          params.elements[el].gj = 0.36 * gj
        end

				# Use largest stringRadius if 2 elements differ.

        if radius > params.nodes[el].stringradius
          params.nodes[el].fc = seg[i].fc
        end
				params.nodes[el].stringradius = max(radius, params.nodes[el].stringradius)
				params.nodes[el+1].fc = seg[i].fc
				params.nodes[el+1].stringradius = radius

				# Pattern always the same, asign 0.5 to this node and 0.5 to next node,
				# during next loop add 0.5 of next element.

				params.nodes[el].weight += eweight / 2.0
				params.nodes[el+1].weight = eweight / 2.0
				params.nodes[el].mass += emass / 2.0
				params.nodes[el+1].mass = emass / 2.0
				params.nodes[el].massinertia += einertialmass / 2.0
				params.nodes[el+1].massinertia = einertialmass / 2.0
				params.nodes[el].axialinertia += eaxialmoment / 2.0
				params.nodes[el+1].axialinertia = eaxialmoment / 2.0
				params.nodes[el].radialinertia += eradialmoment / 2.0
				params.nodes[el+1].radialinertia = eradialmoment / 2.0

				# Set element data if it is not a "(Exc)Stab(ilizer)" or "NBS" or "Bit" or "Rig"

				if typeof(seg[i]) in [Pipe, Collar]
					params.elements[el].fc = seg[i].fc
					params.elements[el].length = elength
					params.elements[el].radius = radius
					params.elements[el].od = seg[i].od
					params.elements[el].id = seg[i].id
					params.elements[el].mass = emass
				end
      end
    end
  end
  typeof(params.nodes)
end


text:
######### StanDataFrames program example  ###########

using CmdStan, StanDataFrames, Test, Statistics

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  "

  observeddata = [
    Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
  ]

  global stanmodel, rc, sim, cnames
  stanmodel = Stanmodel(num_samples=1200, thin=2, name="bernoulli",
    model=bernoullimodel);

  rc, sim, cnames = stan(stanmodel, observeddata, ProjDir, diagnostics=false,
    CmdStanDir=CMDSTAN_HOME);

  @test 0.1 <  mean(sim[:, 8, :]) < 0.5

  dfa = convert_a3d(sim, cnames, Val(:dataframe))

  @test 0.1 <  mean(dfa[1][:theta]) < 0.5

end # cd


text:
# -----------------------------------------------------------------------------
# Parametric type design
# -----------------------------------------------------------------------------

using Dates: Date

# Abstract type hierarchy for personal assets
abstract type Asset end
abstract type Investment <: Asset end
abstract type Equity <: Investment end

# Equity Instruments Types

# Define concrete types
struct Stock <: Equity
    symbol::String
    name::String
end

# Types of stock options
@enum CallPut Call Put

struct StockOption <: Equity
    symbol::String
    type::CallPut
    strike::Float64
    expiration::Date
end

# Trading Types

abstract type Trade end

# Types (direction) of the trade
@enum LongShort Long Short

struct SingleTrade{T <: Investment} <: Trade
    type::LongShort
    instrument::T
    quantity::Int
    price::Float64
end

stock = Stock("AAPL", "Apple Inc")
option = StockOption("AAPLC", Call, 200, Date(2019, 12, 20))
SingleTrade(Long, stock, 100, 188.0)
SingleTrade(Long, option, 100, 3.5)
#= REPL
julia> stock = Stock("AAPL", "Apple Inc")
Stock("AAPL", "Apple Inc")

julia> option = StockOption("AAPLC", Call, 200, Date(2019, 12, 20))
StockOption("AAPLC", Call::CallPut = 0, 200.0, 2019-12-20)

julia> SingleTrade(Long, stock, 100, 188.0)
SingleTrade{Stock}(Long::LongShort = 0, Stock("AAPL", "Apple Inc"), 100, 188.0)

julia> SingleTrade(Long, option, 100, 3.5)
SingleTrade{StockOption}(Long::LongShort = 0, StockOption("AAPLC", Call::CallPut = 0, 200.0, 2019-12-20), 100, 3.5)
=#

SingleTrade{Stock} <: SingleTrade
SingleTrade{StockOption} <: SingleTrade
#=
julia> SingleTrade{Stock} <: SingleTrade
true

julia> SingleTrade{StockOption} <: SingleTrade
true
=#

# -----------------------------------------------------------------------------
# Defining parametric methods
# -----------------------------------------------------------------------------

# Return + or - sign for the direction of trade
function sign(t::SingleTrade{T}) where {T}
    return t.type == Long ? 1 : -1
end

# Calculate payment amount for the trade
function payment(t::SingleTrade{T}) where {T}
    return sign(t) * t.quantity * t.price
end

SingleTrade(Long, stock, 100, 188.0) |> payment
SingleTrade(Long, option, 1, 3.50) |> payment
#= REPL
julia> SingleTrade(Long, stock, 100, 188.0) |> payment
18800.0

julia> SingleTrade(Long, option, 1, 3.50) |> payment
3.5
=#

# Calculate payment amount for option trades (100 shares per contract)
function payment(t::SingleTrade{StockOption})
    return sign(t) * t.quantity * 100 * t.price
end

SingleTrade(Long, stock, 100, 188.0) |> payment
SingleTrade(Long, option, 1, 3.50) |> payment
#= REPL
julia> SingleTrade(Long, stock, 100, 188.0) |> payment
18800.0

julia> SingleTrade(Long, option, 1, 3.50) |> payment
350.0
=#

# -----------------------------------------------------------------------------
# Using multiple parametric type arguments
# -----------------------------------------------------------------------------

struct PairTrade{T <: Investment, S <: Investment} <: Trade
    leg1::SingleTrade{T}
    leg2::SingleTrade{S}
end

payment(t::PairTrade) = payment(t.leg1) + payment(t.leg2)

stock
option
pt = PairTrade(SingleTrade(Long, stock, 100, 188.0), SingleTrade(Short, option, 1, 3.5));
payment(pt)
#= REPL
julia> stock
Stock("AAPL", "Apple Inc")

julia> option
StockOption("AAPLC", Call::CallPut = 0, 200.0, 2019-12-20)

julia> pt = PairTrade(SingleTrade(Long, stock, 100, 188.0), SingleTrade(Short, option, 1, 3.5));

julia> payment(pt)
18450.0
=#


text:
mutable struct LazyDictFromIterator{K, V}
    iter::Any
    state::Union{Nothing, Tuple}
    getter::Function
    items::Dict{K, V}
end

"""
LazyDictFromIterator creates a dictionary from an iterator, but only increments the iterator
and adds items to the dictionary as it needs them. In the worst case it is identical to
creating a dictionary by iterating over the entire list.
Each V should have a K member.

# Arguments
- `K`: type of the dictionary keys
- `V`: type of the dictionary values
- `iter`: any object implementing the Iterator interface
- `getter::Function`: method to call on V to get its K

"""
function LazyDictFromIterator(
                              ::Type{K},
                              ::Type{V},
                              iter::Any,
                              getter::Function,
                             ) where {K, V}
    return LazyDictFromIterator(iter, nothing, getter, Dict{K, V}())
end

"""
Returns the item mapped to key. If the key is already stored then it will be returned
with a dictionary lookup. If it has not been stored then iterate over the list until it is
found.

Returns nothing if key is not found.
"""
function Base.get(container::LazyDictFromIterator, key::K) where K
    if haskey(container.items, key)
        return container.items[key]
    end

    # Exit loop when item is found or throw an exception.
    while true
        if isnothing(container.state)
            result = Base.iterate(container.iter)
        else
            result = Base.iterate(container.iter, container.state)
        end

        if isnothing(result)
            @debug "Key not found" key
            return nothing
        end

        item = result[1]
        container.state = result[2]
        item_key = container.getter(item)

        # Store this item for future lookups.
        container.items[item_key] = item
        if key == item_key
            return item
        end
    end

    @assert false
end

"""Reset the iterator for cases where underlying arrays have changed."""
function reset_iterator(container::LazyDictFromIterator)
    @debug "reset_iterator"
    container.state = nothing
end

"""Replace the iterator, maintaining the cached dict."""
function replace_iterator(container::LazyDictFromIterator, iter)
    @debug "replace_iterator"
    container.state = nothing
    container.iter = iter
end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03rd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 10
    LDA = NMAX
    LDX = NMAX
    LDWORK = 3*NMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    X = Array{Float64,2}(undef, LDX,NMAX)
    BLSIZE = Array{BlasInt,1}(undef, NMAX)
    WR = Array{Float64,1}(undef, NMAX)
    WI = Array{Float64,1}(undef, NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    BWORK = Array{Bool,1}(undef, NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    PMAX = parse(Float64, replace(vs[2],'D'=>'E'))
    TOL = parse(Float64, replace(vs[3],'D'=>'E'))
    JOBX =  vs[4][1]
    SORT =  vs[5][1]
    if ( N<0 || N>NMAX )
        @error "Invalid N=$N"
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    #FOREIGN.dgees!( 'Vectors', 'Not sorted', SELECT, N, A, LDA, SDIM, WR, WI, X, LDX, DWORK, LDWORK, BWORK, INFO )
    Aschur = schur(A[1:N,1:N])
    A[1:N,1:N] .= Aschur.T
    X[1:N,1:N] .= Aschur.Z

# interp call 2

    NBLCKS, INFO = SLICOT.mb03rd!(JOBX, SORT, N, PMAX, A, X, BLSIZE, WR, WI, TOL)
    println(io, "NBLCKS = $NBLCKS")
    @test INFO == 0
    INFO == 0 || return

    if ( INFO!=0 )
        @warn "mb03rd returns INFO=$INFO"
    end
# interp output 1
    println(io, "BLSIZE:")
    _nr = NBLCKS
    show(io, "text/plain", BLSIZE[1:_nr])
    println(io)

# interp output 2
    println(io, "A:")
    _nc = N
    _nr = N
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp output 3
    println(io, "X:")
    _nc = N
    _nr = N
    show(io, "text/plain", X[1:_nr,1:_nc])
    println(io)

    close(f)
end # run_mb03rd()


text:
module RNLA
    using LinearAlgebra

    export rsvd, rpca, rrpca, rid, rcur

    include("rsvd.jl")
    include("randinit.jl")
    include("rqb.jl")
    include("rpca.jl")
    include("sub_iterations.jl")
    # include("rrpca.jl")
    include("rid.jl")
    include("rcur.jl")
    # include("rnmf.jl")
    # include("rdmd.jl")

end # module


text:
module Mergepdfs

#https://stackoverflow.com/questions/64520684/merge-pdf-files-using-python-pypdf2

using Conda
using PyCall

export mergepdfs

export pyModulePyPDF2
pyModulePyPDF2 = PyCall.PyNULL()

#https://anaconda.org/conda-forge/pypdf2

try
    #this line may trigger installation
    copy!(pyModulePyPDF2, PyCall.pyimport("PyPDF2"))
catch
    try
    copy!(pyModulePyPDF2, PyCall.pyimport_conda("PdfFileMerger","PyPDF2","conda-forge"))
    catch
        using Conda
        Conda.add("PyPDF2")
    end
end

function __init__()
	copy!(pyModulePyPDF2, PyCall.pyimport("PyPDF2"))

    py"""
from PyPDF2 import PdfFileMerger
from PyPDF2 import PdfFileReader
from PyPDF2 import PdfFileWriter

from hashlib import md5

from PyPDF2 import PdfFileReader, PdfFileWriter
from PyPDF2.generic import NameObject, DictionaryObject, ArrayObject, NumberObject, ByteStringObject
from PyPDF2.pdf import _alg33, _alg34, _alg35
from PyPDF2.utils import b_

#user_pwd (str) – The user password, which allows for opening and reading the PDF file with the restrictions provided.
#owner_pwd (str) – The owner password, which allows for opening the PDF files without any restrictions. By default, the owner password is the same as the user password.

def encrypt_modification(writer_obj: PdfFileWriter, user_pwd, owner_pwd=None, use_128bit=True):
    import time, random
    if owner_pwd == None:
        owner_pwd = user_pwd
    if use_128bit:
        V = 2
        rev = 3
        keylen = int(128 / 8)
    else:
        V = 1
        rev = 2
        keylen = int(40 / 8)
    # permit copy and printing only:
    P = -44
    O = ByteStringObject(_alg33(owner_pwd, user_pwd, rev, keylen))
    ID_1 = ByteStringObject(md5(b_(repr(time.time()))).digest())
    ID_2 = ByteStringObject(md5(b_(repr(random.random()))).digest())
    writer_obj._ID = ArrayObject((ID_1, ID_2))
    if rev == 2:
        U, key = _alg34(user_pwd, O, P, ID_1)
    else:
        assert rev == 3
        U, key = _alg35(user_pwd, rev, keylen, O, P, ID_1, False)
    encrypt = DictionaryObject()
    encrypt[NameObject("/Filter")] = NameObject("/Standard")
    encrypt[NameObject("/V")] = NumberObject(V)
    if V == 2:
        encrypt[NameObject("/Length")] = NumberObject(keylen * 8)
    encrypt[NameObject("/R")] = NumberObject(rev)
    encrypt[NameObject("/O")] = ByteStringObject(O)
    encrypt[NameObject("/U")] = ByteStringObject(U)
    encrypt[NameObject("/P")] = NumberObject(P)
    writer_obj._encrypt = writer_obj._addObject(encrypt)
    writer_obj._encrypt_key = key

#unmeta = PdfFileReader('my_pdf.pdf')

#writer = PdfFileWriter()
#writer.appendPagesFromReader(unmeta)
#encrypt(writer, '1', '123')

def read_and_encrypt(pdf_inputfile,pdf_outfile,secure_file=False,usrpwd=None,ownrpwd=None):
    #print('1')
    #print(pdf_inputfile)
    pdf = PdfFileReader(pdf_inputfile)
    #print('2')
    pdf_writer = PdfFileWriter()
    #print('3')
    for page in range(pdf.getNumPages()):
        pdfpage = pdf.getPage(page)
        pdf_writer.addPage(pdfpage)

    #print('4')

    with open(pdf_outfile, 'wb') as out:
        if secure_file:
            encrypt_modification(pdf_writer, usrpwd, ownrpwd, use_128bit=True)
        pdf_writer.write(out)

    #print('5')

    return 0
"""

end

function mergepdfs(list_of_pdf_files,resfi)
    @assert all(isfile.(list_of_pdf_files))
    endings = map(x->lowercase.(splitext(x)[2]),list_of_pdf_files)
    @assert all(".pdf".==endings)

    @assert isdir(splitdir(resfi)[1])

    merger = pyModulePyPDF2.PdfFileMerger()
    for fi in list_of_pdf_files
        merger.append(fi)
    end

    isfile(resfi) && rm(resfi)
    merger.write(resfi)
    merger.close()

    return nothing
end

export read_and_encrypt_jl
function read_and_encrypt_jl(pdf_inputfile,pdf_outfile,secure_file=false,usrpwd="",ownrpwd="")
    return py"read_and_encrypt"(pdf_inputfile,pdf_outfile,secure_file=secure_file,usrpwd=usrpwd,ownrpwd=ownrpwd)
end

end #module


text:
using ModelicaScriptingTools: withOMC, testmodel, installAndLoad
using Base.Filesystem
using Test
using OMJulia

projectroot = dirname(@__DIR__)

outdir = joinpath(projectroot, "out")
modeldir = projectroot
# use reference data from SHM project if available
refdir = joinpath(projectroot, "../shm-ref")

if !ispath(outdir)
    mkdir(outdir)
end

withOMC(outdir, modeldir) do omc
    installAndLoad(omc, "Modelica"; version="3.2.3")
    @testset "Simulate examples" begin
        @testset "UnidirectionalConductionExample" begin
            testmodel(omc, "SHMConduction.Examples.ModularExample"; refdir=refdir)
        end
        @testset "PVCExample" begin
            testmodel(omc, "SHMConduction.Examples.PVCExample"; refdir=refdir)
            # we cannot change with_sinus using simflags and -override
            # because it is a structural parameter => create wrapper model
            OMJulia.sendExpression(omc, "model PVCNoSinus extends SHMConduction.Examples.PVCExample(with_sinus=false); end PVCNoSinus;")
            es = OMJulia.sendExpression(omc, "getErrorString()")
            @test isempty(es)
            testmodel(omc, "PVCNoSinus";  override=Dict(
                "stopTime"=>55, "numberOfIntervals"=>27500,
                "variableFilter"=>"con\\.d_interbeat|trigger"
            ), refdir=refdir)
        end
    end
end

max_stars_repo_name: CSchoel/shm-conduction

max_issues_repo_name: CSchoel/shm-conduction

max_forks_repo_name: CSchoel/shm-conduction


text:
using PEPS
using Test

@testset "PEPS.jl" begin
    # Write your own tests here.
end


text: println("Hello Julia");

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))
using Documenter, ClimaCore

format = Documenter.HTML(
    prettyurls = !isempty(get(ENV, "CI", "")),
    collapselevel = 1,
)

makedocs(
    sitename = "ClimaCore.jl",
    strict = false,
    format = format,
    checkdocs = :exports,
    clean = true,
    doctest = true,
    modules = [ClimaCore],
    pages = Any[
        "Home" => "index.md",
        "API" => "api.md",
        "Operators" => "operators.md",
    ],
)

deploydocs(
    repo = "github.com/CliMA/ClimaCore.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)


text:
"""
    struct Instance{J,C,T}
An instance to the GCSPET problem is defined by:
- set of jobs of type `J`
- set of crane starting positions of type `T`
- crane speed of type `T`
- safety distance of type `T`
An instance can be generated through parameters:
- number of jobs
- number of cranes
- load factor
- crane speed      (default = 1)
- safety distance  (default = 1)
"""
struct Instance
    name::String
    Ω::Vector{Job}
    Q::Vector{Crane}
end

name(i::Instance) = i.name
jobs(i::Instance) = i.Ω
cranes(i::Instance) = i.Q
crane_starting_pos(i::Instance) = starting_position.(cranes(i))

"""
    parameters(i::Instance)

Returns a `Dict` containing the names and values of all instance-parameters.
"""
parameters(i::Instance) = Dict(["njobs" => njobs(i), "ncranes" => ncranes(i), "load" => loadfactor(i)])

njobs(i::Instance) = length(jobs(i))
ncranes(i::Instance) = length(cranes(i))
loadfactor(i::Instance) = count(istruck,jobs(i))

xmin(c, safety) = (c-1)*(safety+1)+1
xmax(c, n, njobs, safety) = njobs - (n-c)*(safety+1)
zone(c, ncranes, njobs, safety) = xmin(c, safety):xmax(c, ncranes, njobs, safety)

function Base.show(io::IO, ::MIME"text/plain", i::Instance)
    println(io, "Instance of GCSPET: $(name(i))")
    println(io, "Jobs")
    for j in jobs(i)
        println(io, j)
    end
    for c in cranes(i)
        println(io, c)
    end
end


text:
@testset "872.leaf-similar-trees.jl" begin
    t1 = convert(TreeNode{Int}, [3, 5, 1, 6, 2, 9, 8, nothing, nothing, 7, 4])
    t2 = convert(
        TreeNode{Int},
        [3, 5, 1, 6, 7, 4, 2, nothing, nothing, nothing, nothing, nothing, nothing, 9, 8],
    )
    @test leaf_similar(t1, t2) == true
    @test leaf_similar(TreeNode(1), TreeNode(1)) == true
    @test leaf_similar(TreeNode(1), TreeNode(2)) == false
    @test leaf_similar(TreeNode(1, TreeNode(2)), TreeNode(2, TreeNode(2))) == true
    @test leaf_similar(
        TreeNode(1, TreeNode(2), TreeNode(3)), TreeNode(2, TreeNode(3), TreeNode(2))
    ) == false
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

function test_linearsolver()
    el1 = Quad4([1, 2, 3, 4])
    el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
    el1["temperature thermal conductivity"] = 6.0
    el1["density"] = 36.0

    el2 = Seg2([1, 2])
    el2["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]]
    el2["temperature flux"] = (
        (0.0 => 0.0),
        (1.0 => 600.0)
        )

    field_problem = HeatProblem()
    push!(field_problem, el1)
    push!(field_problem, el2)

    el3 = Seg2([3, 4])
    el3["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
    el3["temperature"] = 0.0

    boundary_problem = DirichletProblem("temperature", 1)

    push!(boundary_problem, el3)

    # Create a solver for a set of problems
    solver = LinearSolver()
    push!(solver, field_problem)
    push!(solver, boundary_problem)

    # Solve problem at time t=1.0 and update fields
    solver(1.0)

    # Postprocess.
    # Interpolate temperature field along boundary of Γ₁ at time t=1.0
    xi = [0.0, -1.0]
    X = el2("geometry", xi, 1.0)
    T = el2("temperature", xi, 1.0)
    info("Temperature at point X = $X is T = $T")
    @test isapprox(T, 100.0)
end

function test_solvers()
    K =  [
    440.0   150.0  -260.0   -30.0    40.0    30.0  -220.0  -150.0
    150.0   440.0    30.0    40.0   -30.0  -260.0  -150.0  -220.0
   -260.0    30.0   440.0  -150.0  -220.0   150.0    40.0   -30.0
    -30.0    40.0  -150.0   440.0   150.0  -220.0    30.0  -260.0
     40.0   -30.0  -220.0   150.0   440.0  -150.0  -260.0    30.0
     30.0  -260.0   150.0  -220.0  -150.0   440.0   -30.0    40.0
   -220.0  -150.0    40.0    30.0  -260.0   -30.0   440.0   150.0
   -150.0  -220.0   -30.0  -260.0    30.0    40.0   150.0   440.0]

    C = 1/3*[
    1.0  0.0  0.0  0.0  0.5  0.0  0.0  0.0
    0.0  1.0  0.0  0.0  0.0  0.5  0.0  0.0
    0.0  0.0  1.0  0.0  0.0  0.0  0.5  0.0
    0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.5
    0.5  0.0  0.0  0.0  1.0  0.0  0.0  0.0
    0.0  0.5  0.0  0.0  0.0  1.0  0.0  0.0
    0.0  0.0  0.5  0.0  0.0  0.0  1.0  0.0
    0.0  0.0  0.0  0.5  0.0  0.0  0.0  1.0]

    f = zeros(8)

    g = 1/100 * [-5.0, 5.0, 10.0, -10.0, -5.0, 5.0, 10.0, -10.0]

    K = sparse(K)
    C = sparse(C)
    f = sparse(f)
    g = sparse(g)

    expected = [-0.1, 0.1, 0.2, -0.2, -0.1, 0.1, 0.2, -0.2]
    u1, la1 = solve(K, f, C, g, Val{:UMFPACK})
    @test isapprox(u1, expected)
    u2, la2 = solve(K, f, C, g, Val{:CHOLMOD})
    @test isapprox(u2, expected)
    # FIXME: how to dynamically include packages only if they are installed?
    #include(Pkg.dir("JuliaFEM"*"/src/petsc.jl"))
    u3, la3 = solve(K, f, C, g, Val{:PETSc_GMRES})
    @test isapprox(u3, expected)
end


text:
# ============================================================
# Structs of Initial & Boundary Conditions
# ============================================================

"""
    mutable struct IB{T} <: AbstractCondition
        fw::Function
        bc::T
    end

Initial & boundary conditions with no distribution function

"""
mutable struct IB{T<:Union{Function,AbstractArray}} <: AbstractCondition
    fw::Function
    bc::T
end

function IB(fw, gas::Scalar)
    γ = gas.a

    bc = function (args...)
        w = fw(args...)
        return ifelse(γ == 0, conserve_prim(w), conserve_prim(w, γ))
    end

    return IB{typeof(bc)}(fw, bc)
end

function IB(fw, gas::AbstractGas)
    bc = function (args...)
        w = fw(args...)
        prim = begin
            if ndims(w) == 1
                conserve_prim(w, gas.γ)
            else
                mixture_conserve_prim(w, gas.γ)
            end
        end

        return prim
    end

    return IB{typeof(bc)}(fw, bc)
end

"""
    mutable struct IB1F{T} <: AbstractCondition
        fw::Function
        ff::Function
        bc::T
    end

Initial & boundary condition with 1 distribution function

"""
mutable struct IB1F{T} <: AbstractCondition
    fw::Function
    ff::Function
    bc::T
end

function IB1F(fw, vs::AbstractVelocitySpace, gas::AbstractProperty)
    bc = function (args...)
        w = fw(args...)
        prim = begin
            if ndims(w) == 1
                conserve_prim(w, gas.γ)
            else
                mixture_conserve_prim(w, gas.γ)
            end
        end

        return prim
    end

    ff = function (args...)
        prim = bc(args...)

        M = begin
            if !isdefined(vs, :v)
                if ndims(prim) == 1
                    maxwellian(vs.u, prim)
                else
                    mixture_maxwellian(vs.u, prim)
                end
            elseif !isdefined(vs, :w)
                if ndims(prim) == 1
                    maxwellian(vs.u, vs.v, prim)
                else
                    mixture_maxwellian(vs.u, vs.v, prim)
                end
            else
                if ndims(prim) == 1
                    maxwellian(vs.u, vs.v, vs.w, prim)
                else
                    mixture_maxwellian(vs.u, vs.v, vs.w, prim)
                end
            end
        end

        return M
    end

    return IB1F{typeof(bc)}(fw, ff, bc)
end

"""
    mutable struct IB2F{T} <: AbstractCondition
        fw::Function
        ff::Function
        bc::T
    end

Initial & boundary condition with 2 distribution functions

"""
mutable struct IB2F{T} <: AbstractCondition
    fw::Function
    ff::Function
    bc::T
end

function IB2F(fw, vs::AbstractVelocitySpace, gas::AbstractProperty)
    bc = function (args...)
        w = fw(args...)
        prim = begin
            if ndims(w) == 1
                conserve_prim(w, gas.γ)
            else
                mixture_conserve_prim(w, gas.γ)
            end
        end

        return prim
    end

    ff = function (args...)
        prim = bc(args...)

        if !isdefined(vs, :v)
            if ndims(prim) == 1
                H = maxwellian(vs.u, prim)
                B = H .* gas.K / 2 / prim[end]
            else
                H = mixture_maxwellian(vs.u, prim)
                B = zero(H)
                for j in axes(B, 2)
                    B[:, j] = H[:, j] * gas.K / (2.0 * prim[end, j])
                end
            end
        elseif !isdefined(vs, :w)
            if ndims(prim) == 1
                H = maxwellian(vs.u, vs.v, prim)
                B = H .* gas.K / 2 / prim[end]
            else
                H = mixture_maxwellian(vs.u, vs.v, prim)
                B = zero(H)
                for j in axes(B, 3)
                    B[:, :, j] = H[:, :, j] * gas.K / (2.0 * prim[end, j])
                end
            end
        else
            if ndims(prim) == 1
                H = maxwellian(vs.u, vs.v, vs.w, prim)
                B = H .* gas.K / 2 / prim[end]
            else
                H = mixture_maxwellian(vs.u, vs.v, vs.w, prim)
                B = zero(H)
                for j in axes(B, 4)
                    B[:, :, :, j] = H[:, :, :, j] * gas.K / (2.0 * prim[end, j])
                end
            end
        end

        return H, B
    end

    return IB2F{typeof(bc)}(fw, ff, bc)
end

"""
    mutable struct IB3F{T} <: AbstractCondition
        fw::Function
        ff::Function
        fE::Function
        fB::Function
        fL::Function
        bc::T
    end

Initial & boundary condition with 3 distribution functions

"""
mutable struct IB3F{T} <: AbstractCondition
    fw::Function
    ff::Function
    fE::Function
    fB::Function
    fL::Function
    bc::T
end

"""
    mutable struct IB4F{T} <: AbstractCondition
        fw::Function
        ff::Function
        fE::Function
        fB::Function
        fL::Function
        bc::T
    end

Initial & boundary condition with 4 distribution functions

"""
mutable struct IB4F{T} <: AbstractCondition
    fw::Function
    ff::Function
    fE::Function
    fB::Function
    fL::Function
    bc::T
end


text:
module Labels

# package code goes here
export Label, varlab, vallab, lblname, defined

struct Label
    var::Dict{Symbol,String}
    val::Dict{Symbol,Dict}
    lblname::Dict{Symbol,Symbol}

    Label(var,val,lblname) = new(var,val,lblname)
    Label() = Label(Dict(),Dict(),Dict())
end

"""
    varlab(l::Label,v::Symbol)

Produces the variable description associated with the `v` variable in the `l` Labels.
"""
function varlab(l::Label,v::Symbol)
    return haskey(l.var,v) ? l.var[v] : string(v)
end

"""
    vallab(l::Label,v::Symbol,val)

Returns the value label associated with `val` value for the `v` variable in the `l` Labels.
"""
function vallab(l::Label,v::Symbol,val)
    lname = lblname(l,v)

    if lname == nothing
        return string(val)
    end

    return haskey(l.val,lname) && haskey(l.val[lname],val) ? l.val[lname][val] : string(val)
end

"""
    lblname(l::Label,v::Symbol)

Finds label name associated with `v` variable in the `l` Labels. If `v` does not have
associated label name, `nothing` is returned.
"""
function lblname(l::Label,v::Symbol)
    if haskey(l.lblname,v)
        return l.lblname[v]
    end
    # nothing is returned here
end
#
# function defined(l::Label,lbl::Symbol)
#     return in(lbl,collect(values(l.lblname)))
# end
#
# function defined(l::Label,v::Symbol)
#     return v in collect(keys(l.var)) || v in collect(keys(l.lblname))
# end

function lblname(l::Label,v::Symbol)
    if haskey(l.lblname,v)
        return l.lblname[v]
    end
end

end # module


text:
using Plots
using Random
Random.seed!(1)

include_ddp()

# Model
include_model("double_integrator")

function f(model::DoubleIntegratorContinuous, x, u, w)
    [x[2]; (1.0 + w[1]) * u[1]]
end

function fd(model::DoubleIntegratorContinuous{Midpoint, FixedTime}, x, u, w, h, t)
	x + h * f(model, x + 0.5 * h * f(model, x, u, w), u, w)
end

model = DoubleIntegratorContinuous{Midpoint, FixedTime}(2, 1, 1)
n = model.n
m = model.m

# Time
T = 11
h = 0.1

# reference
x_ref = [[0.0; 0.0] for t = 1:T]
xT = [0.0; 0.0]
ū_init = [1.0 * randn(model.m) for t = 1:T-1]
u_ref = [zeros(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Objective
Q = [(t < T ? h : 1.0) * (t < T ?
	 Diagonal([1.0; 1.0])
		: Diagonal([1.0; 1.0])) for t = 1:T]
q = [-2.0 * Q[t] * x_ref[t] for t = 1:T]
R = h * [Diagonal(1.0 * ones(model.m)) for t = 1:T-1]
r = [zeros(model.m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
ul = [-5.0]
uu = [5.0]
p = [t < T ? 2 * m : n for t = 1:T]
info_t = Dict(:ul => ul, :uu => uu, :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		c .= [ul - u; u - uu]
	else
		c .= x - cons.con[T].info[:xT]
	end
end

# generate data
N = 1
x_data = zeros(model.n, 0)
u_data = zeros(model.m, 0)

# offset = [[0.0; 0.0], [-2.0; 0.0]]
for i = 1:N
	# Initial conditions, controls, disturbances
	x1 = [1.0; 0.0] + 0.0 * randn(model.n)
	# x1 = [1.0; 0.0] + offset[i]

	# Rollout
	x̄ = rollout(model, x1, copy(ū_init), w, h, T)

	prob = problem_data(model, obj, con_set, copy(x̄), copy(ū_init), w, h, T,
		analytical_dynamics_derivatives = false)

	# Solve
	@time constrained_ddp_solve!(prob,
		linesearch = :armijo,
		max_iter = 1000,
		max_al_iter = 10,
		ρ_init = 1.0,
		ρ_scale = 10.0,
		con_tol = 1.0e-5)

	x, u = current_trajectory(prob)
	x̄, ū = nominal_trajectory(prob)

	x_data = hcat(x_data, hcat(x[1:end-1]...))
	u_data = hcat(u_data, hcat(u...))
end

# Visualize
# plot(hcat([[0.0; 0.0] for t = 1:T]...)',
#     width = 2.0, color = :black, label = "")
# plt = plot(hcat(x...)',
# 	width = 2.0,
# 	color = [:cyan :orange],
# 	label = ["x" "ẋ"],
# 	xlabel = "time step",
# 	ylabel = "state")

# policy
pass_through(x) = x
l = [model.n, 2, 2, model.m] # input -> output layers
activation = [:tanh, :tanh, :pass_through]
np = sum([l[i] * l[i+1] for i = 1:length(l)-1]) + sum(l[2:end])
reg_θ = 1.0e-5

function policy(x, θ)
	W = []
	B = []
	shift = 0

	p = x

	for i = 1:length(l)-1
		s = l[i] * l[i+1]
		w = θ[shift .+ (1:s)]
		shift += s

		b = θ[shift .+ (1:l[i+1])]
		shift += l[i+1]
		# p = reshape(w, l[i+1], l[i]) * p + b

		p = eval(activation[i]).(reshape(w, l[i+1], l[i]) * p + b)
	end

	return p
end

policy(ones(model.n), randn(np))

@variables x_sym[1:model.n], u_sym[1:model.m], θ_sym[1:np]

p_sym = policy(x_sym, θ_sym)
pθ_sym = Symbolics.jacobian(p_sym, θ_sym)

p_func = eval(Symbolics.build_function(p_sym, x_sym, θ_sym)[1])
pθ_func = eval(Symbolics.build_function(pθ_sym, x_sym, θ_sym)[1])

p_func(ones(model.n), randn(np))
pθ_func(ones(model.n), randn(np))

# loss
function err(θ, x, u)
	r = u - policy(x, θ)

	return 0.5 * norm(r)^2.0
end

e_sym = err(θ_sym, x_sym, u_sym)
eθ_sym = Symbolics.gradient(e_sym, θ_sym)
eθ²_sym = Symbolics.sparsehessian(e_sym, θ_sym)

eθ_func = eval(Symbolics.build_function(eθ_sym, θ_sym, x_sym, u_sym)[1])
eθ²_func = eval(Symbolics.build_function(eθ²_sym, θ_sym, x_sym, u_sym)[1])

eθ_func(ones(np), ones(model.n), ones(model.m))
eθ²_func(ones(np), ones(model.n), ones(model.m))

function mse(θ, x, u)
	N = size(x, 2)
	J = 0.0

	for t = 1:N
		J += 0.5 * norm(u[:, t] - p_func(x[:, t], θ))^2.0
	end

	return J / N + reg_θ * transpose(θ) * θ
end

mse(randn(np), x_data, u_data)

p_func(x_data[:, 1], randn(np)) - u_data[:, 1]

function mse_grad(θ, x, u)
	N = size(x, 2)
	grad = zero(θ)

	for t = 1:N
		grad .+= -1.0 * transpose(pθ_func(x[:, t], θ)) * (u[:, t] - p_func(x[:, t], θ))
		# grad .+= eθ_func(θ, x[:, t], u[:, t])
	end

	return grad ./ N + reg_θ * 2.0 * θ
end

norm(ForwardDiff.gradient(a -> mse(a, x_data, u_data), ones(np)) - mse_grad(ones(np), x_data, u_data))

function mse_hess(θ, x, u)
	N = size(x, 2)
	np = length(θ)
	hess = zeros(np, np)

	for t = 1:N
		hess .+= eθ²_func(θ, x[:, t], u[:, t])
	end

	return hess ./ N .+ reg_θ * 2.0
end

mse_hess(randn(np), x_data, u_data)

function gd(θ0, x_data, u_data; iters=1000, verbose=true, α=0.001)
    θ = copy(θ0)

	grad = zero(θ)

    for i = 1:iters
        grad = mse_grad(θ, x_data, u_data)
        θ -= α * grad

		if verbose && (i % 500 == 0)
			println("iter: $i")
			println("obj: ", mse(θ, x_data, u_data))
			println("grad. norm: ", norm(grad))
		end
    end

    return θ
end

function gd_ls(θ0, x_data, u_data; grad_tol = 1.0e-5, iters=100, verbose=true, )
    θ = copy(θ0)
	obj = mse(θ, x_data, u_data)
	grad = mse_grad(θ, x_data, u_data)

	iter = 1
	while norm(grad, 1) > grad_tol && iter < iters
		# obj = mse(θ, x_data, u_data)
		# grad = mse_grad(θ, x_data, u_data)

		# hess = mse_hess(θ, x_data, u_data)
		# reg = -1.0 * minimum(min.(eigen(mse_hess(θ, x_data, u_data)).values, 0.0)) + 1.0e-5

		Δ = grad

		α = 1.0
		while α > 1.0e-8
			θ_cand = θ - α * Δ
			obj_cand = mse(θ_cand, x_data, u_data)
			grad_cand = mse_grad(θ_cand, x_data, u_data)

			cond1 = obj_cand <= obj - 1.0e-4 * α * transpose(Δ) * grad
			# cond2 = transpose(Δ) * grad_cand <= 0.9 * transpose(Δ) * grad
			if cond1 #&& cond2
				θ = θ_cand
				obj = obj_cand
				grad = grad_cand
				break
			else
				α *= 0.5
			end
		end
		# println("alpha: $α")

		if verbose && (iter % 100 == 0)
			if α <= 1.0e-8
				@warn "line search failure"
			end
			println("iter: $iter")
			println("obj: ", mse(θ, x_data, u_data))
			println("grad. norm: ", norm(grad))
		end

		iter += 1
    end

    return θ
end

# -1.0 * minimum(min.(eigen(mse_hess(θ0, x_data, u_data)).values, 0.0))

function newton(θ0, x_data, u_data; grad_tol = 1.0e-5, iters=100, verbose=true, )
    θ = copy(θ0)
	obj = mse(θ, x_data, u_data)
	grad = mse_grad(θ, x_data, u_data)

	iter = 1
	while norm(grad, 1) > grad_tol && iter < iters
		# obj = mse(θ, x_data, u_data)
		# grad = mse_grad(θ, x_data, u_data)

		hess = mse_hess(θ, x_data, u_data)
		reg = -1.0 * minimum(min.(eigen(mse_hess(θ, x_data, u_data)).values, 0.0)) + 1.0e-5

		Δ = (hess + reg * I) \ grad

		α = 1.0
		while α > 1.0e-8
			θ_cand = θ - α * Δ
			obj_cand = mse(θ_cand, x_data, u_data)
			grad_cand = mse_grad(θ_cand, x_data, u_data)

			cond1 = obj_cand <= obj - 1.0e-4 * α * transpose(Δ) * grad
			# cond2 = transpose(Δ) * grad_cand <= 0.9 * transpose(Δ) * grad
			if cond1 #&& cond2
				θ = θ_cand
				obj = obj_cand
				grad = grad_cand
				break
			else
				α *= 0.5
			end
		end
		# println("alpha: $α")

		if verbose && (iter % 1 == 0)
			if α <= 1.0e-8
				@warn "line search failure"
			end
			println("iter: $iter")
			println("obj: ", mse(θ, x_data, u_data))
			println("grad. norm: ", norm(grad))
		end

		iter += 1
    end

    return θ
end

# optimize
θ0 = randn(np)
θ_sol = gd(θ0, x_data, u_data, α=0.05, iters=100000);
# θ_sol = gd(θ_sol, x_data, u_data, α=0.05, iters=1000);

θ_sol = gd_ls(θ0, x_data, u_data, iters=10000);
θ_sol = newton(θ0, x_data, u_data, grad_tol=1.0e-5, iters=1000);

# rollout policy
x_data[:, 1]
# x_roll = [[1.0; 0.0] + offset[2]]
x_roll = [[1.0; 0.0] + 0.00 * randn(model.n)]
u_roll = []
for t = 1:T-1
	push!(u_roll, policy(x_roll[end], θ_sol))
	push!(x_roll, fd(model, x_roll[end],
		u_roll[end],
		# u_data[:, t],
		zeros(model.d), h, t))
end

plt = plot(hcat(x_roll...)',
	width = 2.0,
	color = [:cyan :orange],
	label = ["x" "ẋ"],
	xlabel = "time step",
	ylabel = "state")

x_roll[end]

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule czmq_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("czmq")
JLLWrappers.@generate_main_file("czmq", UUID("bada3277-1da5-58a9-94c6-da212cd39369"))
end  # module czmq_jll


text:
using MAT
using jInv.Mesh
using jInv.LinearSolvers
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.Utils
using ArgParse
using OED

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table s begin

	"--matfile"
		help = "matlab file storing problem description"
	"--exName"
		help = "name of example"
	"--model"
		help = "model function. Choose: 1-> parallelTomo, 2->identity"
		arg_type = Int
		default = 1
	"--IPmode"
		help = "inverse problem model. Choose: 0->unconstrained, 1->eqconstrained, 2->nonneq, 3->box "
		arg_type = Int
		default = 0
	"--lambda"
		help = "regularization parameter for imaging problem"
		default = "1e-1"
	"--maxIter"
		help = "maximum number of iteration"
		arg_type = Int
		default = 20
	"--noiseLevel"
		help = "noiseLevel"
		default = "1e-3"
	"--nex"
		help = "Number of examples. If -1 all examples will be used. Otherwise nEx examples will be randomly choosen"
		arg_type=Int
		default =-1
	"--l"
		help = "Number of projection angles (design parameters)"
		arg_type=Int
		default=2
	"--OEDtype"
		help = "Type of OED problem: 2 -> OEDTypeB (default), 3->OEDTypeC"
		arg_type=Int
		default=2
	"--r"
		help = "Number of projections per angle, default: -1 --> pick all"
		arg_type=Int
		default=-1
	"--pad"
		help = "Padding of domain for imaging matrix"
		arg_type=Int
		default=4
	end
	return parse_args(s)
end

function main()
	parsed_args = parse_commandline()

	matfile = matread(parsed_args["matfile"])
	A       = matfile["A"]
	lambda  = parse(parsed_args["lambda"])
	if haskey(matfile,"L")
            L   = matfile["L"]
            if size(L,2)!=size(A,2); error("user supplied L has wrong size"); end
        else
	    L       = sqrt(lambda)*speye(size(A,2))
        end
	Xtrue   = matfile["x_true"]
	r       = parsed_args["r"]
	l       = parsed_args["l"]
	n       = Int(matfile["n"])
	nex     = parsed_args["nex"]
	OEDtype = parsed_args["OEDtype"]

	Ce = matfile["Ce"]
	ce = Ce*Xtrue
	Ci = matfile["Ci"]
	ci = matfile["ci"]

	if nex!==-1
		idx = randperm(size(Xtrue,2))[1:nex]
		Xtrue = Xtrue[:,idx]
		ci   = ci[:,idx]
		ce   = ce[:,idx]
	end

	Mimg  = getRegularMesh([-1 1 -1 1.],[n,n])
	# get padded domain
	pad   = parsed_args["pad"]
	domain = copy(Mimg.domain);
	domain[1:2:end] -= pad*Mimg.h
	domain[2:2:end] += pad*Mimg.h
	Mgrid = getRegularMesh(domain,[n+2*pad,n+2*pad])

	if r==-1
		r = n+2*pad
	end

	print("Mimg.h=$(Mimg.h)\tvs.\tMgrid.h=$(Mgrid.h)")
	m0 = collect(linspace(1,179,l+2)[2:end-1])
	if OEDtype==2
		# divide problems among available workers
		ids = mod(0:size(Xtrue,2)-1,nworkers())+1
		OEDparam = Array{OEDTypeBParam}(nworkers())
		for k=1:nworkers()
			OEDparam[k] = getOEDTypeBParam(r,l,Xtrue[:,ids.==k],Mimg,Mgrid)
		end
	elseif OEDtype==3
		# divide problems among available workers
		ids = mod(0:size(Xtrue,2)-1,nworkers())+1
		OEDparam = Array{OEDTypeCParam}(nworkers())
		for k=1:nworkers()
			OEDparam[k] = getOEDTypeCParam(r,l,Xtrue[:,ids.==k],Mimg,Mgrid,L=copy(L))
		end
		m0 = [m0; lambda]
		L  = 0*L
	end
	display(OEDparam[1])

	exName  = parsed_args["exName"]

	ipDict  = Dict(0=>:unconstr, 1=>:eqconstr, 2=>:nonnegconstr, 3=>:boxconstr)
	IPmode  = ipDict[parsed_args["IPmode"]]

	maxIter = parsed_args["maxIter"]
	m0 = [ 1.0 89.0; 10 100; 45. 135.; 44.9 45.1; 90. 90]'

	for k=1:size(m0,2)
		mOpt = runOEDTypeB(OEDparam,Xtrue,L,IPmode,vec(m0[:,k]),Ce=Ce,Ci=Ci,ce=ce,ci=ci,maxIter=maxIter,exName=@sprintf "%s-k-%d" exName k )
		println("m0  : $m0")
		println("mOpt: $mOpt")
	end
end

main()


text:
module DashBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Dash Block" => Ahorn.EntityPlacement(
        Maple.DashBlock,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    )
)

Ahorn.editingOptions(entity::Maple.DashBlock) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::Maple.DashBlock) = 8, 8
Ahorn.resizable(entity::Maple.DashBlock) = true, true

Ahorn.selection(entity::Maple.DashBlock) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.DashBlock, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity)

end


text:
#---------------------------------------------------------------------------------------------------
# Low-level multiplication
#---------------------------------------------------------------------------------------------------
export tlmul!, trmul!, tcmul!, tcmul
function tlmul!(
    λ::AbstractVector{<:Number},
    Γ::AbstractArray
)
    α = size(Γ, 1)
    Γ_r = reshape(Γ, α, :)
    lmul!(Diagonal(λ), Γ_r)
end
#---------------------------------------------------------------------------------------------------
function trmul!(
    Γ::AbstractArray,
    λ::AbstractVector{<:Number}
)
    β = size(Γ)[end]
    Γ_r = reshape(Γ, :, β)
    rmul!(Γ_r, Diagonal(λ))
end
#---------------------------------------------------------------------------------------------------
function tcmul!(
    Γ_n::AbstractArray{<:Number, 3},
    cmat::AbstractMatrix,
    Γ::AbstractArray{<:Number, 3}
)
    @tensor Γ_n[:] = cmat[-2,1] * Γ[-1,1,-3]
end
#---------------------------------------------------------------------------------------------------
function tcmul(
    cmat::AbstractMatrix,
    Γ::AbstractArray{<:Number, 3}
)
    @tensor Γ_n[:] := cmat[-2,1] * Γ[-1,1,-3]
    Γ_n
end

#---------------------------------------------------------------------------------------------------
# Grouping
#---------------------------------------------------------------------------------------------------
export tgroup, tgroup!
function tgroup!(
    Γ::AbstractArray{<:Number, 3},
    ΓA::AbstractArray{<:Number, 3},
    ΓB::AbstractArray{<:Number, 3}
)
    α, β = size(ΓA, 1), size(ΓB, 3)
    d1, d2 = size(ΓA, 2), size(ΓB, 2)
    ΓA_r = reshape(ΓA, α * d1, :)
    ΓB_r = reshape(ΓB, :, d2 * β)
    Γ_r = reshape(Γ, α * d1, d2 * β)
    mul!(Γ_r, ΓA_r, ΓB_r)
end
#---------------------------------------------------------------------------------------------------
function tgroup(
    ΓA::AbstractArray{<:Number, 3},
    ΓB::AbstractArray{<:Number, 3}
)
    α, β = size(ΓA, 1), size(ΓB, 3)
    d1, d2 = size(ΓA, 2), size(ΓB, 2)
    ΓA_r = reshape(ΓA, α * d1, :)
    ΓB_r = reshape(ΓB, :, d2 * β)
    reshape(ΓA_r * ΓB_r, α, d1 * d2, β)
end
#---------------------------------------------------------------------------------------------------
function tgroup!(
    Γ::AbstractArray{<:Number, 3},
    ΓA::AbstractArray{<:Number, 3},
    ΓB::AbstractArray{<:Number, 3},
    ΓC::AbstractArray{<:Number, 3}
)
    α, β = size(ΓA, 1), size(ΓC, 3)
    d1, d2, d3 = size(ΓA, 2), size(ΓB, 2), size(ΓC, 2)
    Γ_r = reshape(Γ, α, d1, d2, d3, β)
    @tensor Γ_r[:] = ΓA[-1,-2,1] * ΓB[1,-3,2] * ΓC[2,-4,-5]
end
#---------------------------------------------------------------------------------------------------
function tgroup(
    ΓA::AbstractArray{<:Number, 3},
    ΓB::AbstractArray{<:Number, 3},
    ΓC::AbstractArray{<:Number, 3}
)
    @tensor Γ_r[:] := ΓA[-1,-2,1] * ΓB[1,-3,2] * ΓC[2,-4,-5]
    reshape(Γ_r, size(ΓA, 1), :, size(ΓB, 3))
end

#---------------------------------------------------------------------------------------------------
# Tensor SVD
#---------------------------------------------------------------------------------------------------
export svd_trim, tsvd
function svd_trim(
    mat::AbstractMatrix,
    bound::Integer,
    cutoff::Real,
    renormalize::Bool
)
    res = svd(mat)
    vals = res.S
    len = min(count(i -> (i > cutoff), vals), bound)
    U = res.U[:, 1:len]
    S = vals[1:len]
    V = res.Vt[1:len, :]
    if renormalize
        S ./= norm(S)
    end
    U, S, V
end
#---------------------------------------------------------------------------------------------------
function tsvd(
    T::AbstractArray,
    α::Integer,
    d1::Integer,
    d2::Integer,
    β::Integer;
    bound::Integer=1000,
    cutoff::Real=1e-7,
    renormalize::Bool=false
)
    T_r = reshape(T, α * d1, d2 * β)
    U, S, V = svd_trim(T_r, bound, cutoff, renormalize)
    ΓA = reshape(U, α, d1, :)
    ΓB = reshape(V, :, d2, β)
    ΓA, S, ΓB
end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

mᵢandmₜ(matout, matin, surfacenormal::SVector{N,T}, r::AbstractRay{T,N}) where {T<:Real,N} = dot(surfacenormal, direction(r)) < zero(T) ? (matout, matin) : (matin, matout)

"""
    snell(surfacenormal::AbstractVector{T}, raydirection::AbstractVector{T}, nᵢ::T, nₜ::T) -> Tuple{T,T}

`nᵢ` is the index of refraction on the incidence side of the interface.
`nₜ` is the index of refraction on the transmission side.

Returns `sinθᵢ` and `sinθₜ` according to [Snell's law](https://en.wikipedia.org/wiki/Snell%27s_law).
"""
function snell(surfacenormal::S, raydirection::S, nᵢ::T, nₜ::T) where {T<:Real,S<:AbstractArray{T}}
    # snells law: nᵢsin(θᵢ) = nₜsin(θₜ)
    # sin(θ) = sin(180-θ) so don't have to worry about reversing sign of r

    # prevent div by 0
    if surfacenormal == raydirection || surfacenormal == -raydirection
        return zero(T), zero(T)
    end

    sinθᵢ = norm(cross(surfacenormal, raydirection))
    #only interested in magnitude not direction of the vector normal to surfacenormal and raydirection.
    # If wanted vector then should use -raydirection.

    @assert zero(T) <= sinθᵢ

    sinθₜ = nᵢ / nₜ * sinθᵢ
    @assert zero(T) <= sinθₜ

    return sinθᵢ, sinθₜ
end

function reflectedray(surfacenormal::S, raydirection::S) where {T<:Real,S<:AbstractArray{T}}
    r = raydirection
    nₛ = surfacenormal
    tdot = dot(r, nₛ)

    if tdot == zero(T)
        return -r
    end

    return r - 2 * nₛ * (tdot)
end

function refractedray(incidenceindex::T, transmittedindex::T, surfacenormal::S, raydirection::S) where {T<:Real,S<:AbstractArray{T,1}}

    sinθᵢ, sinθₜ = snell(surfacenormal, raydirection, incidenceindex, transmittedindex)
    # (nᵢ, nₜ) = nᵢandnₜ(incidenceindex, transmittedindex, surfacenormal, raydirection) #redundant computation here checking for nₛ⋅r < 0 in two places.
    # Optimize later if compiler doesn't optimize this away by inlining.

    if (sinθᵢ >= transmittedindex / incidenceindex) # 100% reflectance, zero transmission
        return nothing
    else
        nₛ = surfacenormal
        r = raydirection

        sgn = sign(dot(nₛ, r))
        pₜ = r - nₛ * dot(r, nₛ)
        if sinθᵢ == zero(T)
            refractedparallel = SVector{3,T}(0.0, 0.0, 0.0)
        else
            refractedparallel = pₜ * sinθₜ / sinθᵢ
        end
        refractedperpendicular = sgn * nₛ * sqrt(one(T) - sinθₜ^2)
        refracted = refractedparallel + refractedperpendicular

        return normalize(refracted)
    end
end

"""
    fresnel(nᵢ::T, nₜ::T, sinθᵢ::T, sinθₜ::T) -> Tuple{T,T}

Returns reflectance and tranmission power coefficients according to the [Fresnel equations](https://en.wikipedia.org/wiki/Fresnel_equations).
For geometric ray tracing this coefficient can be used directly to compute intensity on the detector plane. For Huygens phase optics need to take the square root to compute the amplitude.
The power of the transmitted and refracted rays may not sum to one because of the area correction applied to the transmitted component.
The intensity per area can increase or decrease depending on the indices of refraction.

`nᵢ` is the RI of the material which the incident ray travels in, `nₜ` is the RI of the material the transmitted ray travels in.
`sinθᵢ` and `sinθₜ` are the sin of the angles of incidence and transmission respectively.
"""
function fresnel(nᵢ::T, nₜ::T, sinθᵢ::T, sinθₜ::T) where {T<:Real}
    if (sinθᵢ >= nₜ / nᵢ) # 100% reflectance, zero transmission
        return (one(T), zero(T))
    end

    cosθₜ = sqrt(one(T) - sinθₜ^2)
    cosθᵢ = sqrt(one(T) - sinθᵢ^2)

    # nᵢ*sin(θᵢ) = nₜ*sin(θₜ)
    rparallel = ((nᵢ * cosθₜ - nₜ * cosθᵢ) / (nᵢ * cosθₜ + nₜ * cosθᵢ))^2
    tparallel = (2nᵢ * cosθᵢ / (nᵢ * cosθₜ + nₜ * cosθᵢ))^2
    rperp = ((nᵢ * cosθᵢ - nₜ * cosθₜ) / (nᵢ * cosθᵢ + nₜ * cosθₜ))^2
    tperp = (2nᵢ * cosθᵢ / (nᵢ * cosθᵢ + nₜ * cosθₜ))^2

    # transmitted amplitude scale factor
    Tₐ = ((nₜ * cosθₜ) / (nᵢ * cosθᵢ))

    # assuming random polarization
    ravg = 0.5 * (rparallel + rperp)
    tavg = 0.5 * (tparallel + tperp) * Tₐ  #compensate for change in intensity per area caused by the differing angles of incidence and refraction as well as differing indices.

    return (ravg, tavg)
end

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

"""
    processintersection(opticalinterface::OpticalInterface{T}, point::SVector{N,T}, normal::SVector{N,T}, incidentray::OpticalRay{T,N}, temperature::T, pressure::T, ::Bool, firstray::Bool = false) -> Tuple{SVector{N,T}, T, T}

Processes an intersection of an [`OpticalRay`](@ref) with an [`OpticalInterface`](@ref), distinct behaviors must be implemented for each subclass of `OpticalInterface`.

`point` is the 3D intersection point in global space, `normal` is the surface normal at the intersection point.

If `test` is true then the behavior of the ray should be deterministic.
`firstray` indicates that this ray is the first segment of the trace and therefore the origin is not offset.

The values returned are the normalized direction of the ray after the intersection, the _instantaneous_ power of the ray after the intersection and the optical path length of the ray up to the intersection.

`nothing` is returned if the ray should stop here, in order to obtain the correct intensity on the detector through monte carlo integration `nothing` should be returned proportionally to create the correct power distribution.
i.e. If the interface should modulate power to 76% then 24% of calls to this function should return `nothing`.
"""
function processintersection(opticalinterface::FresnelInterface{T}, point::SVector{N,T}, normal::SVector{N,T}, incidentray::OpticalRay{T,N}, temperature::T, pressure::T, test::Bool, firstray::Bool = false) where {T<:Real,N}
    λ = wavelength(incidentray)
    mᵢ, mₜ = mᵢandmₜ(outsidematerialid(opticalinterface), insidematerialid(opticalinterface), normal, incidentray)
    nᵢ = one(T)
    nₜ = one(T)
    α = zero(T)
    if !isair(mᵢ)
        mat = glassforid(mᵢ)::OpticSim.GlassCat.Glass
        nᵢ = index(mat, λ, temperature = temperature, pressure = pressure)::T
        α = absorption(mat, λ, temperature = temperature, pressure = pressure)::T
    end
    if !isair(mₜ)
        nₜ = index(glassforid(mₜ)::OpticSim.GlassCat.Glass, λ, temperature = temperature, pressure = pressure)::T
    end
    (sinθᵢ, sinθₜ) = snell(normal, direction(incidentray), nᵢ, nₜ)
    (powᵣ, powₜ) = fresnel(nᵢ, nₜ, sinθᵢ, sinθₜ)

    incident_pow = power(incidentray)

    # optical distance from ray origin to point of intersection in mm. Compensate for the fact that the ray has been slightly shortened.
    geometricpathlength = norm(point - origin(incidentray)) + (firstray ? zero(T) : T(RAY_OFFSET))
    thisraypathlength = nᵢ * geometricpathlength
    raypathlength = pathlength(incidentray) + thisraypathlength

    # compute updated power based on absorption coefficient of material using Beer's law
    internal_trans = one(T)
    if α > zero(T)
        internal_trans = exp(-α * geometricpathlength)
    end

    # TODO - this is an approximation (total hack) for now until we get better modeling of thin film reflectors, antireflection coatings, etc.
    powᵣ = max(powᵣ, reflectance(opticalinterface)) * internal_trans
    powₜ = powₜ * transmission(opticalinterface) * internal_trans

    # generate new rays using Monte Carlo sampling proportional to power. For most optical surfaces the vast majority of rays will be refracted rays.
    # So could leave this turned on all the time with little impact on performance and get approximate scattering effects for free.
    r = !test * rand()
    # assuming (powᵣ + powₜ) <= 1 (asserted in constructor)
    if r < powₜ
        # refraction
        raydirection = refractedray(nᵢ, nₜ, normal, direction(incidentray))
        raypower = powₜ * incident_pow
    elseif r < (powᵣ + powₜ)
        # reflection
        raypower = powᵣ * incident_pow
        raydirection = reflectedray(normal, direction(incidentray))
    else
        return nothing
    end

    if raydirection === nothing
        return nothing
    else
        return normalize(raydirection), raypower, raypathlength
    end
end


text:
using QuantLib

abstract OptimizationMethod
abstract CostFunction
abstract Constraint

const FINITE_DIFFERENCES_EPSILON = 1e-8

type Projection
  actualParameters::Vector{Float64}
  fixedParameters::Vector{Float64}
  fixParams::BitArray
  numberOfFreeParams::Int
end

function Projection(parameterValues::Vector{Float64}, fixParams::BitArray)
  # get num of free params
  numFree = 0
  for i in fixParams
    if !i
      numFree += 1
    end
  end

  return Projection(parameterValues, parameterValues, fixParams, numFree)
end

function project(proj::Projection, params::Vector{Float64})
  projectedParams = Vector{Float64}(proj.numberOfFreeParams)

  i = 1
  for j = 1:length(proj.fixParams)
    if !proj.fixParams[j]
      projectedParams[i] = params[j]
      i += 1
    end
  end

  return projectedParams
end

function include_params(proj::Projection, params::Vector{Float64})
  y = copy(proj.fixedParameters)
  i = 1
  for j = 1:length(y)
    if !proj.fixParams[j]
      y[j] = params[i]
      i += 1
    end
  end
  return y
end

type NoConstraint <: Constraint end
type PositiveConstraint <: Constraint end
type BoundaryConstraint <: Constraint
  low::Float64
  high::Float64
end

type ProjectedConstraint{C <: Constraint} <: Constraint
  constraint::C
  projection::Projection
end

type EndCriteria
  maxIterations::Int
  maxStationaryStateIterations::Int
  rootEpsilon::Float64
  functionEpsilon::Float64
  gradientNormEpsilon::Float64
end

test{T}(::NoConstraint, ::Vector{T}) = true

test(c::ProjectedConstraint, x::Vector{Float64}) = test(c.constraint, include_params(c.projection, x))

function test(::PositiveConstraint, x::Vector{Float64})
  for i = 1:length(x)
    if x[i] <= 0.0
      return false
    end
  end

  return true
end

function test(c::BoundaryConstraint, x::Vector{Float64})
  for i = 1:length(x)
    if x[i] < c.low || x[i] > c.high
      return false
    end
  end

  return true
end

function update{T}(constraint::Constraint, params::Vector{T}, direction::Vector{Float64}, beta::Float64)
  diff = beta
  new_params = params + diff * direction
  valid = test(constraint, new_params)
  icount = 0
  while !valid
    if (icount > 200)
      error("Can't update parameter vector")
    end

    diff *= 0.5
    icount += 1
    new_params = params + diff * direction
    valid = test(constraint, new_params)
  end

  params += diff * direction
  return params
end

## Cost Function methods ##
function get_jacobin!(cf::CostFunction, jac::Matrix{Float64}, x::Vector{Float64})
  eps_ = FINITE_DIFFERENCES_EPSILON
  xx = zeros(length(x))
  for i = 1:length(x)
    xx[i] += eps_
    fp = QuantLib.func_values(cf, xx)
    xx[i] -= 2.0 * eps_
    fm = QuantLib.func_values(cf, xx)
    for j = 1:length(fp)
      jac[j,i] = 0.5 * (fp[j] - fm[j]) / eps_
    end

    xx[i] = x[i]
  end
  return jac
end


text:
filename() = "day04.txt"

parse_data() = let data = String(read(filename()))
    strip.(split(data, "\n\n"))
end

parse_passport(passport) = Dict(f[1]=>f[2] for f in split.(split(passport, r"[ \n]"), ":"))

# Check and make sure that all required password fields are present
function part1()
    required_fields = ["byr","iyr","eyr","hgt","hcl","ecl","pid"]
    valid_passport(passport) = all(haskey(passport, k) for k in required_fields)
    passports = parse_passport.(parse_data())
    count(valid_passport.(passports))
end

# Check password rules for each field
function part2()
    match_pattern(s, pattern) = match(pattern, s) !== nothing
    chk_byr(s) = match_pattern(s, r"^[12][0-9][0-9][0-9]$") && 1920 <= parse(Int, s) <= 2002
    chk_iyr(s) = match_pattern(s, r"^[12][0-9][0-9][0-9]$") && 2010 <= parse(Int, s) <= 2020
    chk_eyr(s) = match_pattern(s, r"^[12][0-9][0-9][0-9]$") && 2020 <= parse(Int, s) <= 2030
    chk_hgt(s) = let m = match(r"^([1-9][0-9]*)(cm|in)$", s)
        if m === nothing
            false
        elseif m[2] == "cm"
            150 <= parse(Int, m[1]) <= 193
        else # in
            59 <= parse(Int, m[1]) <= 76
        end
    end
    chk_hcl(s) = match_pattern(s, r"^#[0-9a-f]{6}$")
    chk_ecl(s) = s ∈ ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"]
    chk_pid(s) = match_pattern(s, r"^[0-9]{9}$")

    checkers = [chk_byr, chk_iyr, chk_eyr, chk_hgt, chk_hcl, chk_ecl, chk_pid]
    fields   = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]
    chk_all(passport) = all(haskey(passport, f) ? c(passport[f]) : false for (c,f) in zip(checkers, fields))

    passports = parse_passport.(parse_data())
    return count(chk_all.(passports))
end


text:
function formatinfo(s::Spec, x::Ptr)
    width = 2sizeof(x) + 2
    return width + paddingsize(s, width), nothing
end

function formatfield(data::Vector{UInt8}, p::Int, s::Spec, x::Ptr, ::Nothing)
    w = 2sizeof(x) + 2
    pw = paddingwidth(s, w)
    align = default(s.align, ALIGN_RIGHT)
    p = padleft(data, p, s.fill, align, pw)
    p = @copy data p "0x"
    p = hexadecimal(data, p, reinterpret(UInt, x), w - 2, false)
    p = padright(data, p, s.fill, align, pw)
    return p
end


text:
module BanyanDataFrames

using Banyan
using BanyanArrays

using DataFrames, Missings

include("df.jl")
include("gdf.jl")

# Types
export DataFrame, GroupedDataFrame

# I/O
export read_csv, write_csv, read_parquet, write_parquet, read_arrow, write_arrow

# Dataframe properties
export nrow, ncol, size, names, propertynames

# Dataframe filtering
export dropmissing, filter, unique, nonunique

# Dataframe selection and column manipulation
export getindex, setindex!, rename

# Dataframe sorting
export sort

# Dataframe joining
export innerjoin

# Grouped dataframe properties
export length, groupcols, valuecols

# Grouped dataframe methods
export groupby, select, transform, combine, subset

# Missing
export allowmissing, disallowmissing

end # module

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
# commented-out definitions are implemented in C

#abstract Any <: Any
#abstract Type{T}

#abstract Vararg{T}
#Tuple = (Any...)

#type Symbol
#    #opaque
#end

#type TypeName
#    name::Symbol
#end

#type DataType <: Type
#    name::TypeName
#    super::Type
#    parameters::Tuple
#    names::Tuple
#    types::Tuple
#    ctor
#    instance
#    size::Int32
#    abstract::Bool
#    mutable::Bool
#    pointerfree::Bool
#end

#type UnionType <: Type
#    types::Tuple
#end

#None = Union()

#type TypeVar
#    name::Symbol
#    lb::Type
#    ub::Type
#end

#type TypeConstructor
#    parameters::Tuple
#    body
#end

#abstract AbstractArray{T,N}
#abstract DenseArray{T,N} <: AbstractArray{T,N}

#type Array{T,N} <: DenseArray{T,N}
#end

#type Module
#    name::Symbol
#end

#type LambdaStaticData
#    ast::Expr
#    sparams::Tuple
#    tfunc
#    name::Symbol
#    specializations
#    inferred
#    file::Symbol
#    line::Int
#    module::Module
#end

#type Box{T}
#    contents::T
#end

#bitstype {32|64} Ptr{T}

# types for the front end

#type Expr
#    head::Symbol
#    args::Array{Any,1}
#    typ::Any
#end

#immutable LineNumberNode
#    line::Int
#end

#immutable LabelNode
#    label::Int
#end

#immutable GotoNode
#    label::Int
#end

#immutable QuoteNode
#    value
#end

#immutable TopNode
#    name::Symbol
#end

# type Task
#     parent::Task
#     last::Task
#     storage::Any
#     consumers
#     started::Bool
#     done::Bool
#     runnable::Bool
# end

import Core.Intrinsics.ccall

export
    # key types
    Any, DataType, Vararg, ANY, NTuple, None, Top,
    Tuple, Type, TypeConstructor, TypeName, TypeVar, Union, UnionType, Void,
    AbstractArray, DenseArray,
    # special objects
    Box, Function, IntrinsicFunction, LambdaStaticData, Method, MethodTable,
    Module, Nothing, Symbol, Task, Array,
    # numeric types
    Bool, FloatingPoint, Float16, Float32, Float64, Number, Integer, Int, Int8, Int16,
    Int32, Int64, Int128, Ptr, Real, Signed, Uint, Uint8, Uint16, Uint32,
    Uint64, Uint128, Unsigned,
    # string types
    Char, ASCIIString, ByteString, DirectIndexString, String, UTF8String,
    # errors
    BoundsError, DivideError, DomainError, Exception,
    InexactError, InterruptException, MemoryError, OverflowError,
    StackOverflowError, UndefRefError, UndefVarError,
    # AST representation
    Expr, GotoNode, LabelNode, LineNumberNode, QuoteNode, SymbolNode, TopNode,
    GetfieldNode, NewvarNode,
    # object model functions
    apply, fieldtype, getfield, setfield!, yieldto, throw, tuple, is, ===, isdefined,
    # arraylen, arrayref, arrayset, arraysize, tuplelen, tupleref, convert_default,
    # kwcall,
    # type reflection
    issubtype, typeof, isa,
    # typeassert, apply_type,
    # method reflection
    applicable, invoke, method_exists,
    # constants
    JULIA_HOME, nothing, Main,
    # intrinsics module
    Intrinsics
    #ccall, cglobal, abs_float, add_float, add_int, and_int, ashr_int,
    #box, bswap_int, checked_fptosi, checked_fptoui, checked_sadd,
    #checked_smul, checked_ssub, checked_uadd, checked_umul, checked_usub,
    #nan_dom_err, copysign_float, ctlz_int, ctpop_int, cttz_int,
    #div_float, eq_float, eq_int, eqfsi64, eqfui64, flipsign_int, select_value,
    #sqrt_llvm, powi_llvm,
    #fpext, fpiseq, fpislt, fpsiround, fpuiround, fptosi, fptoui,
    #fptrunc, le_float, lefsi64, lefui64, lesif64,
    #leuif64, lshr_int, lt_float, ltfsi64, ltfui64, ltsif64, ltuif64, mul_float,
    #mul_int, ne_float, ne_int, neg_float, neg_int, not_int, or_int, rem_float,
    #sdiv_int, shl_int, sitofp, sle_int, slt_int, smod_int,
    #srem_int, sub_float, sub_int, trunc_int, udiv_int, uitofp,
    #ule_int, ult_int, unbox, urem_int, xor_int, sext_int, zext_int

immutable Nothing; end
const nothing = Nothing()

const (===) = is

abstract Number
abstract Real     <: Number
abstract FloatingPoint <: Real
abstract Integer  <: Real
abstract Signed   <: Integer
abstract Unsigned <: Integer

bitstype 16 Float16 <: FloatingPoint
bitstype 32 Float32 <: FloatingPoint
bitstype 64 Float64 <: FloatingPoint

bitstype 8  Bool <: Integer
bitstype 32 Char <: Integer

bitstype 8   Int8    <: Signed
bitstype 8   Uint8   <: Unsigned
bitstype 16  Int16   <: Signed
bitstype 16  Uint16  <: Unsigned
bitstype 32  Int32   <: Signed
bitstype 32  Uint32  <: Unsigned
bitstype 64  Int64   <: Signed
bitstype 64  Uint64  <: Unsigned
bitstype 128 Int128  <: Signed
bitstype 128 Uint128 <: Unsigned

if is(Int,Int64)
    typealias Uint Uint64
else
    typealias Uint Uint32
end

abstract Exception

type BoundsError        <: Exception end
type DivideError        <: Exception end
type DomainError        <: Exception end
type OverflowError      <: Exception end
type InexactError       <: Exception end
type MemoryError        <: Exception end
type StackOverflowError <: Exception end
type UndefRefError      <: Exception end
type UndefVarError      <: Exception
    var::Symbol
end
type InterruptException <: Exception end

abstract String
abstract DirectIndexString <: String

# simple convert for use by constructors of types in Core
convert(T, x) = convert_default(T, x, convert)

type SymbolNode
    name::Symbol
    typ
    SymbolNode(name::Symbol, t::ANY) = new(name, t)
end

type GetfieldNode
    value
    name::Symbol
    typ
end

immutable ASCIIString <: DirectIndexString
    data::Array{Uint8,1}
end

immutable UTF8String <: String
    data::Array{Uint8,1}
end

typealias ByteString Union(ASCIIString,UTF8String)

include(fname::ByteString) = ccall(:jl_load_, Any, (Any,), fname)

max_stars_repo_name: jey/julia

max_issues_repo_name: alanedelman/julia

max_forks_repo_name: alanedelman/julia


text:
isdefined(Base, :__precompile__) && __precompile__()

module AxisAlgorithms

using WoodburyMatrices
using Compat

export A_ldiv_B_md!,
    A_ldiv_B_md,
    A_mul_B_md!,
    A_mul_B_md,
    A_mul_B_perm!,
    A_mul_B_perm

@doc """
`A_ldiv_B_md(F, src, dim)` solves `F\b` for slices `b` of `src` along dimension `dim`,
storing the result along the same dimension of the output.
Currently, `F` must be an LU-factorized tridiagonal matrix or a Woodbury matrix.
""" ->
A_ldiv_B_md(F, src, dim::Integer) = A_ldiv_B_md!(similar(src), F, src, dim)

@doc """
`A_mul_B_md(M, src, dim)` computes `M*x` for slices `x` of `src` along dimension `dim`,
storing the resulting vector along the same dimension of the output.
`M` must be an `AbstractMatrix`. This uses an in-place naive algorithm.
""" ->
A_mul_B_md(M::AbstractMatrix, src, dim::Integer) = A_mul_B_md!(alloc_matmul(M,src,dim), M, src, dim)

@doc """
`A_mul_B_perm(M, src, dim)` computes `M*x` for slices `x` of `src` along dimension `dim`, storing the
resulting vector along the same dimension of the output.
`M` must be an `AbstractMatrix`. This uses `permutedims` to make dimension
`dim` into the first dimension, performs a standard matrix multiplication, and restores the original
dimension ordering. In many cases, this algorithm exhibits the best cache behavior.
""" ->
A_mul_B_perm(M::AbstractMatrix, src, dim::Integer) = A_mul_B_perm!(alloc_matmul(M,src,dim), M, src, dim)

function alloc_matmul{S,N}(M,src::AbstractArray{S,N},dim)
    sz = [size(src)...]
    sz[dim] = size(M,1)
    T = Base.promote_op(@functorize(*), eltype(M), S) # functorize macro used to support julia 0.4
    Array(T, sz...)::Array{T,N}
end

include("tridiag.jl")
include("matmul.jl")
include("woodbury.jl")

end # module


text:
using Lazy
import Base:convert, colon

takeuntil(pred::Function, l::List) =
    @lazy isempty(l) ? [] :
        pred(first(l)) ? [first(l)] : first(l):takeuntil(pred, tail(l))

colon(x::List, xs::List) = prepend(x, xs)