text:
function complexmodes(z)
    z *= 10
    f = 5e9
    c = 3e8
    μ₀ = 4e-7π
    ϵ₀ = 1e-9/36/π
    a = 6.35e-3
    b = 10e-3
    ϵᵣ₁ = 10
    ϵᵣ₂ = 1
    m = 1

    ω = 2π*f
    k₀ = ω/c
    α = real(z)*k₀
    β = imag(z)*k₀
    γ = α + im*β
    ϵ₁ = ϵ₀*ϵᵣ₁
    ϵ₂ = ϵ₀*ϵᵣ₂
    μ₁ = μ₀
    μ₂ = μ₀
    κ₁ = sqrt(γ^2 + k₀^2*ϵᵣ₁)
    κ₂ = sqrt(γ^2 + k₀^2*ϵᵣ₂)
    η₁ = sqrt(μ₁/ϵ₁)
    η₂ = sqrt(μ₂/ϵ₂)

    Jm_a1 = besselj1(κ₁*a)
    Jm_a2 = besselj1(κ₂*a)
    Ym_a2 = bessely1(κ₂*a)
    Jm_b2 = besselj1(κ₂*b)
    Ym_b2 = bessely1(κ₂*b)
    DJm_a1 = (besselj0(κ₁*a) - besselj(2, κ₁*a))/2
    DJm_a2 = (besselj0(κ₂*a) - besselj(2, κ₂*a))/2
    DJm_b2 = (besselj0(κ₂*b) - besselj(2, κ₂*b))/2
    DYm_a2 = (bessely0(κ₂*a) - bessely(2, κ₂*a))/2
    DYm_b2 = (bessely0(κ₂*b) - bessely(2, κ₂*b))/2

    W = [Jm_a1                 0                       -Jm_a2              -Ym_a2              0                       0;
       0                     Jm_a1/η₁                0                   0                   -Jm_a2/η₂               -Ym_a2/η₂;
       γ*m*Jm_a1/(a*κ₁^2)    -ω*μ₁*DJm_a1/(κ₁*η₁)    -γ*m*Jm_a2/(a*κ₂^2) -γ*m*Ym_a2/(a*κ₂^2) ω*μ₂*DJm_a2/(κ₂*η₂)     ω*μ₂*DYm_a2/(κ₂*η₂);
       -ω*ϵ₁*DJm_a1/κ₁       -m*γ*Jm_a1/(a*κ₁^2*η₁)  ω*ϵ₂*DJm_a2/κ₂      ω*ϵ₂*DYm_a2/κ₂      m*γ*Jm_a2/(a*κ₂^2*η₂)   m*γ*Ym_a2/(a*κ₂^2*η₂);
       0                     0                       Jm_b2               Ym_b2               0                       0;
       0                     0                       γ*m*Jm_b2/(b*κ₂^2)  γ*m*Ym_b2/(b*κ₂^2)  -ω*μ₂*DJm_b2/(κ₂*η₂)    -ω*μ₂*DYm_b2/(κ₂*η₂)]
    w = det(W)
end

R = 1.0
r = 0.15

origcoords = diskdomain(R, r)

# matlab results from https://github.com/PioKow/GRPF for comparison
matlab_zroots = [-0.856115203791905 + 0.000000000114004im,
                 -0.775021521974263 + 0.000000000017321im,
                 -0.703772250164549 - 0.000000000209979im,
                 -0.444429043571261 - 0.000000000174298im,
                 -0.096642302493522 - 0.062923397621120im,
                 -0.096642302349294 + 0.062923397246964im,
                  0.096642302349294 - 0.062923397246964im,
                  0.096642302363736 + 0.062923397315544im,
                  0.444429043571261 + 0.000000000174298im,
                  0.703772250003777 - 0.000000000078642im,
                  0.775021521891359 - 0.000000000075087im,
                  0.856115203873180 - 0.000000000219841im]

matlab_zpoles = [0.000000000022863 - 0.100000000307523im,
                 0.000000000069120 + 0.100000000307523im]

zroots, zpoles = grpf(complexmodes, origcoords)

@test length(zroots) == 12
@test length(zpoles) == 2

@test approxmatch(zroots, matlab_zroots)
@test approxmatch(zpoles, matlab_zpoles)

pzroots, pzpoles, quadrants, phasediffs, tess, g2f = grpf(complexmodes, origcoords, PlotData())

@test approxmatch(pzroots, matlab_zroots)
@test approxmatch(pzpoles, matlab_zpoles)

# Can't test with big values because of special functions (bessel)


text:
"""
    particle dynamics
    - 3D particle subject to contact forces

    - configuration: q = (x, y, z) ∈ R³
    - impact force (magnitude): γ ∈ R₊
    - friction force: β ∈ R⁴₊
        - friction coefficient: μ ∈ R₊

    Discrete Mechanics and Variational Integrators
        pg. 363
"""
mutable struct Particle{T} <: Model{T}
    nq::Int
    nu::Int
    nw::Int
    nc::Int
    m::T # mass
    g::T # gravity
    μ_world::T # friction coefficient
	μ_joint::T

	base::BaseMethods
	dyn::DynamicsMethods

	joint_friction::SVector
end

function lagrangian(model::Particle, q, q̇)
	L = 0.0

	L += 0.5 * model.m * (transpose(q̇) * q̇)
	L -= model.m * model.g * q[3]

	return L
end

function kinematics(::Particle, q)
	return q
end

# mass matrix
function M_func(model::Particle, q)
    m = model.m

    Diagonal(@SVector [m, m, m])
end

# gravity
function C_func(model::Particle, q, q̇)
    m = model.m
    g = model.g

    @SVector [0.0, 0.0, m * g]
end

# signed distance function
function ϕ_func(model::Particle, env::Environment, q)
	SVector{1}(q[3] - env.surf(q[1:2]))
end

# control Jacobian
function B_func(model::Particle, q)
    SMatrix{3, 3}([1.0 0.0 0.0;
                   0.0 1.0 0.0;
                   0.0 0.0 1.0])
end

# disturbance Jacobian
function A_func(model::Particle, q)
    SMatrix{3, 3}([1.0 0.0 0.0;
                   0.0 1.0 0.0;
                   0.0 0.0 1.0])
end

# contact Jacobian
function J_func(model::Particle, env::Environment, q)
    SMatrix{3, 3}([1.0 0.0 0.0;
                   0.0 1.0 0.0;
				   0.0 0.0 1.0])
end

function contact_forces(model::Particle, env::Environment{<:World, LinearizedCone}, γ1, b1, q2, k)
	m = friction_mapping(env)

	SVector{3}(transpose(rotation(env, k[1:2])) * [m * b1; γ1])
end

function contact_forces(model::Particle, env::Environment{<:World, NonlinearCone}, γ1, b1, q2, k)
	m = friction_mapping(env)

	SVector{3}(transpose(rotation(env, k[1:2])) * [b1; γ1])
end

function velocity_stack(model::Particle, env::Environment{<:World, LinearizedCone}, q1, q2, k, h)
	v = J_func(model, env, q2) * (q2 - q1) / h[1]

	v_surf = rotation(env, k[1:2]) * v

	SVector{4}(transpose(friction_mapping(env)) * v_surf[1:2])
end

function velocity_stack(model::Particle, env::Environment{<:World, NonlinearCone}, q1, q2, k, h)
	v = J_func(model, env, q2) * (q2 - q1) / h[1]

	v_surf = rotation(env, k[1:2]) * v

	SVector{2}(v_surf[1:2])
end

# Dimensions
nq = 3              # configuration dimension
nu = 3              # control dimension
nw = 3              # disturbance dimension
nc = 1              # number of contact points

particle = Particle(3, 3, 3, 1, 1.0, 9.81, 1.0, 0.0,
	BaseMethods(), DynamicsMethods(),
	SVector{3}(zeros(3)))

function friction_coefficients(model::Particle)
    return [model.μ_world]
end

function initialize_z!(z, model::Particle, idx::RoboDojo.IndicesZ, q)
    z .= 1.0
    z[idx.q] .= q
end

ext: jl


text:
using EffectiveWaves, Test, LinearAlgebra

@testset "low frequency 3D acoustics" begin

spatial_dim = 3
medium = Acoustic(spatial_dim; ρ=1.2, c=1.5)

# NOTE: we so far only understand the region of the low frequency effective medium when all the particles have the same radius. That is, we get exact matches below (Monte Carlo also confirms this). However, when using species with different radius, it is not clear what effective bounding shape to use.

s1 = Specie(
    Acoustic(spatial_dim; ρ=10.2, c=10.1), Sphere(0.001);
    volume_fraction=0.2
);
s2 = Specie(
    Acoustic(spatial_dim; ρ=0.2, c=4.1), Sphere(0.001);
    volume_fraction=0.15
);
species = [s1,s2]

ωs = LinRange(1e-2,1e-1,3)
tol = 1e-7
basis_order = 1

θ = 0.0
psource = PlaneSource(medium, [sin(θ),0.0,cos(θ)]);
source = plane_source(medium; direction = [sin(θ),0.0,cos(θ)])
material = Material(Sphere(4.0),species);

basis_field_order = 6

ks = ωs ./ medium.c

opts = Dict(
    :tol => tol, :num_wavenumbers => 2,
    :mesh_size => 2.0, :mesh_points => 10,
    :basis_order => basis_order, :basis_field_order => basis_field_order
);
AP_kps = [
    wavenumbers(ω, medium, species;
        symmetry = PlanarAzimuthalSymmetry(),
        numberofparticles = material.numberofparticles,
        opts...
    )
for ω in ωs]

k_effs = [kps[1] for kps in AP_kps]

eff_medium = effective_medium(medium, species; numberofparticles = material.numberofparticles)

k_lows = ωs ./ eff_medium.c

errs = abs.(k_effs - k_lows)
@test errs[1] < errs[2] < errs[3]
@test maximum(errs) < tol
@test errs[1] < 10.0 * tol^2

A_waves = [
    WaveMode(ωs[i], k_effs[i], psource, material;
        basis_order = basis_order, basis_field_order = basis_field_order)
for i in eachindex(k_effs)]

scat_azis = material_scattering_coefficients.(A_waves);

r = maximum(outer_radius.(species))
material_low = Material(Sphere(outer_radius(material.shape) - r),species);
effective_sphere = Particle(eff_medium, material_low.shape)

Linc = basis_field_order + basis_order;

# I can not work out theoretically why this scale_number_density should be needed below.
scale_number_density = 1.0 - 1.0 / material.numberofparticles

errs =  map(eachindex(ωs)) do i
    source_coefficients =  regular_spherical_coefficients(source)(Linc,zeros(3),ωs[i])
    Tmat = MultipleScattering.t_matrix(effective_sphere, medium, ωs[i], Linc)
    Tmat = Tmat ./ scale_number_density
    norm(scat_azis[i] - Tmat * source_coefficients) / norm(scat_azis[i])
end

@test sum(errs .< maximum(outer_radius.(species)) .* abs.(ks) .* 1e-3) == length(ks)
@test errs[1] < errs[2] < errs[3]
@test errs[1] < tol * 1e-3

# Test an effective low frequency small sphere

    s1 = Specie(
        Acoustic(spatial_dim; ρ=10.2, c=10.1), Sphere(0.001);
        volume_fraction=0.15,
        exclusion_distance = 1.567
    )
    species = [s1];
    basis_field_order = 3

    R = 0.01
    material = Material(Sphere(R),species);

    opts = Dict(
        :tol => tol, :num_wavenumbers => 2,
        :mesh_size => 2.0, :mesh_points => 10,
        :basis_order => basis_order, :basis_field_order => basis_field_order
        , :numberofparticles => material.numberofparticles
    );

    kps = wavenumbers(ωs[1], medium, species; symmetry = PlanarAzimuthalSymmetry(), opts...)

    eff_medium = effective_medium(medium, species; numberofparticles = material.numberofparticles)
    k_low = ωs[1] / eff_medium.c

    @test abs(kps[1] - k_low) / abs(k_low) < 1e-10

    wavemode = WaveMode(ωs[1], kps[1], psource, material;
        basis_order = basis_order,
        basis_field_order = basis_field_order
    )

    # I can not work out theoretically why this scale_number_density should be needed below. However, numerically it clearer is needed.
    scale_number_density = 1.0 - 1.0 / material.numberofparticles
    scat_azi = material_scattering_coefficients(wavemode);

    r = maximum(outer_radius.(species))
    material_lows = [Material(Sphere(R + r1),species) for r1 in (-3.0*r):r/20.0:(r)];

    effective_spheres = map(material_lows) do m
        Particle(eff_medium, m.shape)
    end;

    Linc = basis_field_order + basis_order;

    source_coefficients =  regular_spherical_coefficients(source)(Linc,zeros(3),ωs[1])

    errs =  map(eachindex(material_lows)) do i
        Tmat = MultipleScattering.t_matrix(effective_spheres[i], medium, ωs[1], Linc)
        Tmat = Tmat ./ scale_number_density
        norm(scat_azi - Tmat * source_coefficients) / norm(source_coefficients)
    end

    err, i = findmin(errs)
    @test err < 1e-20
    @test outer_radius(material_lows[i].shape) == R - r

# Test the radius for the scattering coefficients seperate from eigenvectors

    wavemodes = [
        EffectiveRegularWaveMode(ωs[1], kps[1], psource, m, wavemode.eigenvectors;
            basis_order = wavemode.basis_order, basis_field_order = wavemode.basis_field_order)
    for m in material_lows]

    # The scattering coefficients depend on the radius of the material
    scat_azis = material_scattering_coefficients.(wavemodes);

    material_low = Material(Sphere(R - r),species);
    effective_sphere = Particle(eff_medium, material_low.shape);

    Tmat = MultipleScattering.t_matrix(effective_sphere, medium, ωs[1], Linc) ./ scale_number_density

    errs =  map(eachindex(material_lows)) do j
        norm(scat_azis[j] - Tmat * source_coefficients) / norm(source_coefficients)
    end;

    err, j = findmin(errs)
    @test err < 1e-20
    # We expect the min error to be for a materail (which contains all particles) with radius R
    @test outer_radius(material_lows[j].shape) == R

# Test radius for boundary condition seperate from the scattering coefficients

    wavemodes = [
        WaveMode(ωs[1], kps[1], psource, m;
            basis_order = basis_order, basis_field_order = basis_field_order)
    for m in material_lows];

    wavemodes_2 = [
        EffectiveRegularWaveMode(ωs[1], kps[1], psource, material, wavemodes[j].eigenvectors;
            basis_order = wavemodes[j].basis_order, basis_field_order = wavemodes[j].basis_field_order)
    for j in eachindex(material_lows)];

    scat_azis = material_scattering_coefficients.(wavemodes_2);

    errs =  map(eachindex(material_lows)) do j
        norm(scat_azis[j] - Tmat * source_coefficients) / norm(source_coefficients)
    end

    # The boundary conditions for this case are not at all sensitive to small changes in the radius. For this reason, here we do things differently.

    # Let us find the case with the radius R we expect to be accurate
    j = findfirst(R .== [outer_radius(m.shape) for m in material_lows])
    err = errs[j]

    @test err < 1e-25
    @test err < errs[1]
    @test err < errs[end]

end

max_stars_repo_name: arturgower/EffectiveWaves


text:
# Autogenerated wrapper script for XNNPACK_jll for x86_64-w64-mingw32
export libxnnpack

using CPUInfo_jll
using PThreadPool_jll
JLLWrappers.@generate_wrapper_header("XNNPACK")
JLLWrappers.@declare_library_product(libxnnpack, "libXNNPACK.dll")
function __init__()
    JLLWrappers.@generate_init_header(CPUInfo_jll, PThreadPool_jll)
    JLLWrappers.@init_library_product(
        libxnnpack,
        "bin\\libXNNPACK.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module ThinModules
using DataStructures

export @thinmod

macro thinmod(module_def)
    if __module__ == Base.__toplevel__
        modname = [module_def.args[2]]
    else
        basedname = find_based(__module__)

        modname = [basedname..., module_def.args[2]]
    end
    context = Context(modname, __module__)

    enqueue!(context.queue, (modname, module_def.args[3], String(__source__.file::Symbol)))
    while !isempty(context.queue)
        (mod, defs, file) = dequeue!(context.queue)
        context.modules_body[mod] = analysis_module!(context, defs, mod, file)
    end

    sorts = topological_sort(context)

    modules = create_modules_hierarchy!(context)
    for m in sorts
        Base.eval(modules[m], context.modules_body[m])
    end
    for (m, fs) in context.file_dependencies
        for f in fs
            @eval modules[m] Base.include_dependency($f)
        end
    end
end
const ModuleName = Array{Symbol}

struct Context
    root_name::ModuleName
    based::Module
    queue::Queue{Tuple{ModuleName,Expr,String}}
    modules_body::Dict{ModuleName,Expr}
    dependencies::Dict{ModuleName,Set{ModuleName}}
    file_dependencies::Dict{ModuleName,Set{String}}
    function Context(root_name::ModuleName, based::Module)
        queue = Queue{Tuple{ModuleName,Expr,String}}()
        modules_body = Dict{ModuleName,Expr}()
        dependencies = Dict{ModuleName,Array{ModuleName}}()
        file_dependencies = Dict{ModuleName,Array{String}}()
        return new(root_name, based, queue, modules_body, dependencies, file_dependencies)
    end
end

function minclude(path::String)
    return Meta.parseall(read(path, String); filename=path)
end

function find_based(m::Module)
    name = [Base.nameof(m)]
    while Base.nameof(Base.parentmodule(m)) != Base.nameof(m)
        push!(name, Base.nameof(Base.parentmodule(m)))
        m = Base.parentmodule(m)
    end
    reverse!(name)
    return name
end

function topological_sort(context::Context)
    indegree = Dict{ModuleName,Int}()
    rev_deps = rev_dependencies(context)
    for m in keys(context.modules_body)
        indegree[m] = 0
    end

    for m in values(rev_deps)
        for d in m
            indegree[d] += 1
        end
    end
    queue = Queue{ModuleName}()
    sorts = ModuleName[]
    for m in keys(context.modules_body)
        if indegree[m] == 0
            enqueue!(queue, m)
        end
    end

    counts = 0
    while !isempty(queue)
        m = dequeue!(queue)
        counts += 1
        push!(sorts, m)
        if haskey(rev_deps, m)
            for d in rev_deps[m]
                indegree[d] -= 1
                if indegree[d] == 0
                    enqueue!(queue, d)
                end
            end
        end
    end
    if counts != length(keys(context.modules_body))
        error("Cyclic dependencies")
    end
    return sorts
end

function create_module!(base::Module, name::Symbol)
    return @eval base module $name end
end

function create_modules_hierarchy!(context::Context)
    mods = Dict{ModuleName,Module}()
    to_create = sort!(collect(keys(context.modules_body)), by=length)
    mods[to_create[1]] = create_module!(context.based, to_create[1][end])
    for mod in to_create[2:end]
        based = mods[mod[1:end - 1]]
        mods[mod] = create_module!(based, mod[end])
    end
    return mods
end
# >---------------------------- import ---------------------------
function resolve_import(expr::Expr, base::ModuleName)
    if expr.head == :(:) || expr.head == :as
        path = expr.args[1].args
    elseif expr.head == :(.)
        path = expr.args
    else
        error("Invalid import expression")
    end
    if path[1] != :(.)
        return convert.(Symbol, path)
    else
        count = 0
        i = 2
        while path[i] == :(.)
            count += 1
            i += 1
        end
        return [base[1:end - count]..., path[count + 2:end]...]
    end
end

function isimport(expr::Expr)
    return expr.head == :using || expr.head == :import
end

function isprefix(prefix::ModuleName, name::ModuleName)
    return length(prefix) <= length(name) && prefix == name[1:length(prefix)]
end

function find_module(context::Context, mod::ModuleName)
    # This is for `import A.x as y`
    # If we find something is not an existing module, we
    # will find its nearest parent module and add it.

    seen = keys(context.modules_body)
    for i in length(mod):-1:1
        if mod[1:i] ∈ seen
            return mod[1:i]
        end
    end
    error("No such module exists.")
end

function rev_dependencies(context::Context)
    rev_deps = Dict{ModuleName,Set{ModuleName}}()
    for m in keys(context.modules_body)
        rev_deps[m] = Set{ModuleName}()
    end
    for (m, deps) in context.dependencies
        for d in deps
            push!(rev_deps[find_module(context, d)], m)
        end
    end
    return rev_deps
end

# <---------------------------- import ---------------------------

function analysis_block!(context::Context, block::Expr, module_name::ModuleName, new_args::Array{Any}, filename::String)
    i = 1
    while i <= length(block.args)
        current = block.args[i]
        if current isa LineNumberNode
            push!(new_args, current)
        elseif current isa Expr
            if current.head == :call && current.args[1] == :include
                to_include_path = abspath(joinpath(dirname(filename), current.args[2]))
                push!(context.file_dependencies[module_name], to_include_path)
                analysis_block!(context, minclude(to_include_path), module_name, new_args, to_include_path)
            elseif isimport(current)
                for imp in current.args
                    nmodule_name = resolve_import(imp, module_name)
                    if isprefix(context.root_name, nmodule_name)
                        push!(context.dependencies[module_name], nmodule_name)
                    end
                end

                push!(new_args, current)
            elseif current.head == :module
                enqueue!(context.queue, ([module_name..., current.args[2]], current.args[3], filename))
            else
                push!(new_args, current)
            end
        end
        i += 1
    end
end

function analysis_module!(context::Context, def::Expr, name::ModuleName, filename::String)
    i = 1
    new_args = Union{Any}[]
    context.file_dependencies[name] = Set{String}()
    context.dependencies[name] = Set{ModuleName}()

    analysis_block!(context, def, name, new_args, filename)

    new_block = Expr(:block, new_args...)
    return new_block
end
end


text:
using Multilane
using MCTS
using POMDPs
using POMDPModelTools
using POMDPSimulators
using Random
using Test
using Cairo

#Set up problem configuration
nb_lanes = 2
pp = PhysicalParam(nb_lanes,lane_length=200.) #2.=>col_length=8
r_crash = -1.
accel_cost = -1e-5
decel_cost = -5e-6
invalid_cost = -1e-5
lineride_cost = -1e-5
lanechange_cost = -2e-5
_discount = 1.
nb_cars=3

rmodel = OriginalRewardModel(r_crash,
                             accel_cost,
                             decel_cost,
                             invalid_cost,
                             lineride_cost,
                             lanechange_cost)

dmodel = IDMMOBILModel(nb_cars, pp)

mdp = OriginalMDP(dmodel, rmodel, _discount, true);

# initialize the solver
# the hyper parameters in MCTS can be tricky to set properly
# n_iterations: the number of iterations that each search runs for
# depth: the depth of the tree (how far away from the current state the algorithm explores)
# exploration constant: this is how much weight to put into exploratory actions.
# A good rule of thumb is to set the exploration constant to what you expect the upper bound on your average
#  expected reward to be.

solver = MCTSSolver(n_iterations=100, depth=10, exploration_constant=1.0, enable_tree_vis=true)

# initialize the policy by passing in your problem and the solver
policy = solve(solver, mdp);

sim = HistoryRecorder(rng=MersenneTwister(9), max_steps=10) # initialize a random number generator

simulate(sim, mdp, policy, rand(sim.rng, initialstate(mdp)))


text:
## (c) 2016 David A. van Leeuwen

## Unit tests for ../src/arithmetic.jl

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

## test arithmetic operations
print("arithmetic, ")

dotops = [:.+, :.-, :.*, :./]

x = NamedArray(randn(5, 10))
@test @inferred -x == -(x.array)
for op in dotops
	@eval begin
		for y in (true, Int8(3), Int16(3), Int32(3), Int64(3), Float32(π), Float64(π), BigFloat(π))
			z = ($op)(x, y)
			@test z.array == ($op)(x.array, y)
			@test names(z) == names(x)
			z = ($op)(y, x)
			@test z.array == ($op)(y, x.array)
			@test names(z) == names(x)
		end
	end
end
for op in dotops
	for T in (Bool, Int8, Int16, Int32, Int64, Float32, Float64)
		@eval begin
			z = ($op)(x, x)
			@test z.array == ($op)(x.array, x.array)
			@test names(z) == names(x)
			@test dimnames(z) == dimnames(x)
			y = rand($T, 5, 10)
			z = ($op)(x, y)
			@test z.array == ($op)(x.array, y)
			@test names(z) == names(x)
			@test dimnames(z) == dimnames(x)
			z = ($op)(y, x)
			@test z.array == ($op)(y, x.array)
			@test names(z) == names(x)
			@test dimnames(z) == dimnames(x)
		end
	end
end

include("init-namedarrays.jl")

@test n / π == π \ n

@test (v - (1:6)).array == v.array - (1:6)
@test ((1:6) - v).array == (1:6) - v.array

## matmul
## Assume dimensions/names are correct
c = NamedArray(Float64, 5, 5)
r5x10 = randn(5, 10)
r10x5 = randn(10, 5)
@test mul!(c, r5x10, r10x5) ≈ r5x10 * r10x5
@test mul!(c, r5x10, adjoint(r5x10)) ≈ r5x10 * r5x10'
@test mul!(c, r5x10, transpose(r5x10)) ≈ r5x10 * r5x10'
@test mul!(c, adjoint(r10x5), r10x5) ≈ r10x5' * r10x5
@test mul!(c, transpose(r10x5), r10x5) ≈ r10x5' * r10x5
@test mul!(c, transpose(r10x5), transpose(r5x10)) ≈ r10x5' * r5x10'

for M in (NamedArray(rand(4)), NamedArray(rand(4,3)))
    @test M * M' ≈ M.array * M.array'
    value = M' * M
    @test (length(value) == 1 ? value[1] : value) ≈ M.array' * M.array # M' * M is always a NamedArray
    @test n * M ≈ n * M.array ≈ n.array * M ≈ n.array * M.array
    ## the first expression dispatches Ac_Mul_Bc!:
    @test isapprox(M' * n.array', M' * n')
    # TODO : @test isapprox(M' * n', M.array' * n')
    # TODO : @test isapprox(M.array' * n', M.array' * n.array')
end

## bug #34
@test unique(names(n * n'))[1] == names(n, 1)
@test unique(names(n' * n))[1] == names(n, 2)

## \
v = NamedArray(randn(2))
m = NamedArray(randn(2, 5))

@test (v \ v)[1] == v.array \ v == v.array \ v.array # v\v returns a NamedMatrix with one element

@test (n \ v).array == n.array \ v.array
@test names(n \ v, 1) == names(n, 2)
@test dimnames(n \ v, 1) == dimnames(n, 2)

@test (v \ n).array == v.array \ n.array
@test names(v \ n, 2) == names(n, 2)
@test dimnames(v \ n, 2) == dimnames(n, 2)

@test (v \ n.array) == v.array \ n.array
@test (n \ m.array).array == n.array \ m.array
@test names(n \ m.array, 1) == names(n, 2)
@test dimnames(n \ m.array, 1) == dimnames(n, 2)

@test (n \ m).array == n.array \ m.array
@test names(n \ m, 1) == names(n, 2)
@test names(n \ m, 2) == names(m, 2)

@test n.array \ v == n.array \ v.array
@test (v.array \ n).array == v.array \ n.array
@test names(v.array \ n, 2) == names(n, 2)
@test dimnames(v.array \ n, 2) == dimnames(n, 2)

@test (n.array \ m).array == n.array \ m.array
@test names(n.array \ m, 2) == names(m, 2)
@test names(n.array \ m, 1) == NamedArrays.defaultnames(n.array, 2)
@test dimnames(n.array \ m, 2) == dimnames(m, 2)

m = NamedArray((x = randn(100, 10); x'x))
for f in (:inv, :sqrt, :pinv, :exp) # TODO : :chol
	@eval fm = ($f)(m)
	@test @eval fm.array == ($f)(m.array)
	@test names(fm) == names(m)
	@test dimnames(fm) == dimnames(m)
end

## tril, triu
m = NamedArray(rand(5,5))
c = copy(m)
tril!(c)
@test tril(m).array == tril(m.array) == c.array
c = copy(m)
triu!(c)
@test triu(m).array == triu(m.array) == c.array

## lufact!
lufn = lu(n)
lufa = lu(n.array)
@test lufn.U == lufa.U
@test lufn.L == lufa.L
@test names(lufn.U, 2) == names(n, 2) # lu doesn't return NamedArrays now
@test dimnames(lufn.U, 2) == dimnames(n, 2)
@test names(lufn.L, 1) == names(n, 1)
@test dimnames(lufn.L, 1) == dimnames(n, 1)
@test lufn.p == lufa.p
@test lufn.P == lufa.P

a = randn(1000,10); s = NamedArray(a'a)

## The necessity for isapprox sugests we don't get BLAS implementations but a fallback...
# TODO: @test isapprox(cholfact(s).factors.array, cholfact(s.array).factors)

@test isapprox(qr(s).factors.array, qr(s.array).factors)
# @test isapprox(qr(s).T, qr(s.array).T)

@test isapprox(eigen(s).values, eigen(s.array).values)
@test eigen!(copy(s)).values == eigen(s).values

@test isapprox(eigen(s).vectors, eigen(s.array).vectors)
@test eigen!(copy(s)).vectors == eigen(s).vectors

@test eigvals(s) == eigvals!(copy(s)) == eigvals(s.array)

@test isapprox(hessenberg(s).factors, hessenberg(s.array).factors)
@test isapprox(hessenberg(s).τ, hessenberg(s.array).τ)
@test isapprox(hessenberg!(copy(s)).factors, hessenberg(s.array).factors)
@test isapprox(hessenberg!(copy(s)).τ, hessenberg(s.array).τ)

s2 = randn(10,10)
for field in [:T, :Z, :values]
	@test isapprox(getproperty(schur(s), field), getproperty(schur(s.array), field))
	@test isapprox(getproperty(schur!(copy(s)), field), getproperty(schur(s.array), field))
	@test isapprox(getproperty(schur(s, s2), field), getproperty(schur(s.array, s2), field))
	@test isapprox(getproperty(schur!(copy(s), copy(s2)), field), getproperty(schur(s.array, s2), field))
end

@test isapprox(svd(s).U, svd(s.array).U)
@test isapprox(svd(s).S, svd(s.array).S)
@test isapprox(svd(s).Vt, svd(s.array).Vt)

@test svdvals(s) == svdvals!(copy(s)) == svdvals(s.array)

# TODO : @test diag(s).array == diag(s.array)
# TODO : @test names(diag(s), 1) == names(s, 1)
# TODO : @test dimnames(diag(s), 1) == dimnames(s, 1)

# TODO : @test diagm(v).array == diagm(v.array)
# TODO : @test names(diagm(v)) == names(v)[[1,1]]

@test cond(n) == cond(n.array)

@test kron(n, n').array == kron(n.array, n.array')

a = randn(10, 10)
b = randn(10, 10)
@test lyap(s, a) == lyap(s.array, a)
@test sylvester(s, a, b).array == sylvester(s.array, a, b)

@test isposdef(s) == isposdef(s.array)


text:
###################################################
## geo/data/filters.jl
## Filtering geo-trips
###################################################

"""
    `removeBadTrips`: filter a `GeoTrip` list to remove outliers
"""
function removeOutliers(trips::GeoData)
    reg = Bool[isRegular(t) for t in trips]
    t = trips[reg]
    @printf("%.2f%% outliers removed\n", 100*(1-length(t)/length(trips)))
    return t
end

"""
    `isRegular`: returns false if the trip is an "outlier"
"""
function isRegular(t::GeoTrip)
    #trips should be > 30s and <3h
    if t.time <= 30. || t.time >= 3*3600.
        return false
    end
    d = tripDistance(t)
    #trips should be > 250m and < 200km
    if d <= 250. || d >= 200_000
        return false
    end
    #  direct speed <= 110 km/h, >= 2km/h
    s = d*3.6/t.time
    if s > 110. || s < 2.
        return false
    end
    return true
end

"""
    `tripDistance`: compute line distance for one GeoTrip
"""
tripDistance(t::GeoTrip) = distanceGeo(t.pLon,t.pLat,t.dLon,t.dLat)

"""
    `inTimeWindow`: keep trips in a certain time window
    - time window is of the form "13:40","13:52" (here minute "52" is NOT included)
    - for GeoTrip object: returns boolean
    - for GeoData object: returns filtered GeoData object
"""
function inTimeWindow(t::GeoTrip, startHour::Int, startMin::Int, endHour::Int, endMin::Int)
    return ((startHour < Dates.hour(t.pTime) ||
        (startHour == Dates.hour(t.pTime) && startMin <= Dates.minute(t.pTime))) &&
    (endHour > Dates.hour(t.pTime) ||
        (endHour   == Dates.hour(t.pTime) && endMin > Dates.minute(t.pTime))))
end

function inTimeWindow(trips::GeoData, startHour::Int, startMin::Int, endHour::Int, endMin::Int)
    mask = falses(length(trips))
    for (i, t) in enumerate(trips)
        if i%10_000 == 0
            @printf("\r%.2f%% trips checked     ",100*i/length(trips))
        end
        mask[i] = inTimeWindow(t, startHour, startMin, endHour, endMin)
    end
    newTrips = trips[mask]
    @printf("\r%.2f%% trips removed\n", 100*(1-length(newTrips)/length(trips)))
    return newTrips
end

function inTimeWindow(trips::GeoData, startTime::AbstractString, endTime::AbstractString)
    s = split(startTime,":")
    e = split(endTime,  ":")
    startHour = parse(Int, s[1])
    endHour   = parse(Int, e[1])
    startMin = parse(Int, s[2])
    endMin   = parse(Int, e[2])
    return inTimeWindow(trips, startHour, startMin, endHour, endMin)
end

"""
    `onlyWeekdays`: keep trips that occur on weekdays
    - for GeoTrip object: returns boolean
    - for GeoData object: returns filtered GeoData object
"""
function onlyWeekdays(t::GeoTrip)
    return Dates.dayofweek(t.pTime) < 6
end

function onlyWeekdays(trips::GeoData)
    mask = falses(length(trips))
    for (i, t) in enumerate(trips)
        if i%10_000 == 0
            @printf("\r%.2f%% trips checked     ",100*i/length(trips))
        end
        mask[i] = onlyWeekdays(t)
    end
    newTrips = trips[mask]
    @printf("\r%.2f%% trips removed\n", 100*(1-length(newTrips)/length(trips)))
    return newTrips
end

"""
    `onlyWeekends`: keep trips that occur on weekends
    - for GeoTrip object: returns boolean
    - for GeoData object: returns filtered GeoData object
"""
function onlyWeekends(t::GeoTrip)
    return Dates.dayofweek(t.pTime) > 5
end

function onlyWeekends(trips::GeoData)
    mask = falses(length(trips))
    for (i, t) in enumerate(trips)
        if i%10_000 == 0
            @printf("\r%.2f%% trips checked     ",100*i/length(trips))
        end
        mask[i] = onlyWeekends(t)
    end
    newTrips = trips[mask]
    @printf("\r%.2f%% trips removed\n", 100*(1-length(newTrips)/length(trips)))
    return newTrips
end

"""
    `inPolygon`: keep trips with pickup and dropoff inside a polygon
    - for one trip: returns boolean
    - for trip list: returns filtered list
"""
function inPolygon(t::GeoTrip, poly::Vector{Tuple{Float32,Float32}})
    return pointInsidePolygon(t.pLon,t.pLat,poly) && pointInsidePolygon(t.dLon,t.dLat,poly)
end

function inPolygon(trips::GeoData, poly::Vector{Tuple{Float32,Float32}})
    mask = falses(length(trips))
    for (i,t) in enumerate(trips)
        if i%10_000 == 0
            @printf("\r%.2f%% trips checked     ",100*i/length(trips))
        end
        mask[i] = inPolygon(t,poly)
    end
    t = trips[mask]
    @printf("\r%.2f%% trips removed\n", 100*(1-length(t)/length(trips)))
    return t
end

"""
    `getPolygon`: given a name, returns the associated geographical polygon
    - For now, just a few so hardcoded, later: perhaps file system
"""
function getPolygon(name::AbstractString)
    if name == "Manhattan"
        return Tuple{Float32,Float32}[(-74.01369f0,40.69977f0), (-74.00597f0,40.702637f0), (-73.99944f0,40.70641f0), (-73.991714f0,40.708492f0), (-73.9761f0,40.71044f0), (-73.96923f0,40.72931f0), (-73.973526f0,40.736073f0), (-73.9615f0,40.75402f0), (-73.941765f0,40.774693f0), (-73.94348f0,40.78223f0), (-73.938156f0,40.78535f0), (-73.93593f0,40.79029f0), (-73.928894f0,40.79432f0), (-73.92872f0,40.803024f0), (-73.93318f0,40.80744f0), (-73.9349f0,40.833942f0), (-73.92134f0,40.85745f0), (-73.91893f0,40.858356f0), (-73.913956f0,40.863678f0), (-73.909706f0,40.872345f0), (-73.91829f0,40.875168f0), (-73.92648f0,40.879192f0), (-73.93344f0,40.87244f0), (-73.933525f0,40.86793f0), (-73.943436f0,40.853584f0), (-73.947945f0,40.85164f0), (-73.94713f0,40.84414f0), (-73.9552f0,40.828682f0), (-73.96091f0,40.8205f0), (-73.97734f0,40.79864f0), (-73.98957f0,40.78077f0), (-73.996994f0,40.770725f0), (-74.00352f0,40.761368f0), (-74.01064f0,40.75103f0), (-74.01532f0,40.719486f0), (-74.01764f0,40.719063f0), (-74.02047f0,40.704067f0)]
    elseif name == "Cambridge"
        return Tuple{Float32,Float32}[(-71.16033f0,42.386353f0),(-71.158035f0,42.386974f0),(-71.15701f0,42.390144f0),(-71.15843f0,42.394585f0),(-71.15651f0,42.39574f0),(-71.1521f0,42.401054f0),(-71.143326f0,42.397266f0),(-71.14316f0,42.39773f0),(-71.14065f0,42.399643f0),(-71.13947f0,42.40012f0),(-71.137825f0,42.400192f0),(-71.13676f0,42.400776f0),(-71.13567f0,42.401546f0),(-71.13397f0,42.40376f0),(-71.1313f0,42.401585f0),(-71.12433f0,42.393604f0),(-71.122856f0,42.39224f0),(-71.12096f0,42.393703f0),(-71.11933f0,42.391636f0),(-71.11706f0,42.38956f0),(-71.11798f0,42.389065f0),(-71.11744f0,42.388367f0),(-71.11818f0,42.388f0),(-71.11549f0,42.38504f0),(-71.11169f0,42.38177f0),(-71.110695f0,42.381683f0),(-71.10785f0,42.379375f0),(-71.108444f0,42.37912f0), (-71.10324f0,42.374348f0),(-71.09753f0,42.374443f0),(-71.09374f0,42.374264f0),(-71.0898f0,42.373955f0),(-71.08607f0,42.37305f0),(-71.08542f0,42.374054f0),(-71.08495f0,42.37346f0),(-71.08441f0,42.373302f0),(-71.08241f0,42.37341f0),(-71.0812f0,42.374058f0),(-71.08036f0,42.37412f0),(-71.079956f0,42.37381f0),(-71.07899f0,42.374172f0),(-71.077545f0,42.373566f0),(-71.075966f0,42.374027f0),(-71.07429f0,42.372536f0),(-71.07267f0,42.373222f0),(-71.07242f0,42.37263f0),(-71.066795f0,42.37184f0),(-71.06405f0,42.369076f0),(-71.06958f0,42.368954f0),(-71.07437f0,42.363743f0),(-71.07485f0,42.3619f0),(-71.07669f0,42.359463f0),(-71.0778f0,42.35864f0),(-71.07809f0,42.357986f0),(-71.08976f0,42.35472f0),(-71.09736f0,42.353386f0),(-71.107574f0,42.35257f0),(-71.10974f0,42.352524f0),(-71.11314f0,42.35298f0),(-71.11584f0,42.354397f0),(-71.117035f0,42.355972f0),(-71.116554f0,42.35941f0),(-71.11706f0,42.3666f0),(-71.11738f0,42.367508f0),(-71.11807f0,42.368336f0),(-71.11947f0,42.36878f0),(-71.12224f0,42.3687f0),(-71.12363f0,42.36908f0),(-71.12459f0,42.369778f0),(-71.12699f0,42.37273f0),(-71.128426f0,42.37367f0),(-71.12929f0,42.37393f0),(-71.130936f0,42.373913f0),(-71.13242f0,42.373425f0),(-71.133194f0,42.37248f0),(-71.13313f0,42.371937f0),(-71.13185f0,42.370285f0),(-71.13181f0,42.369446f0),(-71.132256f0,42.369034f0),(-71.13574f0,42.367092f0),(-71.138954f0,42.364475f0),(-71.14051f0,42.364162f0),(-71.14258f0,42.36476f0),(-71.1436f0,42.364784f0),(-71.14434f0,42.366383f0),(-71.14249f0,42.36688f0),(-71.140976f0,42.369232f0),(-71.14038f0,42.370716f0),(-71.14004f0,42.373566f0),(-71.14396f0,42.374046f0),(-71.14367f0,42.37525f0),(-71.14969f0,42.374554f0),(-71.15439f0,42.37506f0),(-71.15627f0,42.38064f0),(-71.15507f0,42.38021f0),(-71.1552f0,42.38119f0),(-71.15759f0,42.38206f0),(-71.16033f0,42.386353f0)]
    else
        error("Polygon $name not found")
    end
end


text:
module MusicModels

using Flux
using CUDA

device = Flux.cpu

if CUDA.has_cuda()
    device = Flux.gpu
    println("GPU is enabled.")
    CUDA.allowscalar(false)
end

include("pretrain.jl")
include("datasets.jl")
include("MusicTransformer/MusicTransformer.jl")
include("PerformanceRNN/PerformanceRNN.jl")

using .PerformanceRNN
using .MusicTransformer

export generate

function call_registers()
    register_configs(pretrained_configs)
    register_datasets(dataset_configs)
end

@init call_registers()

end


text:
@testset "adt" begin
    @testset "adt List" begin
        # 1. define List
        @data List{T} begin
            Nil()
            Cons(head :: T, tail :: List{T})
        end
        # 2. define interfaces
        len(xs::List{T}) where T = @match xs begin
            Nil{T}() => 0
            Cons{T}(_, tail) => 1 + len(tail)
        end

        @test len(Nil{Any}()) == 0
        xs = Cons(3,Cons(2, Cons(1, Nil{Int}())))
        @test len(xs) == 3
    end
    @testset "adt Arith" begin
        # 1. define Arith
        @data Arith begin
            Num(v :: Int)
            Minus(fst :: Arith, snd :: Arith)
            Add(fst :: Arith, snd :: Arith)
            Mult(fst :: Arith, snd :: Arith)
            Divide(fst :: Arith, snd :: Arith)
        end
        # 2. define interfaces
        function eval_arith(arith :: Arith)
            @match arith begin
                Num(v)       => v
                Add(fst, snd) => eval_arith(fst) + eval_arith(snd)
                Minus(fst, snd) => eval_arith(fst) - eval_arith(snd)
                Mult(fst, snd)   => eval_arith(fst) * eval_arith(snd)
                Divide(fst, snd) => eval_arith(fst) / eval_arith(snd)
            end
        end
        Number = Num
        @test eval_arith(
            Add(Number(1),
                Minus(Number(2),
                    Divide(Number(20),
                            Mult(Number(2),
                                Number(5)))))) == 1
    end

end

@testset "case" begin
    @data CD begin
        D(a, b)
        C{T} :: (a :: Int, b::T) => CD
    end
    @data A begin
        E()
    end
    @test E <: A
    @test fieldnames(D) == (:a, :b)
    @test_throws MethodError C(3.0, :abc)
end

module ADummy
    using MLStyle
end

module BDummy
    using MLStyle
end

@testset "share data with several modules" begin
    @data visible in [ADummy] SSS begin
        SSS_1(Int)
    end
    ADummy.eval(:(SSS_1 = $SSS_1; SSS = $SSS))

    @test ADummy.eval(quote
        @match SSS_1(2) begin
            SSS_1(_) => :ok
        end
    end) == :ok

    BDummy.eval(:(SSS_1 = $SSS_1; SSS = $SSS))

    @test_skip BDummy.eval(quote
        @match SSS_1(2) begin
            SSS_1(_) => :ok
        end
    end)
end


text:
module Inzidenz

end

max_stars_repo_name: EightSQ/inzidenzkarte

max_issues_repo_name: EightSQ/inzidenzkarte

max_forks_repo_name: EightSQ/inzidenzkarte


text:
# consume upto 50% of available memory, across all workers, and across inputs, U and P matrices
def_cache() = max(MAX_BLK_BYTES, floor(Int, Base.Sys.free_memory()/3/2/nworkers()))

zero{T<:AbstractString}(::Type{T}) = convert(T, "")

# TODO: optimize save/load instead of a blind serialize call
function save(model, filename::AbstractString)
    open(filename, "w") do f
        serialize(f, model)
    end
    nothing
end

function load(filename::AbstractString)
    open(filename, "r") do f
        deserialize(f)
    end
end

##
# Types of file specifications

# DLM/CSV file (input only)
type DlmFile <: FileSpec
    name::AbstractString
    dlm::Char
    header::Bool
    quotes::Bool

    function DlmFile(name::AbstractString; dlm::Char=Base.DataFmt.invalid_dlm(Char), header::Bool=false, quotes::Bool=true)
        new(name, dlm, header, quotes)
    end
end

function read_input(fspec::DlmFile)
    # read file and skip the header
    F = readdlm(fspec.name, fspec.dlm, header=fspec.header, quotes=fspec.quotes)
    fspec.header ? F[1] : F
end

# MAT file (input only)
type MatFile <: FileSpec
    filename::AbstractString
    entryname::AbstractString
end
read_input(fspec::MatFile) = read(matopen(fspec.filename), fspec.entryname)

# sparse matrix - in memory (input only)
type SparseMat <: FileSpec
    S::Union{RatingMatrix, SharedRatingMatrix}
end
read_input(fspec::SparseMat) = fspec.S

# sparse matrix blobs (input only)
type SparseBlobs <: FileSpec
    name::AbstractString
    maxcache::Int
    function SparseBlobs(name::AbstractString; maxcache::Int=def_cache())
        new(name, maxcache)
    end
end
read_input(fspec::SparseBlobs) = SparseMatBlobs(fspec.name; maxcache=fspec.maxcache)

# dense matrix blobs (input and output)
type DenseBlobs <: FileSpec
    name::AbstractString
    maxcache::Int
    function DenseBlobs(name::AbstractString; maxcache::Int=def_cache())
        new(name, maxcache)
    end
end
read_input(fspec::DenseBlobs) = DenseMatBlobs(fspec.name; maxcache=fspec.maxcache)

const MAX_BLK_BYTES = 128*1000*1000 #128MB
function _max_items(T::Type, sz::Tuple)
    m,n = sz
    ceil(Int, MAX_BLK_BYTES/sizeof(T)/m)
end
function create{T}(fspec::DenseBlobs, ::Type{T}, sz::Tuple, init::Function, max_items::Int=_max_items(T,sz))
    @logmsg("creating densematarray")
    isdir(fspec.name) || mkdir(fspec.dir)
    m,n = sz
    dm = DenseMatBlobs(T, fspec.name; maxcache=fspec.maxcache)

    startidx = 1
    while startidx <= n
        idxrange = startidx:min(n, startidx + max_items)
        blobsz = (m,length(idxrange))
        M = init(T, blobsz...)
        @logmsg("idxrange: $idxrange, sz: $(size(M))")
        append!(dm, M)
        startidx = last(idxrange) + 1
    end
    dm
end


text:
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]

parr = [
    ["2D state", :η, 12, 12, 0, 12],
    ["2D state", "U[1]", 12, 12, 0, 12],
    ["2D state", "U[2]", 0, 0, 0, 0],
]

explicit = [
    [
        "2D state",
        :η,
        -8.52722969951589915283e-01,
        8.52846676313531282254e-01,
        -2.49578135935735214742e-16,
        6.03454239990563690021e-01,
    ],
    [
        "2D state",
        "U[1]",
        -3.15431401945821825450e+01,
        3.15431401945818628008e+01,
        6.11504145930918957291e-15,
        2.24273815174625497093e+01,
    ],
    [
        "2D state",
        "U[2]",
        -7.62224398365580242501e-13,
        9.72156930292624284356e-13,
        1.39269607441935025982e-14,
        1.95606703846656748360e-13,
    ],
]

refVals = (explicit = (explicit, parr),)

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
####################################
#
# Recursive dynamic extension of SLiDE model
#
####################################

using SLiDE
using CSV
using JuMP
using Complementarity
using DataFrames

############
# LOAD DATA
############

#SLiDE data needs to be built or point to pre-existing build directory
#can pass a name (d, set) = build_data("name_of_build_directory")
!(@isdefined(d_in) && @isdefined(set_in)) && ((d_in, set_in) = build("state_model"))
d = copy(d_in)
set = copy(set_in)

#Specify benchmark year - this is the first solve year where the benchmark replicated
bmkyr = 2016

#Load slide data and time horizon to produce model data and appropriate time-indexed subsets
(sld, set, idx) = model_input(bmkyr, d, set)

########## Model ##########
cge = MCPModel();

##############
# SETS
##############

# Set description
#set[:s] -> sectors
#set[:g] -> goods
#set[:r] -> regions
#set[:m] -> margins
#set[:gm] -> goods_margins

##############
# PARAMETERS
##############

#benchmark values
@NLparameter(cge, ys0[r in set[:r], s in set[:s], g in set[:g]] == sld[:ys0][r,s,g]); # Sectoral supply
@NLparameter(cge, id0[r in set[:r], s in set[:s], g in set[:g]] == sld[:id0][r,s,g]); # Intermediate demand
@NLparameter(cge, ld0[r in set[:r], s in set[:s]] == sld[:ld0][r,s]); # Labor demand
@NLparameter(cge, kd0[r in set[:r], s in set[:s]] == sld[:kd0][r,s]); # Capital demand
@NLparameter(cge, ty0[r in set[:r], s in set[:s]] == sld[:ty0][r,s]); # Production tax (benchmark)
@NLparameter(cge, m0[r in set[:r], g in set[:g]] == sld[:m0][r,g]); # Imports
@NLparameter(cge, x0[r in set[:r], g in set[:g]] == sld[:x0][r,g]); # Exports of goods and services
@NLparameter(cge, rx0[r in set[:r], g in set[:g]] == sld[:rx0][r,g]); # Re-exports of goods and services
@NLparameter(cge, md0[r in set[:r], m in set[:m], g in set[:g]] == sld[:md0][r,m,g]); # Total margin demand
@NLparameter(cge, nm0[r in set[:r], g in set[:g], m in set[:m]] == sld[:nm0][r,g,m]); # Margin demand from national market
@NLparameter(cge, dm0[r in set[:r], g in set[:g], m in set[:m]] == sld[:dm0][r,g,m]); # Margin supply from local market
@NLparameter(cge, s0[r in set[:r], g in set[:g]] == sld[:s0][r,g]); # Aggregate supply
@NLparameter(cge, a0[r in set[:r], g in set[:g]] == sld[:a0][r,g]); # Armington supply
@NLparameter(cge, ta0[r in set[:r], g in set[:g]] == sld[:ta0][r,g]); # Tax net of subsidy rate on intermediate demand
@NLparameter(cge, tm0[r in set[:r], g in set[:g]] == sld[:tm0][r,g]); # Import tariff
@NLparameter(cge, cd0[r in set[:r], g in set[:g]] == sld[:cd0][r,g]); # Final demand
@NLparameter(cge, c0[r in set[:r]] == sld[:c0][r]); # Aggregate final demand
@NLparameter(cge, yh0[r in set[:r], g in set[:g]] == sld[:yh0][r,g]); #Household production
@NLparameter(cge, bopdef0[r in set[:r]] == sld[:bopdef0][r]); #Balance of payments
@NLparameter(cge, hhadj[r in set[:r]] == sld[:hhadj][r]); # Household adjustment
@NLparameter(cge, g0[r in set[:r], g in set[:g]] == sld[:g0][r,g]); # Government demand
@NLparameter(cge, i0[r in set[:r], g in set[:g]] == sld[:i0][r,g]); # Investment demand
@NLparameter(cge, xn0[r in set[:r], g in set[:g]] == sld[:xn0][r,g]); # Regional supply to national market
@NLparameter(cge, xd0[r in set[:r], g in set[:g]] == sld[:xd0][r,g]); # Regional supply to local market
@NLparameter(cge, dd0[r in set[:r], g in set[:g]] == sld[:dd0][r,g]); # Regional demand to local market
@NLparameter(cge, nd0[r in set[:r], g in set[:g]] == sld[:nd0][r,g]); # Regional demand to national market

#counterfactual taxes
@NLparameter(cge, ty[r in set[:r], s in set[:s]] == sld[:ty0][r,s]); # output tax
@NLparameter(cge, ta[r in set[:r], g in set[:g]] == sld[:ta0][r,g]); # Armington tax
@NLparameter(cge, tm[r in set[:r], g in set[:g]] == sld[:tm0][r,g]); # import tariff

# -- Major Assumptions --
# Temporal/Dynamic modifications

@NLparameter(cge, ir == 0.05); # interest rate
@NLparameter(cge, gr == 0.02); # growth rate --- try sector and set[:r] specific
@NLparameter(cge, dr == 0.05); # capital depreciation rate
@NLparameter(cge, thetax == 0.3); # extant production share

#new capital endowment
@NLparameter(cge, ks_n[r in set[:r], s in set[:s]] ==
    value(kd0[r, s])  * (value(dr)+value(gr)) / (1 + value(gr)) );

# mutable old capital endowment
@NLparameter(cge, ks_s[r in set[:r], s in set[:s]] ==
    value(kd0[r, s]) * (1 - value(thetax)) - value(ks_n[r,s]) );

# Mutable total capital endowment - Non-extant capital
@NLparameter(cge, ks_m[r in set[:r]] == sum(value(kd0[r,s]) * (1-value(thetax)) for s in set[:s]));

# Extant capital endowment
@NLparameter(cge, ks_x[r in set[:r], s in set[:s]] ==
    value(kd0[r, s]) * value(thetax) );

# Labor endowment
@NLparameter(cge, le0[r in set[:r], s in set[:s]] == value(ld0[r,s]));

# Benchmark investment supply
@NLparameter(cge, inv0[r in set[:r]] == sum(value(i0[r,g]) for g in set[:g]));

# Benchmark welfare index
@NLparameter(cge, w0[r in set[:r]] == value(inv0[r])+value(c0[r]));

# benchmark value share parameters
@NLparameter(cge, alpha_kl[r in set[:r], s in set[:s]] == ensurefinite(value(ld0[r,s]) / (value(ld0[r,s]) + value(kd0[r,s]))));
@NLparameter(cge, alpha_x[r in set[:r], g in set[:g]] == ensurefinite((value(x0[r,g]) - value(rx0[r,g])) / value(s0[r,g])));
@NLparameter(cge, alpha_d[r in set[:r], g in set[:g]] == ensurefinite(value(xd0[r,g]) / value(s0[r,g])));
@NLparameter(cge, alpha_n[r in set[:r], g in set[:g]] == ensurefinite(value(xn0[r,g]) / value(s0[r,g])));
@NLparameter(cge, theta_n[r in set[:r], g in set[:g]] == ensurefinite(value(nd0[r,g]) / (value(nd0[r,g]) - value(dd0[r,g]))));
@NLparameter(cge, theta_m[r in set[:r], g in set[:g]] == ensurefinite((1+value(tm0[r,g])) * value(m0[r,g]) / (value(nd0[r,g]) + value(dd0[r,g]) + (1 + value(tm0[r,g])) * value(m0[r,g]))));

@NLparameter(cge, theta_inv[r in set[:r], g in set[:g]] == value(i0[r,g]) / value(inv0[r])); # Intermediate input share of investment output

#Substitution and transformation elasticities
@NLparameter(cge, es_va[r in set[:r], s in set[:s]] == SUB_ELAST[:va]); # value-added nest - substitution elasticity
@NLparameter(cge, es_y[r in set[:r], s in set[:s]]  == SUB_ELAST[:y]); # Top-level Y nest (VA,M) - substitution elasticity
@NLparameter(cge, es_m[r in set[:r], s in set[:s]]  == SUB_ELAST[:m]); # Materials nest - substitution elasticity
@NLparameter(cge, et_x[r in set[:r], g in set[:g]]    == TRANS_ELAST[:x]); # Disposition, distribute regional supply to local, national, export - transformation elasticity
@NLparameter(cge, es_a[r in set[:r], g in set[:g]]    == SUB_ELAST[:a]); # Top-level A nest for aggregate demand (set[:m], set[:g]) - substitution elasticity
@NLparameter(cge, es_mar[r in set[:r], g in set[:g]]  == SUB_ELAST[:mar]); # Margin supply - substitution elasticity
@NLparameter(cge, es_d[r in set[:r], g in set[:g]]    == SUB_ELAST[:d]); # Domestic demand aggregation nest (intranational) - substitution elasticity
@NLparameter(cge, es_f[r in set[:r], g in set[:g]]    == SUB_ELAST[:f]); # Domestic and foreign demand aggregation nest (international) - substitution elasticity

@NLparameter(cge, es_inv[r in set[:r]] == 5); # Investment production - substitution elasticity

################
# VARIABLES
################

# specify lower bound
lo = 0.0

# sectors
#@variable(cge, Y[(r, s) in set[:Y]] >= lo, start = 1);
@variable(cge, X[(r, g) in set[:X]] >= lo, start = 1); # Exports
@variable(cge, A[(r, g) in set[:A]] >= lo, start = 1); # Armington
@variable(cge, C[r in set[:r]] >= lo, start = 1); # Consumption
@variable(cge, MS[r in set[:r], m in set[:m]] >= lo, start = 1); # Margin Supply

#commodities:
@variable(cge, PA[(r, g) in set[:PA]] >= lo, start = 1); # Regional market (input)
@variable(cge, PY[(r, g) in set[:PY]] >= lo, start = 1); # Regional market (output)
@variable(cge, PD[(r, g) in set[:PD]] >= lo, start = 1); # Local market price
@variable(cge, PN[g in set[:g]] >= lo, start =1); # National market
@variable(cge, PL[r in set[:r]] >= lo, start = 1); # Wage rate
#@variable(cge, PK[(r, s) in set[:PK]] >= lo, start =1); # Rental rate of capital ###
@variable(cge, PM[r in set[:r], m in set[:m]] >= lo, start =1); # Margin price
@variable(cge, PC[r in set[:r]] >= lo, start = 1); # Consumer price index #####
@variable(cge, PFX >= lo, start = 1); # Foreign exchange

#consumer:
@variable(cge,RA[r in set[:r]]>=lo,start = value(w0[r])) ;

#--- recursive dynamic variable declaration ---
# sectors
@variable(cge,YM[(r,s) in set[:Y]] >= lo, start = (1-value(thetax))); # Mutable production index - replaces Y
@variable(cge,YX[(r,s) in set[:Y]] >= lo, start = value(thetax)); # Extant production index
@variable(cge,INV[r in set[:r]] >= lo, start = 1); # Investment
@variable(cge,W[r in set[:r]] >= lo, start = 1); # Welfare index

# commodities
@variable(cge,RKX[(r,s) in set[:PK]] >= lo, start = 1); # Return to extant capital
@variable(cge,RK[r in set[:r]] >= lo, start = 1); # Return to regional capital
@variable(cge,PINV[r in set[:r]] >= lo, start = 1); # Investment price index
@variable(cge,PW[r in set[:r]] >= lo, start = 1); # Welfare price index

# Reporting variables
@variable(cge,DKM[(r,s) in set[:PK]] >= lo, start = start_value(YM[(r,s)]) * value(kd0[r,s]));

###############################
# -- PLACEHOLDER VARIABLES --
###############################

#----------
### Recursive Model expressions

#Cobb-douglas for mutable/new
@NLexpression(cge, CVAym[r in set[:r], s in set[:s]],
    PL[r]^alpha_kl[r,s] * RK[r]^(1-alpha_kl[r,s]));

#demand for labor in VA
@NLexpression(cge,ALym[r in set[:r], s in set[:s]],
    ld0[r,s] * CVAym[r,s] / PL[r]);

#demand for capital in VA
@NLexpression(cge,AKym[r in set[:r],s in set[:s]],
    kd0[r,s] * CVAym[r,s] / RK[r]);

@NLexpression(cge,DINV[r in set[:r], g in set[:g]],
    (i0[r,g] * (PINV[r]/(haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0)^es_inv[r])));

#----------

#CES function for output demand - including
# exports (absent of re-exports) times the price for foreign exchange,
# region's supply to national market times the national market price
# regional supply to local market times domestic price
@NLexpression(cge,RX[r in set[:r],g in set[:g]],
  (alpha_x[r,g]*PFX^5+alpha_n[r,g]*PN[g]^5+alpha_d[r,g]*(haskey(PD.lookup[1], (r, g)) ? PD[(r, g)] : 1.0)^5)^(1/5) );

#demand for exports via demand function
@NLexpression(cge,AX[r in set[:r],g in set[:g]], (x0[r,g] - rx0[r,g])*(PFX/RX[r,g])^4 );

#demand for contribution to national market
@NLexpression(cge,AN[r in set[:r],g in set[:g]], xn0[r,g]*(PN[g]/(RX[r,g]))^4 );

#demand for regionals supply to local market
@NLexpression(cge,AD[r in set[:r],g in set[:g]],
  xd0[r,g] * ((haskey(PD.lookup[1], (r, g)) ? PD[(r, g)] : 1.0) / (RX[r,g]))^4 );

  ###

# CES function for tradeoff between national and domestic market
@NLexpression(cge,CDN[r in set[:r],g in set[:g]],
  (theta_n[r,g]*PN[g]^(1-2)+(1-theta_n[r,g])*(haskey(PD.lookup[1], (r, g)) ? PD[(r, g)] : 1.0)^(1-2))^(1/(1-2)) );

# CES function for tradeoff between domestic consumption and foreign exports
# recall tm in the import tariff thus tm / tm0 is the relative change in import tariff rates
@NLexpression(cge,CDM[r in set[:r],g in set[:g]],
  ((1-theta_m[r,g])*CDN[r,g]^(1-4)+theta_m[r,g]*(PFX*(1+tm[r,g])/(1+tm0[r,g]))^(1-4))^(1/(1-4)) );

# set[:r] demand from the national market <- note nesting of CDN in CDM
@NLexpression(cge,DN[r in set[:r],g in set[:g]],
  nd0[r,g]*(CDN[r,g]/PN[g])^2*(CDM[r,g]/CDN[r,g])^4 );

# region demand from local market <- note nesting of CDN in CDM
@NLexpression(cge,DD[r in set[:r],g in set[:g]],
  dd0[r,g]*(CDN[r,g]/(haskey(PD.lookup[1], (r, g)) ? PD[(r, g)] : 1.0))^2*(CDM[r,g]/CDN[r,g])^4 );

# import demand
@NLexpression(cge,MD[r in set[:r],g in set[:g]],
  m0[r,g]*(CDM[r,g]*(1+tm[r,g])/(PFX*(1+tm0[r,g])))^4 );

# final demand
@NLexpression(cge,CD[r in set[:r],g in set[:g]],
  cd0[r,g]*PC[r] / (haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0) );

###############################
# -- Zero Profit Conditions --
###############################

#----------
#Recursive  --- update to Y
@mapping(cge,profit_ym[(r, s) in set[:Y]],
# cost of intermediate demand
    sum((haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0) * id0[r,g,s] for g in set[:g])
# cost of labor inputs
    + PL[r] * ALym[r,s]
# cost of capital inputs
    + RK[r] * AKym[r,s]
    -
# revenue from sectoral supply (take note of r/s/g indices on ys0)
    sum((haskey(PY.lookup[1], (r, g)) ? PY[(r, g)] : 1.0)  * ys0[r,s,g] for g in set[:g]) * (1-ty[r,s])
);

@mapping(cge,profit_yx[(r, s) in set[:Y]],
# cost of intermediate demand
    sum((haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0) * id0[r,g,s] for g in set[:g])
# cost of labor inputs
    + PL[r] * ld0[r,s]
# cost of capital inputs
    + (haskey(RKX.lookup[1], (r, s)) ? RKX[(r,s)] : 1.0) * kd0[r,s]
    -
# revenue from sectoral supply (take note of r/s/g indices on ys0)
    sum((haskey(PY.lookup[1], (r, g)) ? PY[(r, g)] : 1.0)  * ys0[r,s,g] for g in set[:g]) * (1-ty[r,s])
);

#----------

@mapping(cge,profit_x[(r, g) in set[:X]],
# output 'cost' from aggregate supply
    (haskey(PY.lookup[1], (r, g)) ? PY[(r, g)] : 1.0) * s0[r,g]
    - (
# revenues from foreign exchange
        PFX * AX[r,g]
# revenues from national market
        + PN[g] * AN[r,g]
# revenues from domestic market
        + (haskey(PD.lookup[1], (r, g)) ? PD[(r, g)] : 1.0) * AD[r,g]
    )
);

@mapping(cge,profit_a[(r, g) in set[:A]],
# costs from national market
    PN[g] * DN[r,g]
# costs from domestic market
    + (haskey(PD.lookup[1], (r, g)) ? PD[(r, g)] : 1.0) * DD[r,g]
# costs from imports, including import tariff
    + PFX * (1+tm[r,g]) * MD[r,g]
# costs of margin demand
    + sum(PM[r,m] * md0[r,m,g] for m in set[:m])
    - (
# revenues from regional market based on armington supply
       (haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0) * (1-ta[r,g]) * a0[r,g]
# revenues from re-exports
        + PFX * rx0[r,g]
    )
);

@mapping(cge, profit_c[r in set[:r]],
# costs of inputs - computed as final demand times regional market prices
    sum((haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0) * CD[r,g] for g in set[:g])
    -
# revenues/benefit computed as CPI * reference consumption
    PC[r] * c0[r]
);

@mapping(cge, profit_inv[r in set[:r]],
# inputs to investment
    sum((haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0) * DINV[r,g] for g in set[:g])
    -
# Investment
    PINV[r] * inv0[r]
);

@mapping(cge, profit_w[r in set[:r]],
# inputs to welfare index
    PINV[r]*inv0[r] + PC[r]*c0[r]
    -
# Welfare
    PW[r]*w0[r]
);

@mapping(cge,profit_ms[r in set[:r], m in set[:m]],
# provision of set[:m] to national market
    sum(PN[gm]   * nm0[r,gm,m] for gm in set[:gm])
# provision of set[:m] to domestic market
    + sum((haskey(PD.lookup[1], (r, gm)) ? PD[(r, gm)] : 1.0) * dm0[r,gm,m] for gm in set[:gm])
    -
# total margin demand
    PM[r,m] * sum(md0[r,m,gm] for gm in set[:gm])
);

###################################
# -- Market Clearing Conditions --
###################################

#----------
#Recursive dynamics mkt clearance

@mapping(cge,market_rk[r in set[:r]],
# mutable capital supply
    ks_m[r]
    -
# mutable capital demand
    sum((haskey(YM.lookup[1], (r, s)) ? YM[(r, s)] : 1.0) * AKym[r,s] for s in set[:s])
);

@mapping(cge,market_rkx[(r, s) in set[:PK]],
# extant capital supply
    (ks_x[r,s])
    -
# extant capital demand
    (haskey(YX.lookup[1], (r, s)) ? YX[(r, s)] : 1.) * kd0[r,s]
);

@mapping(cge,market_pa[(r, g) in set[:PA]],
# absorption or supply
    (haskey(A.lookup[1], (r, g)) ? A[(r, g)] : 1.) * a0[r,g]
    - (
# government demand (exogenous)
        g0[r,g]
# demand for investment
        + INV[r]*DINV[r,g]
# final demand
        + C[r] * CD[r,g]
# intermediate demand
        + sum((haskey(YM.lookup[1], (r, s)) ? YM[(r, s)] : 1) * id0[r,g,s] for s in set[:s] if (r,s) in set[:Y])
        + sum((haskey(YX.lookup[1], (r, s)) ? YX[(r, s)] : 1) * id0[r,g,s] for s in set[:s] if (r,s) in set[:Y])
    )
);

@mapping(cge,market_py[(r, g) in set[:PY]],
# sectoral supply
    sum((haskey(YM.lookup[1], (r, s)) ? YM[(r, s)] : 0.0) *ys0[r,s,g] for s in set[:s])
    + sum((haskey(YX.lookup[1], (r, s)) ? YX[(r, s)] : 0.0) *ys0[r,s,g] for s in set[:s])
# household production (exogenous)
    + yh0[r,g]
    -
# aggregate supply (akin to market demand)
    (haskey(X.lookup[1], (r, g)) ? X[(r, g)] : 1) * s0[r,g]
);

@mapping(cge,market_pl[r in set[:r]],
# supply of labor
    sum(le0[r,s] for s in set[:s])
    -
# demand for labor in all set[:s]
    (
        sum((haskey(YM.lookup[1], (r, s)) ? YM[(r, s)] : 1) * ALym[r,s] for s in set[:s])
        + sum((haskey(YX.lookup[1], (r, s)) ? YX[(r, s)] : 1) * ld0[r,s] for s in set[:s])
    )
);
#----------

@mapping(cge,market_pd[(r, g) in set[:PD]],
# aggregate supply
    (haskey(X.lookup[1], (r, g)) ? X[(r, g)] : 1)  * AD[r,g]
    - (
# demand for local market
        (haskey(A.lookup[1], (r, g)) ? A[(r, g)] : 1.) * DD[r,g]
# margin formulation from local market
        + sum(MS[r,m] * dm0[r,g,m] for m in set[:m] if (g in set[:gm] ) )
    )
);

@mapping(cge,market_pn[g in set[:g]],
# supply to the national market
    sum((haskey(X.lookup[1], (r, g)) ? X[(r, g)] : 1)  * AN[r,g] for r in set[:r])
    - (
# demand from the national market
        sum((haskey(A.lookup[1], (r, g)) ? A[(r, g)] : 1.) * DN[r,g] for r in set[:r])
# margin formulation from the national market
        + sum(MS[r,m] * nm0[r,g,m] for r in set[:r] for m in set[:m] if (g in set[:gm]) )
    )
);

@mapping(cge,market_pm[r in set[:r], m in set[:m]],
# margin supply
    MS[r,m] * sum(md0[r,m,gm] for gm in set[:gm])
    -
# margin demand
    sum((haskey(A.lookup[1], (r, g)) ? A[(r, g)] : 1.0) * md0[r,m,g] for g in set[:g])
);

@mapping(cge,market_pc[r in set[:r]],
# Consumption supply
    C[r] * c0[r]
    -
# Consumption demand
    W[r] * c0[r]
);

@mapping(cge, market_pinv[r in set[:r]],
# investment supply
    INV[r]*inv0[r]
    -
# investment demanded
    W[r]*inv0[r]
);

@mapping(cge, market_pw[r in set[:r]],
# welfare supply
    W[r]*w0[r]
    -
# full consumption welfare / utility
    RA[r] / PW[r]
);

@mapping(cge,market_pfx,
# balance of payments (exogenous)
    sum(bopdef0[r] for r in set[:r])
# supply of exports
    + sum((haskey(X.lookup[1], (r, g)) ? X[(r, g)] : 1.0)  * AX[r,g] for r in set[:r] for g in set[:g])
# supply of re-exports
    + sum((haskey(A.lookup[1], (r, g)) ? A[(r, g)] : 1.0) * rx0[r,g] for r in set[:r] for g in set[:g] if (r,g) in set[:A])
    -
# import demand
    sum((haskey(A.lookup[1], (r, g)) ? A[(r, g)] : 1.0) * MD[r,g] for r in set[:r] for g in set[:g] if (r,g) in set[:A])
);

#----------
#Income balance update for recursive dynamics
@mapping(cge,income_ra[r in set[:r]],
# consumption/utility
    RA[r]
    -
    (
# labor income
        PL[r] * sum(le0[r,s] for s in set[:s])
# capital income
        + RK[r] * ks_m[r]
        +sum((haskey(RKX.lookup[1], (r, s)) ? RKX[(r,s)] : 1.0) * ks_x[r,s] for s in set[:s])
# provision of household supply
        + sum( (haskey(PY.lookup[1], (r, g)) ? PY[(r, g)] : 1.0) * yh0[r,g] for g in set[:g])
# revenue or costs of foreign exchange including household adjustment
        + PFX * (bopdef0[r] + hhadj[r])
# government provision
        - sum((haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0) * (g0[r,g]) for g in set[:g])
# import taxes - assumes lumpsum recycling
        + sum((haskey(A.lookup[1], (r, g)) ? A[(r, g)] : 1.) * MD[r,g] * PFX * tm[r,g] for g in set[:g] if (r,g) in set[:A])
# taxes on intermediate demand - assumes lumpsum recycling
        + sum((haskey(A.lookup[1], (r, g)) ? A[(r, g)] : 1.) * a0[r,g]*(haskey(PA.lookup[1], (r, g)) ? PA[(r, g)] : 1.0)*ta[r,g] for g in set[:g] if (r,g) in set[:A])
# production taxes - assumes lumpsum recycling
        + sum( (haskey(YM.lookup[1], (r, s)) ? YM[(r, s)] : 1) * ys0[r,s,g] * ty[r,s] for s in set[:s], g in set[:g])
        + sum( (haskey(YX.lookup[1], (r, s)) ? YX[(r, s)] : 1) * ys0[r,s,g] * ty[r,s] for s in set[:s], g in set[:g])
    )
);

#----------
# Reporting
@mapping(cge, DKMdef[(r,s) in set[:PK]],
    DKM[(r,s)]
    -
    YM[(r,s)]*AKym[r,s]
);

####################################
# -- Complementarity Conditions --
####################################

# define complementarity conditions
# note the pattern of ZPC -> primal variable  &  MCC -> dual variable (price)
#@complementarity(cge,profit_y,Y);
@complementarity(cge,profit_ym,YM);
@complementarity(cge,profit_yx,YX);
@complementarity(cge,profit_x,X);
@complementarity(cge,profit_a,A);
@complementarity(cge,profit_c,C);
@complementarity(cge,profit_ms,MS);
@complementarity(cge,market_pa,PA);
@complementarity(cge,market_py,PY);
@complementarity(cge,market_pd,PD);
@complementarity(cge,market_pn,PN);
@complementarity(cge,market_pl,PL);
#@complementarity(cge,market_pk,PK);
@complementarity(cge,market_rk,RK);
@complementarity(cge,market_rkx,RKX);
@complementarity(cge,market_pm,PM);
@complementarity(cge,market_pc,PC);
@complementarity(cge,market_pfx,PFX);
@complementarity(cge,income_ra,RA);

#----------
#Recursive Dynamics
@complementarity(cge,profit_inv,INV);
@complementarity(cge,profit_w,W);
@complementarity(cge,market_pinv,PINV);
@complementarity(cge,market_pw,PW);

#Reporting
@complementarity(cge,DKMdef,DKM);

####################
# -- Model Solve --
####################

#set up the options for the path solver
PATHSolver.options(convergence_tolerance=1e-6, output=:yes, time_limit=3600, cumulative_iteration_limit=0)

# solve the model
status = solveMCP(cge)

# Pre-loop calculations
ktot_mx = Dict((r,bmkyr) => value(ks_m[r]) + sum(value(ks_x[r,s]) for s in set[:s])
    for r in set[:r])

srv = 1-value(dr)

# Begin loop
# !!!! consider defining empty dictionaries outside of the loop to store values for postproc
# Been tested from 2017 to 2050 in 1 year increments
for t in 2017:2020
    @info("Begin loop prior to $t Solve.")

    # !!!! Should growth rate be in here? Benchmark cannot replicate if not here.
    # new mutable capital
    newkm = Dict((r,t) => (ktot_mx[r,bmkyr]*(1+value(gr)) - ktot_mx[r,bmkyr]*srv) * result_value(INV[r])
        for r in set[:r]);

    # mutable capital demand from previous solve
    dkml = Dict((r,s) => (haskey(DKM.lookup[1], (r,s)) ? result_value(DKM[(r,s)]) : 0.0)
        for r in set[:r], s in set[:s]);

    # Total mutable/putty capital
    # !!!! Should old putty be frozen?
    # ktot_m = Dict((r,t) => newkm[r,t] + value(ks_m[r])*srv - sum(value(thetax)*dkml[r,s]*srv for s in set[:s])
    #     for r in set[:r]);
    ktot_m = Dict((r,t) => newkm[r,t] + value(ks_m[r])*srv
        for r in set[:r]);

    for r in set[:r]
        set_value(ks_m[r], ktot_m[r,t]);
    end

    yxl = Dict((r,s) => (haskey(YX.lookup[1], (r,s)) ? result_value(YX[(r,s)]) : 0.0)
        for r in set[:r], s in set[:s]);

    # Total extant/clay capital
    # !!!! Should old clay increase by frozen putty each period?
    # ktot_x = Dict((r,s,t) => srv*yxl[r,s]*value(kd0[r,s]) + value(thetax)*dkml[r,s]*srv
    #     for r in set[:r], s in set[:s]);
    ktot_x = Dict((r,s,t) => srv*yxl[r,s]*value(kd0[r,s])
        for r in set[:r], s in set[:s]);

    for r in set[:r], s in set[:s]
        set_value(ks_x[r,s], ktot_x[r,s,t]);
    end

    for r in set[:r], s in set[:s]
    #update labor endowments --- separate parameters for labor endowments versus demand
        set_value(le0[r,s], (1 + value(gr)) * value(le0[r,s]));
    end

    #update balance of payments and household adjustment
    for r in set[:r]
        set_value(bopdef0[r], (1 + value(gr)) * value(bopdef0[r]));
        set_value(hhadj[r], (1 + value(gr)) * value(hhadj[r]));
    end

    #update government and household production
    for r in set[:r], g in set[:g]
        set_value(g0[r,g], (1 + value(gr)) * value(g0[r,g]));
        set_value(yh0[r,g], (1 + value(gr)) * value(yh0[r,g]));
    end

    #update all model variable start values to previous period solution value
    set_start_value.(all_variables(cge), result_value.(all_variables(cge)));

    #update consumption start value
    for r in set[:r]
        set_start_value(C[r], result_value(C[r])*(1+value(gr)));
        set_start_value(INV[r], result_value(INV[r])*(1+value(gr)));
        set_start_value(W[r], result_value(W[r])*(1+value(gr)));
        set_start_value(RA[r], start_value(W[r])*value(w0[r]));
    end

    #update exports start value
    for (r,g) in set[:X]
        set_start_value(X[(r,g)], result_value(X[(r,g)])*(1+value(gr)));
    end

    #update armington start value
    for (r,g) in set[:A]
        set_start_value(A[(r,g)], result_value(A[(r,g)])*(1+value(gr)));
    end

    #update margin supply start value
    for r in set[:r], m in set[:m]
        set_start_value(MS[r,m], result_value(MS[r,m])*(1+value(gr)));
    end

    #update output variable start values
    # !!!! Eventually you could get so much clay, that YM level value is negative (Divide by kd0?)
    for (r,s) in set[:Y]
        set_start_value(YX[(r,s)], value(ks_x[r,s])/value(kd0[r,s])); # !!!! divide by kd0? Sensitive to thetax and dr
        set_start_value(YM[(r,s)], (start_value(C[r]) - start_value(YX[(r,s)])));
        set_start_value(DKM[(r,s)], start_value(YM[(r,s)])*value(kd0[r,s]));
    end

    #update value shares
    for r in set[:r], s in set[:s]
        set_value(alpha_kl[r,s], ensurefinite(value(ld0[r,s])/(value(ld0[r,s]) + value(kd0[r,s]))));
    end

    #update value shares
    for r in set[:r], g in set[:g]
        set_value(alpha_x[r,g], ensurefinite((value(x0[r, g]) - value(rx0[r, g])) / value(s0[r, g])));
        set_value(alpha_d[r,g], ensurefinite((value(xd0[r,g])) / value(s0[r, g])));
        set_value(alpha_n[r,g], ensurefinite(value(xn0[r,g]) / (value(s0[r, g]))));
        set_value(theta_n[r,g], ensurefinite(value(nd0[r, g]) / (value(nd0[r, g]) + value(dd0[r, g]))));
        set_value(theta_m[r,g], ensurefinite((1+value(tm0[r, g])) * value(m0[r, g]) / (value(nd0[r, g]) + value(dd0[r, g]) + (1 + value(tm0[r, g])) * value(m0[r, g]))));
    end

    #set up the options for the path solver
    #PATHSolver.options(convergence_tolerance=1e-6, output=:yes, time_limit=3600, cumulative_iteration_limit=10000)
    PATHSolver.options(convergence_tolerance=1e-6, output=:yes, time_limit=3600, cumulative_iteration_limit=0)

    # solve next period
    status = solveMCP(cge)
end

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
using AlgebraPDF
using Test

g1 = FGauss((μ1=1.1,σ1=0.1))
g2 = FGauss((μ2=2.2,σ2=0.2))
g3 = FGauss((μ3=3.3,σ3=0.3))
d = AlgebraPDF.FSum([g1,g2,g3], (α1=0.1,α2=0.2,α3=0.3))
d2 = AlgebraPDF.FSum([g1,g2,g3], Ext(α1=0.1,α2=0.2,α3=0.3)) |> x->fixpar(x, :α1)

@testset "Sum of regular functions" begin
    @test length(d) == 3
    @test abs(func(d, 1.5)) < 0.01
    @test abs(func(d, 1.5; p=pars(d)+(μ1=1.5,))) > 0.01
    @test func(d, 1.5; p=pars(d)+(α1=1500,)) > 0.01
    #
    @test abs(d(1.1)-d(2.2)) < 0.01
    @test abs(d(2.2)-d(3.3)) < 0.01

    keys(fixedpars(d2)) == (:α1,)
    @test d2(1.1) == d(1.1)
end
#

ng1 = Normalized(g1, (1,4))
ng2 = Normalized(g2, (1,4))
ng3 = Normalized(g3, (1,4))
nd = FSum([ng1,ng2,ng3], (α1=0.1,α2=0.2,α3=0.3))
#
@testset "Sum of normalized functions" begin

    @test typeof(d) <: FSum{T} where T<:AbstractFunctionWithParameters
    @test typeof(nd) <: FSum{T} where T<:AbstractPDF

    @test func(nd, 1.1) isa Number
    @test func(nd, [1.1,1,1])  isa Vector
    @test func(nd, 1:0.1:2)  isa Vector

    @test lims(nd) == lims(ng1)

    @test nd(1.1) > d(1.1) # due to normalization
end

@testset "Addition: FSum + FSum" begin
    s12 = FSum([g1, g2], (α1=1.1, α2=2.2))
    s3 = FSum([g3], (α3=3.3,))
    s123 = FSum([g1, g2, g3], (α1=1.1, α2=2.2, α3=3.3))

    @test s12 + s3 == s123

    s12_e = FSum([g1, g2], Ext(α1=1.1, α2=2.2))
    s3_e = FSum([g3], Ext(α3=3.3,))
    s123_e = FSum([g1, g2, g3], Ext(α1=1.1, α2=2.2, α3=3.3))

    @test s12_e + s3_e == s123_e
    @test s12_e + s3 == s123_e
    @test s12 + s3_e == s123_e
end

@testset "Multiplication: c * FSum, FSum * c " begin
    @test d == (α1=0.1,) * g1 + (α2=0.2,) * g2 + (α3=0.3,) * g3
    x2 = (Ext(α1=0.1,) * g1 + (α2=0.2,) * g2 + (α3=0.3,) * g3) |> d->fixpar(d,:α1)
    @test typeof(d2) == typeof(x2) && d2.fs == x2.fs && d2.αs == x2.αs
    @test nd == (α1=0.1,) * ng1 + (α2=0.2,) * ng2 + (α3=0.3,) * ng3
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

matrixcompletion_insts = [
    [(k, d) for d in vcat(2, 5:5:max_d)] # includes compile run
    for (k, max_d) in ((10, 45), (20, 30))
    ]

insts = OrderedDict()
insts["nat"] = (nothing, matrixcompletion_insts)
insts["extEP"] = (:ExpPSD, matrixcompletion_insts)
insts["extSEP"] = (:SOCExpPSD, matrixcompletion_insts)
return (MatrixCompletionJuMP, insts)


text:
# utility functions

function make_blocks(nobs, nblocks)
    maxbs, reminder = divrem(nobs, nblocks)

    res = fill(maxbs, nblocks)
    if reminder > 0
        res[1:reminder] .= maxbs + 1
    end
    res
end

# the truncated newton method for matrix inversion
# adapted from https://en.wikipedia.org/wiki/Conjugate_gradient_method
# solves Ax = b for x, overwrites x
# stops if the error is < ɛ or after reaching max_iter
function truncated_newton!(A::Matrix{T}, b::Vector{T},
                              x::Vector{T}, ɛ::T, max_iter::Int) where T
    r = b - A * x
    p = deepcopy(r)
    Ap = deepcopy(r)
    rsold = dot(r, r)

    n = length(r)

    for i in 1:max_iter
        # Ap[:] = A * p
        mul!(Ap, A, p)
        α = rsold / dot(p, Ap)
        # x += α * p
        BLAS.axpy!(α, p, x)
        # r -= α * Ap
        BLAS.axpy!(-α, Ap, r)
        rsnew = dot(r, r)
        rsnew < ɛ &&  break
        β = rsnew / rsold
        # p[:] = r + β * p
        BLAS.scal!(n, β, p, 1)
        BLAS.axpy!(1, r, p)
        rsold = rsnew
    end
    return x
end

function unmerge_values(data::Vector{T}, assignments::Vector{Int}, n_clusters::Int) where {T <: AbstractFloat}
    ŷ = similar(data)
    elems = map(i -> sum(assignments .== i), 1:n_clusters)
    Idxs = ones(Int64, n_clusters)
    @views @inbounds for i in 1:n_clusters - 1
      # Loop on previous
        Idxs[i + 1] = sum(elems[1:i]) + 1
    end
    @views @inbounds for i in 1:length(data)
        current_cluster = assignments[i]
        search_index = Idxs[current_cluster]

        ŷ[i] = data[search_index]
      # Update corresponding cluster
        Idxs[current_cluster] += 1
    end
    return ŷ
end


text:
@testset "CARLIN algorithm: logistic model" begin
    r = -0.5
    K = 0.8
    F1 = hcat(r)
    F2 = hcat(-r/K)
    X0 = 0.47 .. 0.53
    prob = @ivp(CanonicalQuadraticForm(F1, F2), x(0) ∈ X0)

    sol = solve(prob, T=10.0, alg=CARLIN(N=3, bloat=true))
    @test dim(sol) == 1
    @test ρ([1.0], sol(10.0)) < 0.2

    sol = solve(prob, T=10.0, alg=CARLIN(N=3, bloat=false))
    @test dim(sol) == 1
    @test ρ([1.0], sol(10.0)) < 0.01

    # Use the compressed Kronecker form
    sol = solve(prob, T=10.0, alg=CARLIN(N=3, bloat=false, compress=true))
    @test dim(sol) == 1
    @test ρ([1.0], sol(10.0)) < 0.01
end


text:
using PlotlyJS

function random_line()
    n = 400
    rw() = cumsum(randn(n))
    trace1 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#1f77b4", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="#1f77b4", width=1))
    trace2 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#9467bd", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="rgb(44, 160, 44)", width=1))
    trace3 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#bcbd22", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="#bcbd22", width=1))
    layout = Layout(autosize=false, width=500, height=500,
                    margin=attr(l=0, r=0, b=0, t=65))
    plot([trace1, trace2, trace3], layout)
end

function topo_surface()
    z = Vector[[27.80985, 49.61936, 83.08067, 116.6632, 130.414, 150.7206, 220.1871,
                156.1536, 148.6416, 203.7845, 206.0386, 107.1618, 68.36975, 45.3359,
                49.96142, 21.89279, 17.02552, 11.74317,   14.75226, 13.6671, 5.677561,
                3.31234, 1.156517, -0.147662],
               [27.71966, 48.55022, 65.21374, 95.27666, 116.9964, 133.9056, 152.3412,
                151.934, 160.1139, 179.5327, 147.6184, 170.3943, 121.8194, 52.58537,
                33.08871, 38.40972, 44.24843, 69.5786, 4.019351, 3.050024, 3.039719,
                2.996142, 2.967954, 1.999594],
               [30.4267, 33.47752, 44.80953, 62.47495, 77.43523, 104.2153, 102.7393, 137.0004,
                186.0706, 219.3173, 181.7615, 120.9154, 143.1835, 82.40501, 48.47132,
                74.71461, 60.0909, 7.073525, 6.089851, 6.53745, 6.666096, 7.306965, 5.73684,
                3.625628],
               [16.66549, 30.1086, 39.96952, 44.12225, 59.57512, 77.56929, 106.8925,
                166.5539, 175.2381, 185.2815, 154.5056, 83.0433, 62.61732, 62.33167,
                60.55916, 55.92124, 15.17284, 8.248324, 36.68087, 61.93413, 20.26867,
                68.58819, 46.49812, 0.2360095],
               [8.815617, 18.3516, 8.658275, 27.5859, 48.62691, 60.18013, 91.3286,
                145.7109, 116.0653, 106.2662, 68.69447, 53.10596, 37.92797, 47.95942,
                47.42691, 69.20731, 44.95468, 29.17197, 17.91674, 16.25515, 14.65559,
                17.26048, 31.22245, 46.71704],
               [6.628881, 10.41339, 24.81939, 26.08952, 30.1605, 52.30802, 64.71007,
                76.30823, 84.63686, 99.4324, 62.52132, 46.81647, 55.76606, 82.4099,
                140.2647, 81.26501, 56.45756, 30.42164, 17.28782, 8.302431, 2.981626,
                2.698536, 5.886086, 5.268358],
               [21.83975, 6.63927, 18.97085, 32.89204, 43.15014, 62.86014, 104.6657,
                130.2294, 114.8494, 106.9873, 61.89647, 55.55682, 86.80986, 89.27802,
                122.4221, 123.9698, 109.0952, 98.41956, 77.61374, 32.49031, 14.67344,
                7.370775, 0.03711011, 0.6423392],
               [53.34303, 26.79797, 6.63927, 10.88787, 17.2044, 56.18116, 79.70141,
                90.8453, 98.27675, 80.87243, 74.7931, 75.54661, 73.4373, 74.11694, 68.1749,
                46.24076, 39.93857, 31.21653, 36.88335, 40.02525, 117.4297, 12.70328,
                1.729771, 0],
               [25.66785, 63.05717, 22.1414, 17.074, 41.74483, 60.27227, 81.42432, 114.444,
                102.3234, 101.7878, 111.031, 119.2309, 114.0777, 110.5296, 59.19355,
                42.47175, 14.63598, 6.944074, 6.944075, 27.74936, 0, 0, 0.09449376, 0.07732264],
               [12.827, 69.20554, 46.76293, 13.96517, 33.88744, 61.82613, 84.74799,
                121.122, 145.2741, 153.1797, 204.786, 227.9242, 236.3038, 228.3655,
                79.34425, 25.93483, 6.944074, 6.944074, 6.944075, 7.553681, 0, 0, 0, 0],
               [0, 68.66396, 59.0435, 33.35762, 47.45282, 57.8355, 78.91689, 107.8275,
                168.0053, 130.9597, 212.5541, 165.8122, 210.2429, 181.1713, 189.7617,
                137.3378, 84.65395, 8.677168, 6.956576, 8.468093, 0, 0, 0, 0],
               [0, 95.17499, 80.03818, 59.89862, 39.58476, 50.28058, 63.81641, 80.61302,
                66.37824, 198.7651, 244.3467, 294.2474, 264.3517, 176.4082, 60.21857,
                77.41475, 53.16981, 56.16393, 6.949235, 7.531059, 3.780177, 0, 0, 0],
               [0, 134.9879, 130.3696, 96.86325, 75.70494, 58.86466, 57.20374, 55.18837,
                78.128, 108.5582, 154.3774, 319.1686, 372.8826, 275.4655, 130.2632, 54.93822,
                25.49719, 8.047439, 8.084393, 5.115252, 5.678269, 0, 0, 0],
               [0, 48.08919, 142.5558, 140.3777, 154.7261, 87.9361, 58.11092, 52.83869,
                67.14822, 83.66798, 118.9242, 150.0681, 272.9709, 341.1366, 238.664, 190.2,
                116.8943, 91.48672, 14.0157, 42.29277, 5.115252, 0, 0, 0],
               [0, 54.1941, 146.3839, 99.48143, 96.19411, 102.9473, 76.14089, 57.7844,
                47.0402, 64.36799, 84.23767, 162.7181, 121.3275, 213.1646, 328.482,
                285.4489, 283.8319, 212.815, 164.549, 92.29631, 7.244015, 1.167, 0, 0],
               [0, 6.919659, 195.1709, 132.5253, 135.2341, 89.85069, 89.45549, 60.29967,
                50.33806, 39.17583, 59.06854, 74.52159, 84.93402, 187.1219, 123.9673,
                103.7027, 128.986, 165.1283, 249.7054, 95.39966, 10.00284, 2.39255, 0, 0],
               [0, 21.73871, 123.1339, 176.7414, 158.2698, 137.235, 105.3089, 86.63255, 53.11591,
                29.03865, 30.40539, 39.04902, 49.23405, 63.27853, 111.4215, 101.1956,
                40.00962, 59.84565, 74.51253, 17.06316, 2.435141, 2.287471, -0.0003636982, 0],
               [0, 0, 62.04672, 136.3122, 201.7952, 168.1343, 95.2046, 58.90624, 46.94091,
                49.27053, 37.10416, 17.97011, 30.93697, 33.39257, 44.03077, 55.64542,
                78.22423, 14.42782, 9.954997, 7.768213, 13.0254, 21.73166, 2.156372,
                0.5317867],
               [0, 0, 79.62993, 139.6978, 173.167, 192.8718, 196.3499, 144.6611, 106.5424,
                57.16653, 41.16107, 32.12764, 13.8566, 10.91772, 12.07177, 22.38254,
                24.72105, 6.803666, 4.200841, 16.46857, 15.70744, 33.96221, 7.575688,
                -0.04880907],
               [0, 0, 33.2664, 57.53643, 167.2241, 196.4833, 194.7966, 182.1884, 119.6961,
                73.02113, 48.36549, 33.74652, 26.2379, 16.3578, 6.811293, 6.63927, 6.639271,
                8.468093, 6.194273, 3.591233, 3.81486, 8.600739, 5.21889, 0],
               [0, 0, 29.77937, 54.97282, 144.7995, 207.4904, 165.3432, 171.4047, 174.9216,
                100.2733, 61.46441, 50.19171, 26.08209, 17.18218, 8.468093, 6.63927,
                6.334467, 6.334467, 5.666687, 4.272203, 0, 0, 0, 0],
               [0, 0, 31.409, 132.7418, 185.5796, 121.8299, 185.3841, 160.6566, 116.1478,
                118.1078, 141.7946, 65.56351, 48.84066, 23.13864, 18.12932, 10.28531,
                6.029663, 6.044627, 5.694764, 3.739085, 3.896037, 0, 0, 0],
               [0, 0, 19.58994, 42.30355, 96.26777, 187.1207, 179.6626, 221.3898, 154.2617,
                142.1604, 148.5737, 67.17937, 40.69044, 39.74512, 26.10166, 14.48469,
                8.65873, 3.896037, 3.571392, 3.896037, 3.896037, 3.896037, 1.077756, 0],
               [0.001229679, 3.008948, 5.909858, 33.50574, 104.3341, 152.2165, 198.1988,
                191.841, 228.7349, 168.1041, 144.2759, 110.7436, 57.65214, 42.63504,
                27.91891, 15.41052, 8.056102, 3.90283, 3.879774, 3.936718, 3.968634,
                0.1236256, 3.985531, -0.1835741],
               [0, 5.626141, 7.676256, 63.16226, 45.99762, 79.56688, 227.311, 203.9287,
                172.5618, 177.1462, 140.4554, 123.9905, 110.346, 65.12319, 34.31887,
                24.5278, 9.561069, 3.334991, 5.590495, 5.487353, 5.909499, 5.868994,
                5.833817, 3.568177]]
    trace = surface(z=z)
    layout = Layout(title="Mt. Bruno Elevation", autosize=false, width=500,
                    height=500, margin=attr(l=65, r=50, b=65, t=90))
    plot(trace, layout)
end

function multiple_surface()
    z1 = Vector[[8.83, 8.89, 8.81, 8.87, 8.9, 8.87],
                [8.89, 8.94, 8.85, 8.94, 8.96, 8.92],
                [8.84, 8.9, 8.82, 8.92, 8.93, 8.91],
                [8.79, 8.85, 8.79, 8.9, 8.94, 8.92],
                [8.79, 8.88, 8.81, 8.9, 8.95, 8.92],
                [8.8, 8.82, 8.78, 8.91, 8.94, 8.92],
                [8.75, 8.78, 8.77, 8.91, 8.95, 8.92],
                [8.8, 8.8, 8.77, 8.91, 8.95, 8.94],
                [8.74, 8.81, 8.76, 8.93, 8.98, 8.99],
                [8.89, 8.99, 8.92, 9.1, 9.13, 9.11],
                [8.97, 8.97, 8.91, 9.09, 9.11, 9.11],
                [9.04, 9.08, 9.05, 9.25, 9.28, 9.27],
                [9, 9.01, 9, 9.2, 9.23, 9.2],
                [8.99, 8.99, 8.98, 9.18, 9.2, 9.19],
                [8.93, 8.97, 8.97, 9.18, 9.2, 9.18]]
    z2 = map(x->x+1, z1)
    z3 = map(x->x-1, z1)
    trace1 = surface(z=z1, colorscale="Viridis")
    trace2 = surface(z=z2, showscale=false, opacity=0.9, colorscale="Viridis")
    trace3 = surface(z=z3, showscale=false, opacity=0.9, colorscale="Viridis")
    plot([trace1, trace2, trace3])
end

function clustering_alpha_shapes()
    @eval using DataFrames, RDatasets, Colors

    # load data
    iris = dataset("datasets", "iris")
    nms = unique(iris[:Species])
    colors = [RGB(0.89, 0.1, 0.1), RGB(0.21, 0.50, 0.72), RGB(0.28, 0.68, 0.3)]

    data = GenericTrace[]

    for (i, nm) in enumerate(nms)
        df = iris[iris[:Species] .== nm, :]
        x=df[:SepalLength]
        y=df[:SepalWidth]
        z=df[:PetalLength]
        trace = scatter3d(;name=nm, mode="markers",
                           marker_size=3, marker_color=colors[i], marker_line_width=0,
                           x=x, y=y, z=z)
        push!(data, trace)

        cluster = mesh3d(;color=colors[i], opacity=0.3, x=x, y=y, z=z)
        push!(data, cluster)
    end

    # notice the nested attrs to create complex JSON objects
    layout = Layout(width=800, height=550, autosize=false, title="Iris dataset",
                    scene=attr(xaxis=attr(gridcolor="rgb(255, 255, 255)",
                                          zerolinecolor="rgb(255, 255, 255)",
                                          showbackground=true,
                                          backgroundcolor="rgb(230, 230,230)"),
                               yaxis=attr(gridcolor="rgb(255, 255, 255)",
                                           zerolinecolor="rgb(255, 255, 255)",
                                           showbackground=true,
                                           backgroundcolor="rgb(230, 230,230)"),
                               zaxis=attr(gridcolor="rgb(255, 255, 255)",
                                           zerolinecolor="rgb(255, 255, 255)",
                                           showbackground=true,
                                           backgroundcolor="rgb(230, 230,230)"),
                               aspectratio=attr(x=1, y=1, z=0.7),
                               aspectmode = "manual"))
    plot(data, layout)
end

function scatter_3d()
    @eval using Distributions
    Σ = fill(0.5, 3, 3) + Diagonal([0.5, 0.5, 0.5])
    obs1 = rand(MvNormal(zeros(3), Σ), 200)'
    obs2 = rand(MvNormal(zeros(3), 0.5Σ), 100)'

    trace1 = scatter3d(;x=obs1[:, 1], y=obs1[:, 2], z=obs1[:, 3],
                        mode="markers", opacity=0.8,
                        marker_size=12, marker_line_width=0.5,
                        marker_line_color="rgba(217, 217, 217, 0.14)")

    trace2 = scatter3d(;x=obs2[:, 1], y=obs2[:, 2], z=obs2[:, 3],
                        mode="markers", opacity=0.9,
                        marker=attr(color="rgb(127, 127, 127)",
                                    symbol="circle", line_width=1.0,
                                    line_color="rgb(204, 204, 204)"))

    layout = Layout(margin=attr(l=0, r=0, t=0, b=0))

    plot([trace1, trace2], layout)
end

function trisurf()
    facecolor = repeat([
    	"rgb(50, 200, 200)",
    	"rgb(100, 200, 255)",
    	"rgb(150, 200, 115)",
    	"rgb(200, 200, 50)",
    	"rgb(230, 200, 10)",
    	"rgb(255, 140, 0)"
    ], inner=[2])

    t = mesh3d(
        x=[0, 0, 1, 1, 0, 0, 1, 1],
        y=[0, 1, 1, 0, 0, 1, 1, 0],
        z=[0, 0, 0, 0, 1, 1, 1, 1],
        i=[7, 0, 0, 0, 4, 4, 2, 6, 4, 0, 3, 7],
        j=[3, 4, 1, 2, 5, 6, 5, 5, 0, 1, 2, 2],
        k=[0, 7, 2, 3, 6, 7, 1, 2, 5, 5, 7, 6],
        facecolor=facecolor)

    plot(t)
end

function meshcube()
    t = mesh3d(
        x=[0, 0, 1, 1, 0, 0, 1, 1],
        y=[0, 1, 1, 0, 0, 1, 1, 0],
        z=[0, 0, 0, 0, 1, 1, 1, 1],
        i=[7, 0, 0, 0, 4, 4, 6, 6, 4, 0, 3, 2],
        j=[3, 4, 1, 2, 5, 6, 5, 2, 0, 1, 6, 3],
        k=[0, 7, 2, 3, 6, 7, 1, 1, 5, 5, 7, 6],
        intensity=range(0, stop=1, length=8),
        colorscale=[
            [0, "rgb(255, 0, 255)"],
            [0.5, "rgb(0, 255, 0)"],
            [1, "rgb(0, 0, 255)"]
        ]
    )
    plot(t)
end


text:
function show(io::IO, g::Graph)
    if nv(g) == 0
        print(io, "empty undirected graph")
    else
        print(io, "{$(nv(g)), $(ne(g))} undirected graph")
    end
end

function Graph(n::Int)
    fadjlist = Vector{Vector{Int}}()
    sizehint!(fadjlist,n)
    for i = 1:n
        # sizehint!(i_s, n/4)
        # sizehint!(o_s, n/4)
        push!(fadjlist, Vector{Int}())
    end
    return Graph(1:n, 0, fadjlist)
end

Graph() = Graph(0)

function Graph{T<:Real}(adjmx::AbstractMatrix{T})
    dima,dimb = size(adjmx)
    isequal(dima,dimb) || error("Adjacency / distance matrices must be square")
    issymmetric(adjmx) || error("Adjacency / distance matrices must be symmetric")

    g = Graph(dima)
    for i in find(triu(adjmx))
        ind = ind2sub((dima,dimb),i)
        add_edge!(g,ind...)
    end
    return g
end

function Graph(g::DiGraph)
    gnv = nv(g)

    edgect = 0
    newfadj = deepcopy(g.fadjlist)
    for i in 1:gnv
        for j in badj(g,i)
            if (_insert_and_dedup!(newfadj[i], j))
                edgect += 2     # this is a new edge only in badjlist
            else
                edgect += 1     # this is an existing edge - we already have it
            end
        end
    end

    return Graph(vertices(g), edgect ÷ 2, newfadj)
end

"""Returns the backwards adjacency list of a graph.
For each vertex the Array of `dst` for each edge eminating from that vertex.

NOTE: returns a reference, not a copy. Do not modify result.
"""
badj(g::Graph) = fadj(g)
badj(g::Graph, v::Int) = fadj(g, v)

"""Returns the adjacency list of a graph.
For each vertex the Array of `dst` for each edge eminating from that vertex.

NOTE: returns a reference, not a copy. Do not modify result.
"""
adj(g::Graph) = fadj(g)
adj(g::Graph, v::Int) = fadj(g, v)

function copy(g::Graph)
    return Graph(g.vertices, g.ne, deepcopy(g.fadjlist))
end

==(g::Graph, h::Graph) =
    vertices(g) == vertices(h) &&
    ne(g) == ne(h) &&
    fadj(g) == fadj(h)

"Returns `true` if `g` is a `DiGraph`."
is_directed(g::Graph) = false

function has_edge(g::Graph, e::Edge)
    u, v = e
    u > nv(g) || v > nv(g) && return false
    if degree(g,u) > degree(g,v)
        u, v = v, u
    end
    return length(searchsorted(fadj(g,u), v)) > 0
end

function add_edge!(g::Graph, e::Edge)

    s, d = e
    (s in vertices(g) && d in vertices(g)) || return false
    inserted = _insert_and_dedup!(g.fadjlist[s], d)
    if inserted
        g.ne += 1
    end
    if s != d
        inserted = _insert_and_dedup!(g.fadjlist[d], s)
    end
    return inserted
end

function rem_edge!(g::Graph, e::Edge)
    i = searchsorted(g.fadjlist[src(e)], dst(e))
    length(i) > 0 || return false   # edge not in graph
    i = i[1]
    deleteat!(g.fadjlist[src(e)], i)
    if src(e) != dst(e)     # not a self loop
        i = searchsorted(g.fadjlist[dst(e)], src(e))[1]
        deleteat!(g.fadjlist[dst(e)], i)
    end
    g.ne -= 1
    return true # edge successfully removed
end

"""Add a new vertex to the graph `g`."""
function add_vertex!(g::Graph)
    g.vertices = 1:nv(g)+1
    push!(g.fadjlist, Vector{Int}())

    return true
end

"""Return the number of edges (both ingoing and outgoing) from the vertex `v`."""
degree(g::Graph, v::Int) = indegree(g,v)

doc"""Density is defined as the ratio of the number of actual edges to the
number of possible edges. This is $|v| |v-1|$ for directed graphs and
$(|v| |v-1|) / 2$ for undirected graphs.
"""
density(g::Graph) = (2*ne(g)) / (nv(g) * (nv(g)-1))


text:
module ChaoticEncryption

include("lorenz_key.jl")
include("logistic_key.jl")
include("substitution.jl")

export lorenz_key, logistic_key
export substitution_encryption, substitution_decryption

end


text:
using TFRecord
using Test

@testset "TFRecord.jl" begin

    n = 10
    f1 = rand(Bool, n)
    f2 = rand(1:5, n)
    f3 = rand(("cat", "dog", "chicken", "horse", "goat"), n)
    f4 = rand(Float32, n)

    TFRecord.write(
        "example.tfrecord",
        (
            Dict(
                "feature1" => f1[i],
                "feature2" => f2[i],
                "feature3" => f3[i],
                "feature4" => f4[i],
            )
            for i in 1:n
        )
    )

    reader = TFRecord.read("example.tfrecord")

    for (i, example) in enumerate(reader)
        @test example.features.feature["feature1"].int64_list.value[] == Int(f1[i])
        @test example.features.feature["feature2"].int64_list.value[] == f2[i]
        @test String(example.features.feature["feature3"].bytes_list.value[]) == f3[i]
        @test example.features.feature["feature4"].float_list.value[] == f4[i]
    end
end


text:
function field_cache(v::Number,x)
  nx = length(x)
  c = zeros(typeof(v),nx)
  CachedArray(c)
end

function evaluate_field!(c,v::Number,x)
  nx = length(x)
  setsize!(c,(nx,))
  r = c.array
  for i in eachindex(x)
    @inbounds r[i] = v
  end
  r
end

function field_gradient(v::Number)
  T = typeof(v)
  E = eltype(T)
  zero(E)
end

function field_cache(v::AbstractArray{<:Number},x)
  nx = length(x)
  sv = size(v)
  s = (nx,sv...)
  c = zeros(eltype(v),s)
  CachedArray(c)
end

function evaluate_field!(c,v::AbstractArray{<:Number},x)
  nx = length(x)
  sv = size(v)
  s = (nx,sv...)
  setsize!(c,s)
  cis = CartesianIndices(v)
  r = c.array
  for i in eachindex(x)
    for ci in cis
      @inbounds r[i,ci] = v[ci]
    end
  end
  r
end

function field_gradient(v::AbstractArray{<:Number})
  T = eltype(v)
  E = eltype(T)
  Fill(zero(E),size(v))
end


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11

function predict(Theta1, Theta2, X)
  # Useful values
  m = size(X, 1)
  num_labels = size(Theta2, 1)

  # You need to return the following variables correctly
   p = PREDICT(Theta1, Theta2, X)
  h1 = sigmoid([ones(m, 1) X] * Theta1')
  h2 = sigmoid([ones(m, 1) h1] * Theta2')

  for i in 1:m
    p[i] = findmax(h2[i, :])[2]
  end
  return p
end


text:
type AsymFGMPC <: ParamPC
    θ::Float64

    function AsymFGMPC(θ::Float64)
        if (θ <= 0) | (θ >= 1)
            error("θ parameter must be between 0 and 1")
        end
        new(θ)
    end

end

function params(cop::AsymFGMPC)
    return (cop.θ,)
end

######################
## copula functions ##
######################

function cdf(cop::AsymFGMPC, u1::FloatVec, u2::FloatVec)
    θ = params(cop)[1]
    val = u1.*u2 + θ*u1.*u2.*(1 - u1 - u2 + u1.*u2).*u2.*(1 - u1)
    return val
end

function pdf(cop::AsymFGMPC, u1::FloatVec, u2::FloatVec)
    θ = params(cop)[1]
    val = u1.*(θ*u1.*(6 - 9u2).*u2 + θ*(12u2 - 8).*u2) +
    	θ*(2 - 3u2).*u2 + 1
    return val
end

function hfun(cop::AsymFGMPC, u1::FloatVec, u2::FloatVec)
    θ = params(cop)[1]
    val = -3θ*u1.^3.*u2.^2 + 2θ*u1.^3.*u2 +
         6θ*u1.^2.*u2.^2 - 4θ*u1.^2.*u2 -
         3θ*u1.*u2.^2 + 2θ*u1.*u2 + u1
    return val
end

function vfun(cop::AsymFGMPC, u2::FloatVec, u1::FloatVec)
    θ = params(cop)[1]

    val = -3θ*u1.^2.*u2.^3 + 3θ*u1.^2.*u2.^2 + 4θ*u1.*u2.^3 -
         4θ*u1.*u2.^2 - θ*u2.^3 + θ*u2.^2 + u2
    return val
end


text:
# Without print statement
5
x = 5
x + 1

# With print statement
println(5)
x = 5
println(x + 1)

max_stars_repo_path: Chapter 2/Exercise 2-1.jl

max_stars_repo_name: utilForever/ThinkJulia

max_issues_repo_path: Chapter 2/Exercise 2-1.jl

max_issues_repo_name: utilForever/ThinkJulia

max_forks_repo_path: Chapter 2/Exercise 2-1.jl

max_forks_repo_name: utilForever/ThinkJulia


text:
function getAreaKernel{T <: FloatingPoint}(:: Type{T})
        nType = T == Float64 ? "double" : "float"
        nPi = T == Float64 ? "M_PI" : "M_PI_F"
        nPi4 = T == Float64 ? "M_PI_4" : "M_PI_4_F"

        return "
        #if defined(cl_khr_fp64)  // Khronos extension available?
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
        #elif defined(cl_amd_fp64)  // AMD extension available?
        #pragma OPENCL EXTENSION cl_amd_fp64 : enable
        #endif
        #define number $nType
        #define number4 $(nType)4
        #define PI $nPi
        #define PI_4 $nPi4

        __kernel void area(
                      __global const number *d,
                      __global number4 *dout,
                      const number LONG,
                      const number SHORT) {

            int i = get_global_id(0);

            const number dir = d[i];

            const number area1 = LONG*SHORT/2 * (   PI_4/2-dir - atan((SHORT-LONG)*sin(2*(  PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(  PI_4/2-dir)))));
            const number area2 = LONG*SHORT/2 * ( 3*PI_4/2-dir - atan((SHORT-LONG)*sin(2*(3*PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(3*PI_4/2-dir)))));
            const number area3 = LONG*SHORT/2 * ( 5*PI_4/2-dir - atan((SHORT-LONG)*sin(2*(5*PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(5*PI_4/2-dir)))));
            const number area4 = LONG*SHORT/2 * ( 7*PI_4/2-dir - atan((SHORT-LONG)*sin(2*(7*PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(7*PI_4/2-dir)))));
            const number area5 = LONG*SHORT/2 * ( 9*PI_4/2-dir - atan((SHORT-LONG)*sin(2*(9*PI_4/2-dir)) / (SHORT+LONG + (SHORT-LONG)*cos(2*(9*PI_4/2-dir)))));

            const number lps = SHORT * LONG * PI;

            number4 a;

            a.s0 = (area2-area1)/lps;
            a.s1 = (area3-area2)/lps;
            a.s2 = (area4-area3)/lps;
            a.s3 = (area5-area4)/lps;

            dout[i] = a;
        }
"
end

function areaCL!{T <: FloatingPoint}(
    d_buff,
    area_buff,
    long :: Real,
    d1 :: Int, d2 :: Int,
    ctx, queue, program, :: Type{T})

    k_a = Kernel(program, "area")

    short = one(T)
    call(queue, k_a, d1*d2, nothing, d_buff, area_buff, convert(T, long), short)
end


text:
"""
    Drawing3D

Overarching type for drawings displayable on `Axes3D`.
"""
abstract type Drawing3D <: Drawing end

# NOTE missing inf or nan NOT allowed
# NOTE one at a time, syntax for multiple would be confusing
mutable struct Scatter3D{T} <: Drawing3D
    data       ::T  # data container
    linestyle  ::LineStyle     # line style (color, width, ...)
    markerstyle::MarkerStyle   # marker style (color, size, ...)
    label      ::String        # plot labels (to go in the legend)
end
Scatter3D(d) = default(Scatter3D, d)

#
#
# mutable struct Surface <: Drawing3D
#
# end