text:
using Base: product
using SupplyChainOptimization
#using Cbc
using JuMP
using Test

include("UFLlib.jl")

Seattle = Location(47.608013, -122.335167)

include("UnitTests.jl")

function create_model_storage_customer()
    #storage -> customer
    sc = SupplyChain()

    product = Product("p1")
    add_product!(sc, product)

    c = Customer("c1", Seattle)
    add_customer!(sc, c)
    add_demand!(sc, c, product; demand=[100.0])

    storage = Storage("s1", Seattle; fixed_cost=1000.0, opening_cost=10.0, closing_cost=10.0, initial_opened=true)
    add_storage!(sc, storage)
    add_product!(storage, product; initial_inventory=100)

    add_lane!(sc, Lane(storage, c; unit_cost=1.0))

    return sc
end

function create_test_model2()
    #supplier -> storage -> customer
    sc = SupplyChain()

    product = Product("p1")
    add_product!(sc, product)

    c = Customer("c1", Seattle)
    add_customer!(sc, c)
    add_demand!(sc, c, product; demand=[100.0])

    storage = Storage("s1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_storage!(sc, storage)
    add_product!(storage, product)

    supplier = Supplier("supplier1", Seattle)
    add_supplier!(sc, supplier)
    add_product!(supplier, product; unit_cost=0.0, maximum_throughput=Inf)

    add_lane!(sc, Lane(storage, c; unit_cost=1.0))
    add_lane!(sc, Lane(supplier, storage; unit_cost=1.0))

    return sc, product, supplier
end

function create_test_model3()
    #plant -> storage -> customer
    sc = SupplyChain()

    product = Product("p1")
    add_product!(sc, product)

    customer = Customer("c1", Seattle)
    add_customer!(sc, customer)
    add_demand!(sc, customer, product; demand=[100.0])

    storage = Storage("s1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_storage!(sc, storage)
    add_product!(storage, product)

    plant = Plant("plant1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_plant!(sc, plant)
    add_product!(plant, product; bill_of_material=Dict{Product, Float64}(), unit_cost=1, maximum_throughput=Inf)

    add_lane!(sc, Lane(storage, customer; unit_cost=1.0))
    add_lane!(sc, Lane(plant, storage; unit_cost=1.0))

    return sc, product, plant
end

function create_test_model4()
    #supplier -> plant -> storage -> customer
    sc = SupplyChain()

    product1 = add_product!(sc, Product("p1"))
    product2 = add_product!(sc, Product("p2"))

    supplier = Supplier("supplier1", Seattle)
    add_supplier!(sc, supplier)
    add_product!(supplier, product1; unit_cost=0.0, maximum_throughput=Inf)

    customer = Customer("c1", Seattle)
    add_customer!(sc, customer)
    add_demand!(sc, customer, product2; demand=[100.0])

    storage = Storage("s1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_storage!(sc, storage)
    add_product!(storage, product2)

    plant = Plant("plant1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_plant!(sc, plant)
    add_product!(plant, product2; bill_of_material=Dict{Product, Float64}(product1 => 1), unit_cost=1, maximum_throughput=Inf)

    add_lane!(sc, Lane(storage, customer; unit_cost=1.0))
    add_lane!(sc, Lane(plant, storage; unit_cost=1.0))
    add_lane!(sc, Lane(supplier, plant; unit_cost=1.0))

    return sc, product2, plant
end

function create_test_model5()
    #supplier -> plant -> storage -> customer
    sc = SupplyChain(2)

    product1 = add_product!(sc, Product("p1"))
    product2 = add_product!(sc, Product("p2"))

    supplier = add_supplier!(sc, Supplier("supplier1", Seattle))
    add_product!(supplier, product1; unit_cost=0.0, maximum_throughput=Inf)

    customer = add_customer!(sc, Customer("c1", Seattle))
    add_demand!(sc, customer, product2; demand=[100.0, 100.0])

    storage = add_storage!(sc, Storage("s1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false))
    add_product!(storage, product2)
    plant = add_plant!(sc, Plant("plant1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false))
    add_product!(plant, product2; bill_of_material=Dict{Product, Float64}(product1 => 1), unit_cost=1, maximum_throughput=Inf)
    add_lane!(sc, Lane(storage, customer; unit_cost=1.0))
    add_lane!(sc, Lane(plant, storage; unit_cost=1.0))
    add_lane!(sc, Lane(supplier, plant; unit_cost=1.0))

    return sc, product2, plant
end

function create_test_model6()
    #supplier -> plant -> storage -> customer
    sc = SupplyChain(2)

    product1 = Product("p1")
    product2 = Product("p2")
    add_product!(sc, product1)
    add_product!(sc, product2)

    supplier = Supplier("supplier1", Seattle)
    add_supplier!(sc, supplier)
    add_product!(supplier, product1; unit_cost=0.0, maximum_throughput=Inf)

    plant = Plant("plant1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_plant!(sc, plant)
    add_product!(plant, product2; bill_of_material=Dict{Product, Float64}(product1 => 1), unit_cost=1, maximum_throughput=Inf)
    add_lane!(sc, Lane(supplier, plant; unit_cost=1.0))

    for i in 1:500
        customer = Customer("c$i", Seattle)
        add_customer!(sc, customer)
        add_demand!(sc, customer, product2; demand=[100.0, 100.0])
    end

    for i in 1:50
        storage = Storage("s$i", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
        add_storage!(sc, storage)
        add_product!(storage, product2)
        add_lane!(sc, Lane(plant, storage; unit_cost=1.0))
    end

    for customer in sc.customers, storage in sc.storages
        add_lane!(sc, Lane(storage, customer; unit_cost=1.0))
    end

    return sc, product2, plant
end

function create_test_model7()
    #storage -> customer
    sc = SupplyChain()

    product = Product("p1")
    add_product!(sc, product)

    c = Customer("c1", Seattle)
    add_customer!(sc, c)
    add_demand!(sc, c, product; demand=[100.0], service_level=0.0)

    storage = Storage("s1", Seattle; fixed_cost=1000.0, opening_cost=10.0, closing_cost=10.0, initial_opened=true)
    add_storage!(sc, storage)
    add_product!(storage, product; initial_inventory=100)

    add_lane!(sc, Lane(storage, c; unit_cost=1.0))

    return sc
end

function create_test_broken_model()
    #supplier -> plant -> storage -> customer
    sc = SupplyChain(2)

    product1 = Product("p1")
    product2 = Product("p2")
    add_product!(sc, product1)
    add_product!(sc, product2)

    supplier = Supplier("supplier1", Seattle)
    add_supplier!(sc, supplier)
    add_product!(supplier, product1; unit_cost=0.0, maximum_throughput=Inf)

    customer = Customer("c1", Seattle)
    add_customer!(sc, customer)
    add_demand!(sc, customer, product2; demand=[100.0, 100.0])

    storage = Storage("s1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_storage!(sc, storage)
    add_product!(storage, product2)
    plant = Plant("plant1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_plant!(sc, plant)

    add_lane!(sc, Lane(storage, customer; unit_cost=1.0))
    add_lane!(sc, Lane(plant, storage; unit_cost=1.0))
    add_lane!(sc, Lane(supplier, plant; unit_cost=1.0))

    return sc, product2, plant
end

function create_test_infeasible_model()
    #supplier -> plant -> storage -> customer
    sc = SupplyChain(2)

    product1 = add_product!(sc, Product("p1"))
    product2 = add_product!(sc, Product("p2"))

    supplier = Supplier("supplier1", Seattle)
    add_supplier!(sc, supplier)
    add_product!(supplier, product1; unit_cost=0.0, maximum_throughput=Inf)

    customer = Customer("c1", Seattle)
    add_customer!(sc, customer)
    add_demand!(sc, customer, product2; demand=[100.0, 100.0])

    storage = Storage("s1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_storage!(sc, storage)
    add_product!(storage, product2)

    plant = Plant("plant1", Seattle; fixed_cost=1000.0, opening_cost=500.0, closing_cost=500.0, initial_opened=false)
    add_plant!(sc, plant)

    add_lane!(sc, Lane(storage, customer; unit_cost=1.0))
    add_lane!(sc, Lane(plant, storage; unit_cost=1.0))
    add_lane!(sc, Lane(supplier, plant; unit_cost=1.0))

    return sc, product2, plant
end

@testset "Happy Path" begin

@test !isnothing(create_model_storage_customer())

@test begin
    sc = create_model_storage_customer()
    SupplyChainOptimization.optimize_network!(sc)

    get_total_costs(sc) == 1100 &&
    get_total_fixed_costs(sc) == 1000 &&
    get_total_transportation_costs(sc) == 100
end

@test begin
    sc = create_model_storage_customer()
    SupplyChainOptimization.optimize_network!(sc)

    get_shipments(sc, first(sc.storages), first(sc.products)) == 100 &&
    is_opened(sc, first(sc.storages)) &&
    !is_opening(sc, first(sc.storages)) &&
    !is_closing(sc, first(sc.storages)) &&
    get_inventory_at_start(sc, first(sc.storages), first(sc.products), 1) == 100 &&
    get_inventory_at_end(sc, first(sc.storages), first(sc.products), 1) == 0
end

@test begin
    sc = create_model_storage_customer()
    SupplyChainOptimization.optimize_network!(sc)
    plot_costs(sc)
    true
end

@test begin
    sc = create_model_storage_customer()
    SupplyChainOptimization.optimize_network!(sc)
    plot_flows(sc)
    true
end

@test begin
    sc = create_model_storage_customer()
    SupplyChainOptimization.optimize_network!(sc)
    plot_network(sc)
    true
end

@test begin
    sc, product, supplier = create_test_model2()
    SupplyChainOptimization.optimize_network!(sc)
    get_total_costs(sc) == 1000 + 500 + 200 &&
    get_shipments(sc, supplier, product, 1) == 100
end

@test begin
    sc, product, plant = create_test_model3()
    SupplyChainOptimization.optimize_network!(sc)
    get_total_costs(sc) == 3300 &&
    get_production(sc, plant, product, 1) == 100
end

@test begin
    sc, product2, plant = create_test_model4()
    SupplyChainOptimization.optimize_network!(sc)
    get_total_costs(sc) == 3400 &&
    get_production(sc, plant, product2, 1) == 100
end

@test begin
    sc, product2, plant = create_test_model5()
    SupplyChainOptimization.optimize_network!(sc)
    get_total_costs(sc) == 1500 + 3000 + 600 + 200 && get_production(sc, plant, product2, 1)  == 200
end

@test begin
    sc = create_test_model7()
    SupplyChainOptimization.optimize_network!(sc)
    get_total_costs(sc) == 10
end
end

@testset "Infeasible" begin
    @test  begin
        sc, product2, plant = create_test_infeasible_model()
        SupplyChainOptimization.optimize_network!(sc)
        status = termination_status(sc.optimization_model)
        status == JuMP.MathOptInterface.INFEASIBLE
    end
end

@testset "Invalid" begin
    #@test_throws ArgumentError begin
    #    sc, product2, plant = create_test_broken_model()
    #    SupplyChainOptimization.optimize!(sc)
    #    status = termination_status(sc.optimization_model)
    #end
end

@testset "Scaling" begin
    @test begin
        sc, product2, plant = create_test_model6()
        SupplyChainOptimization.optimize_network!(sc)
        true
    end
end


text:
struct SMC{R} <: AbstractMCMC.AbstractSampler
    nparticles::Int
    resampler::R
end

"""
    SMC(n[, resampler = ResampleWithESSThreshold()])
    SMC(n, [resampler = resample_systematic, ]threshold)

Create a sequential Monte Carlo (SMC) sampler with `n` particles.

If the algorithm for the resampling step is not specified explicitly, systematic resampling
is performed if the estimated effective sample size per particle drops below 0.5.
"""
SMC(nparticles::Int) = SMC(nparticles, ResampleWithESSThreshold())

# Convenient constructors with ESS threshold
function SMC(nparticles::Int, resampler, threshold::Real)
    return SMC(nparticles, ResampleWithESSThreshold(resampler, threshold))
end
SMC(nparticles::Int, threshold::Real) = SMC(nparticles, resample_systematic, threshold)

struct SMCSample{P,W,L}
    trajectories::P
    weights::W
    logevidence::L
end

function AbstractMCMC.sample(model::AbstractMCMC.AbstractModel, sampler::SMC; kwargs...)
    return AbstractMCMC.sample(Random.GLOBAL_RNG, model, sampler; kwargs...)
end

function AbstractMCMC.sample(
    rng::Random.AbstractRNG, model::AbstractMCMC.AbstractModel, sampler::SMC; kwargs...
)
    if !isempty(kwargs)
        @warn "keyword arguments $(keys(kwargs)) are not supported by `SMC`"
    end

    # Create a set of particles.
    particles = ParticleContainer(
        [Trace(model, TracedRNG()) for _ in 1:(sampler.nparticles)], TracedRNG()
    )

    # Perform particle sweep.
    logevidence = sweep!(rng, particles, sampler.resampler)

    return SMCSample(collect(particles), getweights(particles), logevidence)
end

struct PG{R} <: AbstractMCMC.AbstractSampler
    """Number of particles."""
    nparticles::Int
    """Resampling algorithm."""
    resampler::R
end

"""
    PG(n, [resampler = ResampleWithESSThreshold()])
    PG(n, [resampler = resample_systematic, ]threshold)

Create a Particle Gibbs sampler with `n` particles.

If the algorithm for the resampling step is not specified explicitly, systematic resampling
is performed if the estimated effective sample size per particle drops below 0.5.
"""
PG(nparticles::Int) = PG(nparticles, ResampleWithESSThreshold())

# Convenient constructors with ESS threshold
function PG(nparticles::Int, resampler, threshold::Real)
    return PG(nparticles, ResampleWithESSThreshold(resampler, threshold))
end
PG(nparticles::Int, threshold::Real) = PG(nparticles, resample_systematic, threshold)

struct PGState{T}
    trajectory::T
end

struct PGSample{T,L}
    trajectory::T
    logevidence::L
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG, model::AbstractMCMC.AbstractModel, sampler::PG; kwargs...
)
    # Create a new set of particles.
    particles = ParticleContainer(
        [Trace(model, TracedRNG()) for _ in 1:(sampler.nparticles)], TracedRNG()
    )

    # Perform a particle sweep.
    logevidence = sweep!(rng, particles, sampler.resampler)

    # Pick a particle to be retained.
    trajectory = rand(rng, particles)

    return PGSample(trajectory, logevidence), PGState(trajectory)
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::AbstractMCMC.AbstractModel,
    sampler::PG,
    state::PGState;
    kwargs...,
)
    # Create a new set of particles.
    nparticles = sampler.nparticles
    x = map(1:nparticles) do i
        if i == nparticles
            # Create reference trajectory.
            forkr(state.trajectory)
        else
            Trace(model, TracedRNG())
        end
    end
    particles = ParticleContainer(x, TracedRNG())

    # Perform a particle sweep.
    logevidence = sweep!(rng, particles, sampler.resampler, particles.vals[nparticles])

    # Pick a particle to be retained.
    newtrajectory = rand(rng, particles)

    return PGSample(newtrajectory, logevidence), PGState(newtrajectory)
end

ext: jl


text:
using Documenter, ANaiveJuliaTutorial

makedocs(;
    modules=[ANaiveJuliaTutorial],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/singularitti/ANaiveJuliaTutorial.jl/blob/{commit}{path}#L{line}",
    sitename="ANaiveJuliaTutorial.jl",
    authors="Qi Zhang <singularitti@outlook.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/singularitti/ANaiveJuliaTutorial.jl",
)


text:
module CUDA

using GPUCompiler

using GPUArrays

using LLVM
using LLVM.Interop
using Core: LLVMPtr

using Adapt

using Requires

using LinearAlgebra

using BFloat16s

using Memoize

## source code includes

include("pointer.jl")

# core library
include("../lib/utils/APIUtils.jl")
include("../lib/cudadrv/CUDAdrv.jl")

# essential stuff
include("initialization.jl")
include("state.jl")
include("debug.jl")

# binary dependencies
include("../deps/discovery.jl")
include("../deps/compatibility.jl")
include("../deps/bindeps.jl")

# device functionality (needs to be loaded first, because of generated functions)
include("device/pointer.jl")
include("device/array.jl")
include("device/intrinsics.jl")
include("device/llvm.jl")
include("device/runtime.jl")
include("device/texture.jl")

include("pool.jl")

# compiler libraries
include("../lib/cupti/CUPTI.jl")
include("../lib/nvtx/NVTX.jl")
export CUPTI, NVTX

# compiler implementation
include("compiler/gpucompiler.jl")
include("compiler/execution.jl")
include("compiler/exceptions.jl")
include("compiler/reflection.jl")

# array abstraction
include("array.jl")
include("gpuarrays.jl")
include("utilities.jl")
include("texture.jl")

# array libraries
include("../lib/complex.jl")
include("../lib/library_types.jl")
include("../lib/cublas/CUBLAS.jl")
include("../lib/cusparse/CUSPARSE.jl")
include("../lib/cusolver/CUSOLVER.jl")
include("../lib/cufft/CUFFT.jl")
include("../lib/curand/CURAND.jl")
include("../lib/cudnn/CUDNN.jl")
include("../lib/cutensor/CUTENSOR.jl")
export CUBLAS, CUSPARSE, CUSOLVER, CUFFT, CURAND, CUDNN, CUTENSOR

# integrations and specialized functionality
include("indexing.jl")
include("broadcast.jl")
include("mapreduce.jl")
include("accumulate.jl")
include("reverse.jl")
include("linalg.jl")
include("iterator.jl")
include("statistics.jl")
include("random.jl")
include("sorting.jl")

# other libraries
include("../lib/nvml/NVML.jl")
const has_nvml = NVML.has_nvml
export NVML, has_nvml

include("deprecated.jl")

include("precompile.jl")
_precompile_()

end


text:
const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

for sym in names(@__MODULE__, all=true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_") || startswith(sym_string, "@_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") &&
         Base.isidentifier(sym_string[2:end])))
       continue
    end
    @eval export $sym
end

max_stars_repo_name: Fuinn/mos-interface-jl

max_issues_repo_name: Fuinn/mos-interface-jl

max_forks_repo_name: Fuinn/mos-interface-jl


text:
# ---
# title: 1346. Check If N and Its Double Exist
# id: problem1346
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/check-if-n-and-its-double-exist/description/>
# hidden: true
# ---
#
# Given an array `arr` of integers, check if there exists two integers `N` and
# `M` such that `N` is the double of `M` ( i.e. `N = 2 * M`).
#
# More formally check if there exists two indices `i` and `j` such that :
#
#   * `i != j`
#   * `0 <= i, j < arr.length`
#   * `arr[i] == 2 * arr[j]`
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [10,2,5,3]
#     Output: true
#     Explanation: N = 10 is the double of M = 5,that is, 10 = 2 * 5.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [7,1,14,11]
#     Output: true
#     Explanation: N = 14 is the double of M = 7,that is, 14 = 2 * 7.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [3,1,7,11]
#     Output: false
#     Explanation: In this case does not exist N and M, such that N = 2 * M.
#
#
#
#
# **Constraints:**
#
#   * `2 <= arr.length <= 500`
#   * `-10^3 <= arr[i] <= 10^3`
#
#
## @lc code=start
using LeetCode

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


text:
## ExoplanetsSysSim/src/planet.jl
## (c) 2015 Eric B. Ford

# export Planet

struct Planet
  radius::Float64       # solar radii
  mass::Float64         # solar masses
  id::Int64              # id number (for purposes of tracking or grouping planets)
end

function Planet(radius::Float64, mass::Float64; id::Int64=0)
  pl = Planet(radius, mass, id)
end

#Planet() = Planet(0.0, 0.0)  # Commented out, so don't accidentally have invalid Planets

function test_planet_constructors(sim_param::SimParam)
  #blank = Planet()
  earth = Planet(0.0091705248,3e-6)
end


text:
stress_scale = 1.0e+5
strain_scale = 1

# tid = parse(Int64, ARGS[1])
force_scale = 5.0
fiber_size = 2

testtype = "NeuralNetwork2D"
nntype = "piecewise"
include("nnutil.jl")

H0 = [1.0406424793819175e6 209077.08366547766         0.0
      209077.08366547766   1.0411467691352057e6       0.0
      0.0                  0.0                   419057.32049008965]/stress_scale

T = 200.0
NT = 200

nx, ny = 10, 5
porder = 2
nxf, nyf =10*fiber_size,5*fiber_size
Lx, Ly = 1.0, 0.5
prop_dummy = Dict("name"=> "PlaneStress","rho"=> 4.5, "E"=> 1e+6, "nu"=> 0.2)  #dummy

#scale  length uses cm, stress uses     GPa,     time uses  s,
#       current is 10cm    and      10-4GPa              ms
scales = [10.0, 1.0e-4, 1.0e-3]

function BuildDomain(nx::Int64, ny::Int64, Lx::Float64, Ly::Float64, tid::Int64, porder::Int64, force_scale::Float64, prop::Dict{String, Any})
    nodes, EBC, g, gt, FBC, fext, ft = BoundaryCondition(tid, nx, ny,porder, Lx, Ly; force_scale=force_scale)

ndofs=2
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords,elnodes, prop, 3))
    end
end

domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext), ft, gt
end
# write  to a file named: testnn$(idx)_$(nntype)_from$(restart_id)_test$(tid).txt
# the displacement "disp"
# in the file, the first line has
# Time    ux_tr    uy_tr    ux_tb    uy_tb
# Each tid, we save ux uy at right-top corner and right-bottom corner
# and the stress history "stress"

function TestNN(tid::Int64, nx::Int64, ny::Int64, prop::Dict{String, Any}, file_name::String, restart_id::Int64)
    printstyled("tid=$tid\n", color=:green)

domain, ft, gt = BuildDomain(nx, ny, Lx, Ly, tid, porder, force_scale, prop)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

ts = LinRange(0, T, NT+1)
adaptive_solver_args = Dict("Newmark_rho"=> 0.0,
                          "Newton_maxiter"=>10,
                          "Newton_Abs_Err"=>1e-4,
                          "Newton_Rel_Err"=>1e-6,
                          "damped_Newton_eta" => 1.0)

globdat, domain, ts = AdaptiveSolver("NewmarkSolver", globdat, domain, T, NT, adaptive_solver_args)

Disp = zeros(Float64, length(ts), 5)
Disp[:,1] = ts
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], nx*porder+1, ny*porder+1)[end,end] for i = 1:length(domain.history["state"])]
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], nx*porder+1, ny*porder+1)[end,end] for i = 1:length(domain.history["state"])]
Disp[:,2], Disp[:,3] = ux, uy
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], nx*porder+1, ny*porder+1)[div((nx*porder+1),2),end] for i = 1:length(domain.history["state"])]
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], nx*porder+1, ny*porder+1)[div((nx*porder+1),2),end] for i = 1:length(domain.history["state"])]
Disp[:,4], Disp[:,5]  = ux, uy

#file = matopen("Plot/test_nntrain$(idx)_$(nntype)_from$(restart_id)_test$(tid).txt", "w")
file = matopen(file_name, "w")
write(file, "state", Disp)
write(file, "full_stress_history", domain.history["stress"])
write(file, "full_state_history", domain.history["state"])
close(file)
end

function Reference(tid::Int64, nx::Int64, ny::Int64)
    full_state_history, full_fext_history = read_data("$(@__DIR__)/Data/order$porder/$(tid)_$(force_scale)_$(fiber_size).dat")
    strain, stress = read_strain_stress("Data/order$porder/$(tid)_$(force_scale)_$(fiber_size).dat")
    N = length(full_state_history)
    ts = LinRange(0, T, N)
    Disp = zeros(Float64, N, 5)
    Disp[:,1] = ts
    ux = [reshape(full_state_history[i][1:(nx*porder+1)*(ny*porder+1)], nx*porder+1, ny*porder+1)[end,end] for i = 1:N]
    uy = [reshape(full_state_history[i][(nx*porder+1)*(ny*porder+1)+1:end], nx*porder+1, ny*porder+1)[end,end] for i = 1:N]
    Disp[:,2], Disp[:,3] = ux, uy
    ux = [reshape(full_state_history[i][1:(nx*porder+1)*(ny*porder+1)], nx*porder+1, ny*porder+1)[div((nx*porder+1),2),end] for i = 1:N]
    uy = [reshape(full_state_history[i][(nx*porder+1)*(ny*porder+1)+1:end], nx*porder+1, ny*porder+1)[div((nx*porder+1),2),end] for i = 1:N]
    Disp[:,4], Disp[:,5]  = ux, uy

    file = matopen("Plot/reference$(tid).txt", "w")
    write(file, "state", Disp)
    write(file, "full_stress_history", stress)
    write(file, "full_state_history", full_state_history)
    close(file)
end

function Plot(tid::Int64, restart_id::Int64)
    close("all")
    L_scale, t_scale = scales[1], scales[3]

    # Reference
    vars = matread("Plot/reference$(tid).txt")
    disp_ref = vars["state"]
    markevery=5

    plot(disp_ref[1:markevery:end, 1]*t_scale, disp_ref[1:markevery:end, 2]*L_scale, "or", fillstyle="none", label="Reference")
    plot(disp_ref[1:markevery:end, 1]*t_scale, disp_ref[1:markevery:end, 3]*L_scale, "oy", fillstyle="none")
    plot(disp_ref[1:markevery:end, 1]*t_scale, disp_ref[1:markevery:end, 4]*L_scale, "ob", fillstyle="none")
    plot(disp_ref[1:markevery:end, 1]*t_scale, disp_ref[1:markevery:end, 5]*L_scale, "og", fillstyle="none")

    # NNTrain
    vars = matread("Plot/test_nntrain$(idx)_$(nntype)_from$(restart_id)_test$(tid).txt")
    disp_test = vars["state"]
    plot(disp_test[:, 1]*t_scale, disp_test[:, 2]*L_scale, "r",  label="Indirect")
    plot(disp_test[:, 1]*t_scale, disp_test[:, 3]*L_scale, "y")
    plot(disp_test[:, 1]*t_scale, disp_test[:, 4]*L_scale, "b")
    plot(disp_test[:, 1]*t_scale, disp_test[:, 5]*L_scale, "g")

    # NNLearn
    vars = matread("Plot/test_nnlearn$(idx)_$(nntype)_test$(tid).txt")
    #use fint, strain, stress, for debugging purpose
    disp_test = vars["state"]
    plot(disp_test[:, 1]*t_scale, disp_test[:, 2]*L_scale, "--r",  label="Direct")
    plot(disp_test[:, 1]*t_scale, disp_test[:, 3]*L_scale, "--y")
    plot(disp_test[:, 1]*t_scale, disp_test[:, 4]*L_scale, "--b")
    plot(disp_test[:, 1]*t_scale, disp_test[:, 5]*L_scale, "--g")

    xlabel("Time (s)")
    ylabel("Displacement (cm)")
    PyPlot.tight_layout()
    savefig("Plot/plate_plasticity_disp_nn$(idx)_$(nntype)_from$(restart_id)_test$(tid).pdf")
    mpl.save("Plot/plate_plasticity_disp_nn$(idx)_$(nntype)_from$(restart_id)_test$(tid).tex")
end

function PlotStress(tid::Int64, nx::Int64, ny::Int64, nxf::Int64, nyf::Int64, fiber_size::Int64, porder::Int64,
    Lx::Float64, Ly::Float64, force_scale::Float64, restart_id::Int64)

    #visualize exact solution
    close("all")
    domain, _, _ = BuildDomain(nxf, nyf, Lx, Ly, tid, porder, force_scale, prop_dummy)
    vars = matread("Plot/reference$(tid).txt")
    ts, stress, full_state_history = vars["state"][:,1], vars["full_stress_history"], vars["full_state_history"]
    # stress size is NT
    # full_state_history size is NT+1
    frame = Int64(NT/2)
    vmin, vmax = visσ(domain, nxf, nyf,  stress[frame], full_state_history[frame+1]; scaling = scales)

    xlabel("X (cm)")
    ylabel("Y (cm)")
    PyPlot.tight_layout()

    savefig("Plot/plate_plasticity_stress_reference$(tid).pdf")

    #visualize nntrain solution
    close("all")
    domain,_,_ = BuildDomain(nx, ny, Lx, Ly, tid, porder, force_scale, prop_dummy)
    vars = matread("Plot/test_nntrain$(idx)_$(nntype)_from$(restart_id)_test$(tid).txt")
    ts, stress, full_state_history = vars["state"][:,1], vars["full_stress_history"], vars["full_state_history"]
    # stress size is NT
    # full_state_history size is NT+1
    frame = 0
    for iframe = 1 : length(ts)
        if ts[iframe] > T/2.0
            frame = ((ts[iframe] - T/2.0) > (T/2.0 - ts[iframe - 1]) ? iframe-1  : iframe)
            break
        end
    end
    @show "frame is ", frame, " ,time is ", ts[frame]

    visσ(domain, nx, ny,  stress[frame-1], full_state_history[frame], vmin, vmax; scaling = scales)
    xlabel("X (cm)")
    ylabel("Y (cm)")
    PyPlot.tight_layout()
    savefig("Plot/plate_plasticity_stress_test_nntrain$(idx)_$(nntype)_from$(restart_id)_test$(tid).pdf")

    #visualize nnlearn solution
    close("all")
    domain,_,_ = BuildDomain(nx, ny, Lx, Ly, tid, porder, force_scale, prop_dummy)
    vars = matread("Plot/test_nnlearn$(idx)_$(nntype)_test$(tid).txt")
    ts, stress, full_state_history = vars["state"][:,1], vars["full_stress_history"], vars["full_state_history"]
    # stress size is NT
    # full_state_history size is NT+1
    frame = 0
    for iframe = 1 : length(ts)
        if ts[iframe] > T/2.0
            frame = ((ts[iframe] - T/2.0) > (T/2.0 - ts[iframe - 1]) ? iframe-1  : iframe)
            break
        end
    end
    @show "frame is ", frame, " ,time is ", ts[frame]

    visσ(domain, nx, ny,  stress[frame-1], full_state_history[frame], vmin, vmax; scaling = scales)
    xlabel("X (cm)")
    ylabel("Y (cm)")
    PyPlot.tight_layout()
    savefig("Plot/plate_plasticity_stress_test_nnlearn$(idx)_$(nntype)_test$(tid).pdf")

end

#####################################################
GENERATE_DATA = false
PLOT = true
restart_id_list = [3]
tid_list = [106, 206, 300]
train_id = 50

if GENERATE_DATA
for restart_id in restart_id_list
    # Reference
    for tid in tid_list
        Reference(tid, nxf, nyf)
    end

    # NNTrain
    s = ae_to_code("Data/$(nntype)/NN_Train_$(idx)_from_$(restart_id)_$(H_function)_ite$(train_id).mat", nntype)
    eval(Meta.parse(s))
    prop = Dict("name"=> testtype, "rho"=> 4.5, "nn"=>post_nn)
    for tid in tid_list
        file_name = "Plot/test_nntrain$(idx)_$(nntype)_from$(restart_id)_test$(tid).txt"
        TestNN(tid, nx, ny, prop, file_name, restart_id)
    end

    # NNLearn
    s = ae_to_code("Data/$(nntype)/NNLearn_$(idx)_$(H_function)_ite$(train_id).mat", nntype)
    eval(Meta.parse(s))
    prop = Dict("name"=> testtype, "rho"=> 4.5, "nn"=>post_nn)
    for tid in tid_list
        file_name = "Plot/test_nnlearn$(idx)_$(nntype)_test$(tid).txt"
        TestNN(tid, nx, ny, prop, file_name, restart_id)
    end

end
end

if PLOT
for tid in tid_list
    for restart_id in restart_id_list
        Plot(tid, restart_id)
        PlotStress(tid, nx, ny, nxf, nyf, fiber_size, porder, Lx, Ly, force_scale, restart_id)
    end
end
end


text:
const C = CartesianIndex{2}

function day25()
    part = [0, 0]

    lines = filter(s -> !isempty(strip(s)), strip.(readlines("AoCdata/AoC_2021_day25.txt")))
    nrows, ncols = length(lines), length(first(lines))
    mat = fill('.', nrows, ncols)
    for i in 1:nrows
        mat[i, :] .= collect(lines[i])
    end
    tmat, moved, steps = deepcopy(mat), 1, 0
    while moved > 0
        moved = 0
        for c in CartesianIndices(mat)
            if mat[c] == '>'
                dest = C(c[1], mod1(c[2] + 1, ncols))
                if mat[dest] == '.'
                    tmat[c] = '.'
                    tmat[dest] = '>'
                    moved += 1
                end
            end
        end
        mat = deepcopy(tmat)
        for c in CartesianIndices(mat)
            if mat[c] == 'v'
                dest = C(mod1(c[1] + 1, nrows), c[2])
                if mat[dest] == '.'
                    tmat[c] = '.'
                    tmat[dest] = 'v'
                    moved += 1
                end
            end
        end
        mat = deepcopy(tmat)
        steps += 1
    end

    part[1] = steps

    return part
end

part = day25()
println("Part 1: ", part[1])


text:
using SpikingNetworks;@unit

parseParameters"""
N=4000
area = 20000*um^2
Cm = 1*ufarad*cm^(-2)*area
gl = 5e-5*siemens*cm^(-2)*area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm^(-2)*area
g_kd = 30*msiemens*cm^(-2)*area
VT = -63*mV
taue = 5*ms
taui = 10*ms
Ee = 0*mV
Ei = -80*mV"""

parseParameters"""
N=4000;M=4*N÷5
we=zeros(N,N)
wi=zeros(N,N)
we[1:N,1:M]=(rand(N,M).<0.02)*6*nS
wi[1:N,M+1:N]=(rand(N,N-M).<0.02)*67*nS"""

parseEquations"""
dv/dt=(gl*(El-v)+ge*(Ee-v)+gi*(Ei-v)-g_na*m^3*h*(v-ENa)-g_kd*n^4*(v-EK))/Cm
dm/dt=0.32*(13-v+VT)/(exp((13-v+VT)/4)-1)*(1-m)-0.28*(v-VT-40)/(exp((v-VT-40)/5)-1)*m
dn/dt=0.032*(15-v+VT)/(exp((15-v+VT)/5)-1)*(1-n)-0.5*exp((10-v+VT)/40)*n
dh/dt=0.128*exp((17-v+VT)/18)*(1-h)-4/(1+exp((40-v+VT)/5))*h
dge/dt=-ge/taue
dgi/dt=-gi/taui"""

net=init"""
N=N
v=El+5*(randn(N)-1)
m=0*ones(N)
n=0*ones(N)
h=1*ones(N)
ge=(randn(N)*1.5+4)*10*nS
gi=(randn(N)*12+20)*10*nS
we=we
wi=wi"""

spike"""
v.>-20
3"""

synapse"""
ge+=we
gi+=wi"""

stateMonitor"[1,10,100]"

@time run!(net,100/net.dt)
using MatlabPlot
mplot(net.vRecord')
# workspace()


text:
## Panel helpers
#============================================#

# Evaluation of velocities on points
source_velocity(σ1, σ2, panel :: AbstractPanel2D, x, y) = panel_velocity(linear_source_velocity_a, linear_source_velocity_b, σ1, σ2, panel, x, y)
vortex_velocity(γ1, γ2, panel :: AbstractPanel2D, x, y) = panel_velocity(linear_vortex_velocity_a, linear_vortex_velocity_b, γ1, γ2, panel, x, y)

total_velocity(σ_1j :: Real, σ_2j :: Real, γ_1j :: Real, γ_2j :: Real, panel_j :: AbstractPanel2D, x, y) = source_velocity(σ_1j, σ_2j, panel_j, x, y) .+ vortex_velocity(γ_1j, γ_2j, panel_j, x, y)

# Evaluation of velocities on panels
source_velocity(σ_1j, σ_2j, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) = panel_velocity(linear_source_velocity_a, linear_source_velocity_b, σ_1j, σ_2j, panel_j, panel_i)
vortex_velocity(γ_1j, γ_2j, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) = panel_velocity(linear_vortex_velocity_a, linear_vortex_velocity_b, γ_1j, γ_2j, panel_j, panel_i)

total_velocity(σ_1j, σ_2j, γ_1j, γ_2j, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) = source_velocity(σ_1j, σ_2j, panel_j, panel_i) .+ vortex_velocity(γ_1j, γ_2j, panel_j, panel_i)

## Matrix construction
#============================================#

linear_vortex_matrix(panels_1 :: Vector{<: AbstractPanel2D}, panels_2 :: Vector{<: AbstractPanel2D}) = two_point_matrix(linear_vortex_velocity_a, linear_vortex_velocity_b, panels_1, panels_2)
linear_source_matrix(panels_1 :: Vector{<: AbstractPanel2D}, panels_2 :: Vector{<: AbstractPanel2D}) = two_point_matrix(linear_source_velocity_a, linear_source_velocity_b, panels_1, panels_2)

kutta(panels_1 :: Vector{<: AbstractPanel2D}, panels_2 :: Vector{<: AbstractPanel2D}) = [ 1 zeros(length(panels_1) - 1)' 1 zeros(length(panels_2) - length(panels_1))' ]
source_influence_matrix(panels_1 :: Vector{<: AbstractPanel2D}, panels_2 :: Vector{<: AbstractPanel2D}) = [ linear_source_matrix(panels_1, panels_2); kutta(panels_1, panels_2) ]
vortex_influence_matrix(panels_1 :: Vector{<: AbstractPanel2D}, panels_2 :: Vector{<: AbstractPanel2D}) = [ linear_vortex_matrix(panels_1, panels_2); kutta(panels_1, panels_2) ]

# Diagonal cases
kutta(panels :: Vector{<: AbstractPanel2D}) = [ 1 zeros(length(panels) - 1)' 1 ]
source_influence_matrix(panels :: Vector{<: AbstractPanel2D}) = [ linear_source_matrix(panels, panels); kutta(panels) ]
vortex_influence_matrix(panels :: Vector{<: AbstractPanel2D}) = [ linear_vortex_matrix(panels, panels); kutta(panels) ]

# Constant-strength sources
source_velocity(σ1, panel :: AbstractPanel2D, x, y) = panel_velocity(constant_source_velocity, σ1, panel, x, y)
source_velocity(σ_j, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) = panel_velocity(constant_source_velocity, σ_j, panel_j, panel_i)

constant_source_influence_coefficient(panel_j, panel_i) = ifelse(panel_i === panel_j, 0.5, influence_coefficient(constant_source_velocity, panel_normal, panel_j, panel_i))

constant_source_matrix(panels :: Vector{<: Panel2D}) = [ constant_source_influence_coefficient(panel_j, panel_i) for (panel_j, panel_i) in product(panels, panels)]

constant_source_boundary_condition(panels :: Vector{<: Panel2D}, u) = -dot.(Ref(u), panel_normal.(panels))

# Generic influence coefficient computation on panels
#============================================#

influence_coefficient(velocity_func :: F, angle_func, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) where F <: Function = dot(panel_velocity(velocity_func, 1., panel_j, panel_i), angle_func(panel_i))

neumann_influence_matrix(func, angle_func, panel_is :: Vector{<: AbstractPanel2D}, panel_js :: Vector{<: AbstractPanel2D}) = [ influence_coefficient(func, angle_func, panel_j, panel_i) for (panel_i, panel_j) in product(panel_is, panel_js) ]

function two_point_matrix(vel_a :: F, vel_b :: G, panels_1 :: Vector{<: AbstractPanel2D}, panels_2 :: Vector{<: AbstractPanel2D}) where F where G
	N 	  = min(length(panels_1), length(panels_2))
	inf_a = neumann_influence_matrix(vel_a, panel_normal, panels_1, panels_2)
	inf_b = neumann_influence_matrix(vel_b, panel_normal, panels_1, panels_2)
	[ inf_a zeros(N) ] + [ zeros(N) inf_b ]
end

neumann_boundary_condition(panels :: Vector{<: AbstractPanel2D}, u) = [ -dot.(Ref(u), panel_normal.(panels)); 0 ]

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
module linux
	using sys

	@pkgconf begin
		using libc.POSIX

		c``
		c"""
			#include <asm-generic/types.h>
			#include <asm/types.h>
			#include <linux/types.h>
			// add other linux/*.h here...
		"""ji
	end
end


text:
function UnivariatePeriodicInterpolator(arg0::UnivariateInterpolator, arg1::jdouble)
    return UnivariatePeriodicInterpolator((UnivariateInterpolator, jdouble), arg0, arg1)
end

function UnivariatePeriodicInterpolator(arg0::UnivariateInterpolator, arg1::jdouble, arg2::jint)
    return UnivariatePeriodicInterpolator((UnivariateInterpolator, jdouble, jint), arg0, arg1, arg2)
end

function interpolate(obj::UnivariatePeriodicInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "interpolate", UnivariateFunction, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end


text:
dict = Dict('a' => 97, 'b' => 98) # list keys/values
# Dict{Char,Int64} with 2 entries:
#   'b' => 98
#   'a' => 97

dict = Dict(c => Int(c) for c = 'a':'d') # dict comprehension
# Dict{Char,Int64} with 4 entries:
#   'b' => 98
#   'a' => 97
#   'd' => 100
#   'c' => 99

dict['é'] = 233; dict # add an element
# Dict{Char,Int64} with 3 entries:
#   'b' => 98
#   'a' => 97
#   'é' => 233

emptydict = Dict() # create an empty dict
# Dict{Any,Any} with 0 entries

dict["a"] = 1 # type mismatch
# ERROR: MethodError: Cannot `convert` an object of type String to an object of type Char

typeof(dict) # type is infered correctly
# Dict{Char,Int64}

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function solve_nde(nde, NN, T₀, alg, nde_params)
    nn_weights, _ = Flux.destructure(NN)
    # See: https://github.com/SciML/DiffEqFlux.jl/blob/449efcecfc11f1eab65d0e467cf57db9f5a5dbec/src/neural_de.jl#L68
    return solve(nde, alg, reltol=1e-4, u0=T₀, p=[nn_weights; nde_params],
                 sense=InterpolatingAdjoint(autojacvec=ZygoteVJP(), checkpointing=true))
end

function solve_nde(ds, NN, NDEType, algorithm, T_scaling, wT_scaling; T₀=nothing)

    T = ds["T"]
    wT = ds["wT"]
    Nz = size(T, 3)
    zc = znodes(T)
    zf = znodes(wT)

    H = abs(zf[1]) # Domain depth/height

    Δẑ = diff(zc)[1] / H  # Non-dimensional grid spacing
    Dzᶠ = Dᶠ(Nz, Δẑ) # Differentiation matrix operator

    nde_params = FreeConvectionNDEParameters(ds, T_scaling, wT_scaling)
    nde = NDEType(NN, ds)

    if isnothing(T₀)
        T₀ = T_scaling.(interior(T)[1, 1, :, 1])
    end

    T = solve_nde(nde, NN, T₀, algorithm, nde_params) |> Array

    bottom_flux, top_flux, _ = nde_params

    Nz, Nt = size(T)
    wT = zeros(Nz+1, Nt)

    for n in 1:Nt
        T_n = T[:, n]

        wT_NN_interior = NN(T_n)
        wT_NN_n = cat(bottom_flux, wT_NN_interior, top_flux, dims=1)

        if NDEType == FreeConvectionNDE
            @. wT[:, n] = wT_NN_n
        elseif NDEType == ConvectiveAdjustmentNDE
            ∂T∂z_n = Dzᶠ * T_n
            K∂T∂z_n = min.(0, 10 * ∂T∂z_n)
            @. wT[:, n] = wT_NN_n - K∂T∂z_n
        end
    end

    return (T=inv(T_scaling).(T), wT=inv(wT_scaling).(wT))
end


text:
## General test for orthogonal polynomials
import SpecialPolynomials: domain

T = Float64
Ps = (
    Chebyshev,
    OrthonormalChebyshev,
    ChebyshevU,
    Laguerre{0},
    Laguerre{1 / 2},
    OrthonormalLaguerre{1 / 2},
    Hermite,
    ChebyshevHermite,
    Jacobi{1 / 2,1 / 2},
    Jacobi{1 / 4,3 / 4},
    Jacobi{-1 / 2,1 / 2},
    Jacobi{1 / 2,-1 / 2},
    OrthonormalJacobi{1 / 2,1 / 2},
    Legendre,
    OrthonormalLegendre,
    Gegenbauer{1 / 2},
    OrthonormalGegenbauer{1 / 2},
    Bessel{3 / 2}, # Bessel{1} is an issue
    Bessel{1 / 2},
)
DPs = (Charlier{1 / 2}, Meixner{1 / 2,1 / 2}, Krawchouk{1 / 2,10}, Hahn{1 / 4,1 / 2,10})

@testset "Construction" begin
    for P in union(Ps, DPs)
        x = variable(P)
        # basic recurrence

        for n in 1:6
            @test basis(P, n + 1) ≈
                  (SP.An(P, n) * x + SP.Bn(P, n)) * basis(P, n) -
                  SP.Cn(P, n) * basis(P, n - 1)
        end

        # leading term and ratios
        x = variable(Polynomial)
        for n in 0:5
            @test basis(P, n)(x)[end] ≈ SP.kn(P, n)
        end
        for n in 0:5
            @test SP.k1k0(P, n) ≈ SP.kn(P, n + 1) / SP.kn(P, n)
        end
        for n in 1:5
            @test SP.k1k_1(P, n) ≈ SP.kn(P, n + 1) / SP.kn(P, n - 1)
        end

        # basis
        x = variable(P)
        p0, p1, p2, p3, p4, p5, p6 = ps = basis.(P, 0:6)
        p = P([1, 2, 3, 4, 5, 6])
        for x in range(max(-1, first(domain(P))), stop=min(1, last(domain(P))), length=10)
            @test p(x) ≈ 1p0(x) + 2p1(x) + 3p2(x) + 4p3(x) + 5p4(x) + 6p5(x)
        end

        ## Different type
        T = Float64
        @test eltype(P(ones(Int32, 4))) == Int32
        @test eltype(P{T}(ones(Int32, 4))) == T
        @test degree(P(1) - 1) == -1
        @test degree(zero(P)) == -1
        @test degree(one(P)) == 0
        @test degree(variable(P)) == 1

        # variable and P() constructor for `x` in basis
        @test degree(variable(P)) == 1
        @test variable(P)(1) ≈ 1
        @test degree(P()) == 1
        @test P()(1) ≈ 1
        @test variable(P, :y) == P(:y)

        # test  ignore constant's symbol
        @test zero(P, :x) == zero(P, :y)
        @test one(P, :x) == one(P, :y)
        @test !(variable(P, :x) == variable(P, :y))

        @test !(zero(P, :x) === zero(P, :y))
        @test !(one(P, :x) === one(P, :y))
        @test !(variable(P, :x) === variable(P, :y))

        @test zero(P, :x) ≈ zero(P, :y)
        @test one(P, :x) ≈ one(P, :y)
        @test (variable(P, :x) ≈ variable(P, :x))
        @test_throws ArgumentError variable(P, :x) ≈ variable(P, :y)
    end
end

@testset "Structural equations" begin
    for P in Ps
        P <: Hermite && continue #  Hermite  has issues
        for i in 2:5
            ps = basis.(P, (i + 1):-1:(i - 1))
            dps = derivative.(ps)
            pᵢ, dpᵢ = ps[2], dps[2]

            a, b, c, d, e = SP.abcde(P)
            x = variable(P)
            σ = a * x^2 + b * x + c

            # x⋅p_n   = [an, bn, cn]    ⋅ [p_{n+1}, p_n, p_{n-1}]     #  Eqn (7)
            as = SP.an(P, i), SP.bn(P, i), SP.cn(P, i)
            @test x * pᵢ ≈ sum(a * p for (a, p) in zip(as, ps))

            # σ⋅p'_n  = [αn, βn, γn]    ⋅  [p_{n+1}, p_n, p_{n-1}]    # Eqn (9), n ≥ 1
            αs = SP.αn(P, i), SP.βn(P, i), SP.γn(P, i)
            @test σ * dpᵢ ≈ sum(a * p for (a, p) in zip(αs, ps))

            # p_n    = [ân, b̂n, ĉn]    ⋅  [p'_{n+1}, p'_n, p'_{n-1}] # Eqn (19)
            âs = SP.ân(P, i), SP.b̂n(P, i), SP.ĉn(P, i)
            @test pᵢ ≈ sum(a * dp for (a, dp) in zip(âs, dps))

            # x⋅p'_n  = [αᴵn, βᴵn, γᴵn] ⋅  [p'_{n+1}, p'_n, p'_{n-1}] # Eqn  (14) with  α^*, β^*,  γ^*
            @test (a=a, b=b, c=c, d=d + 2a, e=e + b) == SP.abcdeᴵ(P)
            αᴵs = SP.αᴵn(P, i), SP.βᴵn(P, i), SP.γᴵn(P, i)
            @test x * dpᵢ ≈ sum(a * dp for (a, dp) in zip(αᴵs, dps))
        end
    end

    for P in DPs
        for i in 1:5
            ps = basis.(P, (i + 1):-1:(i - 1))
            dps = SP.:∇ₓ.(ps)
            Δps = SP.:Δₓ.(ps)
            pᵢ, dpᵢ, Δpᵢ = ps[2], dps[2], Δps[2]

            a, b, c, d, e = SP.abcde(P)
            x = variable(P)
            σ = a * x^2 + b * x + c

            # x⋅p_n   = [an, bn, cn]    ⋅ [p_{n+1}, p_n, p_{n-1}]     #  Eqn (7)
            as = SP.an(P, i), SP.bn(P, i), SP.cn(P, i)
            @test x * pᵢ ≈ sum(a * p for (a, p) in zip(as, ps))

            αs = SP.αn(P, i), SP.βn(P, i), SP.γn(P, i)
            @test σ * dpᵢ ≈ sum(a * p for (a, p) in zip(αs, ps))

            âs = SP.ân(P, i), SP.b̂n(P, i), SP.ĉn(P, i)  #  Eqn 20
            @test pᵢ ≈ sum(a * dp for (a, dp) in zip(âs, Δps))

            αᴵs = SP.αᴵn(P, i), SP.βᴵn(P, i), SP.γᴵn(P, i)   # Eqn (16)
            @test x * Δpᵢ ≈ sum(a * dp for (a, dp) in zip(αᴵs, Δps))
        end
    end
end

@testset "Conversion" begin
    for (PP, Q) in ((Ps, Polynomial), (DPs, FallingFactorial))
        for P in PP
            for _ in 1:10
                ps = rand(1.0:8, 5)
                p = P(ps)
                @test convert(P, convert(Q, p)) ≈ p
            end
        end
    end
    # through evaluation
    for PP in (Ps, DPs)
        for P in PP
            for Q in PP
                p = P([1, 2, 3])
                x = variable(Q)
                @test p(x) ≈ convert(Q, p)
            end
        end
    end

    x = variable(Polynomial)
    # Connections: groupings have same sigma = a⋅x² + b⋅x + c
    for Ps in (
        (
            Chebyshev,
            ChebyshevU,
            Gegenbauer{1 / 2},
            Gegenbauer{1},
            Jacobi{1 / 4,1 / 4},
            Jacobi{3 / 4,3 / 4},
        ),
        (Laguerre{0}, Laguerre{1 / 2}, Laguerre{1}),
        (Bessel{1 / 4}, Bessel{1 / 2}, Bessel{3 / 4}, Bessel{2}),
    )
        for P in Ps
            for Q in Ps
                for n in 2:5
                    q = basis(Q, n)
                    # various ways to convert to P

                    @test q(variable(P))(variable()) ≈ q(variable())  # evaluation
                    @test SP._convert_cop(P, q)(variable()) ≈ q(variable())  # structural  equations
                    @test convert(P, convert(Polynomial, q))(variable()) ≈ q(variable()) # through Standard  basis
                    @test convert(P, q)(x) ≈ q(x)
                end
            end
        end
    end

    # Connections via FastTransforms
    x = variable()
    for P in Ps
        for Q in Ps
            p = P([1.0, 2.0, 3.0])
            @test convert(Q, p)(x) ≈ p(x)
        end
    end
end

@testset "Evaluation" begin
    for P in Ps
        (SP.ismonic(P) || SP.isorthonormal(P)) && continue
        n, x = 4, 1.23
        p = basis(P, n)
        # compares  clenshaw  to  hypergeometric
        @test p(x) ≈ Basis(P, n)(x)
    end
end

@testset "Arithmetic" begin
    x = variable(Polynomial{Float64})

    for P in union(Ps, DPs)
        pNULL = zero(P)
        p0, p1, p2, p3, p4 = ps = basis.(P, 0:4)

        as = [1, 2, 3, 4, 5]
        @test sum(a * pᵢ for (a, pᵢ) in zip(as, ps)) ≈ P(as)

        x = variable(Polynomial)
        for p in ps
            for q in ps
                for op in (+, *)
                    @test op(p, q)(x) ≈ op(p(x), q(x))
                end
            end
        end

        for _ in 1:10
            p, q = P(rand(5)), P(rand(5))
            for op in (+, *)
                @test (op(p, q))(x) ≈ op(p(x), q(x))
            end
        end
    end
end

@testset "Elementwise operations" begin
    for P in union(Ps, DPs)
        xs = [1, 2, 3, 4]
        p = P(xs)
        @test -p == P(-xs)
        ys = copy(xs ./ 1) # might need float
        ys[1] += one(P)[0]
        @test p + 1 ≈ p + P(1)
        @test 2p == P(2xs)
    end
end

# too slow
@testset "Orthogonality" begin
    n = 5
    for P in Ps
        a, b = extrema(P)
        (isinf(a) || isinf(b)) && continue
        for i in 2:n
            for j in (i + 1):n
                val = SP.innerproduct(P, Polynomials.basis(P, i), Polynomials.basis(P, j))
                @test abs(val) <= 1e-4
            end
        end
    end

    for P in (Krawchouk{1 / 2,n}, Hahn{1 / 4,1 / 2,n})
        for i in 2:n
            for j in (i + 1):n
                val = SP.innerproduct(P, Polynomials.basis(P, i), Polynomials.basis(P, j))
                @test abs(val) <= 1e-12
            end
        end
    end
end

@testset "divrem" begin
    for P in union(Ps, DPs)
        ps, qs = rand(1:6, 4), rand(1:6, 3)
        p, q = P(ps), P(qs)
        a, b = divrem(p, q)
        @test p ≈ q * a + b
    end
end

@testset "Derivatives and integrals" begin
    for P in Ps
        Q = Polynomial
        x = variable(Q)

        for i in 1:6
            p = basis(P, i)
            @test derivative(p)(x) ≈ derivative(p(x))
        end

        p = P([1.0, 2, 3, 4, 3, 2, 1])
        q = convert(Polynomial, p)
        @test maximum(
            abs,
            [
                integrate(p, a, a + 1 / 2) - integrate(q, a, a + 1 / 2) for
                a in range(0, stop=1 / 2, length=10)
            ],
        ) <= sqrt(eps())

        for _ in 1:10
            ps = rand(1:10, 5)
            p = P(ps)
            @test (derivative ∘ derivative)(p) ≈ derivative(p, 2)

            q = (derivative ∘ integrate)(p) - p
            q = chop(q, atol=1e-8)
            @test degree(q) <= 0

            q = (integrate ∘ derivative)(p) - p
            q = chop(q, atol=1e-8)
            @test degree(q) <= 0
        end
    end

    # special  case 𝐍  cases
    𝐍 = 9
    for P in (Krawchouk{1 / 2,𝐍}, Krawchouk{1 / 4,𝐍}, Hahn{1 / 2,1 / 2,𝐍})
        for i in 0:𝐍
            p = basis(P, i)
            (p(𝐍) - p(0)) <= 1e-12 && continue
            @test sum(SP.∇ₓ(p)(j) for j in 1:𝐍) ≈ p(𝐍) - p(0)
            @test sum(SP.Δₓ(p)(j) for j in 0:(𝐍 - 1)) ≈ p(𝐍) - p(0)
        end
    end
end

@testset "roots" begin
    for P in union(Ps, DPs)
        P <: Bessel && continue
        for n in 6:2:10
            rts = roots(convert(Polynomial, basis(P, n)))
            evals = eigvals(SP.jacobi_matrix(P, n))
            @test maximum(abs, rts - evals) <= 1e-4
        end
    end

    # comrade matrix approach compared to conversion
    for P in
        (Legendre, MonicLegendre, Chebyshev, MonicChebyshev, ChebyshevU, MonicChebyshevU)
        for T in (Float64, Complex{Float64})
            ps = vcat(rand(T, 4), 1)
            p = P(ps)
            q = convert(Polynomial, p)
            λs = sort(norm.(roots(p)))
            γs = sort(norm.(roots(q)))
            @test maximum(abs, λs - γs) < sqrt(eps())
        end
    end

    # comrade pencil
    for P in (Legendre, Chebyshev, ChebyshevU)
        p = P(rand(5))
        λ = rand()
        C₀, C₁ = SP.comrade_pencil(p)
        M = λ * C₁ - C₀
        L, Ũ = SP.comrade_pencil_LU(p)(λ)
        @test det(Ũ) ≈ 1
        Ũ[end, end] *= p(λ)
        @test L * Ũ ≈ M
    end
end

@testset "fitting" begin
    f(x) = exp(-x) * cospi(x)

    # Interpolating
    for P in Ps
        P <: SP.AbstractOrthogonalPolynomial || continue
        dom = domain(P)
        (isinf(first(dom)) || isinf(last(dom))) && continue
        q = fit(P, f, 10)
        @test maximum(abs, q(x) - f(x) for x in range(0, stop=1 / 2, length=10)) <= 1e-4
    end

    # least squares
    for P in (Chebyshev,)
        q = fit(Val(:lsq), P, f, 10)
        @test maximum(abs, q(x) - f(x) for x in range(0, stop=1 / 2, length=10)) <= 1e-4
    end

    # series
    for P in (Chebyshev,)
        q = fit(Val(:series), P, f)
        @test maximum(abs, q(x) - f(x) for x in range(0, stop=1 / 2, length=10)) <=
              sqrt(eps())
    end
end

@testset "quadrature" begin
    f(x) = x^7
    n = 4

    for P in Ps
        P <: SP.AbstractOrthogonalPolynomial || continue
        !all(isfinite.(extrema(P))) && continue
        q = sum(f(tau) * w for (tau, w) in zip(SP.gauss_nodes_weights(P, n)...))
        p = SP.innerproduct(P, f, one)
        @test abs(p - q) <= 10sqrt(eps(T))
    end
end

@testset "Showing" begin
    @test sprint(show, Hermite([1, -2, 3])) == "Hermite(1⋅H₀(x) - 2⋅H₁(x) + 3⋅H₂(x))"

    ps = [1.0 + 2.0im, -0.5, 3.7im, 0.1 + 0.01im]     # issue 26 with complex coefficients
    p = Hermite(ps)
    @test sprint(show, p) ==
          "Hermite((1.0 + 2.0im)⋅H₀(x) - 0.5⋅H₁(x) + 3.7im⋅H₂(x) + (0.1 + 0.01im)⋅H₃(x))"
end


text:
struct ConstantLayer
    w
end

ConstantLayer(N::Integer) = ConstantLayer(zeros(N))

(m::ConstantLayer)(x) = m.w

Flux.@functor ConstantLayer

mutable struct A2TNetwork
  base::Chain
  attn::Chain
  solutions::Array{Function} # function takes in state and outputs vector of action values
end

function (m::A2TNetwork)(input)
    b = m.base(input) #output is (Na, b)
    w = m.attn(input) #output is (Nt+1, b)
    B, Nt = size(input, 2), size(w,1) - 1
    qs = Zygote.ignore(() -> [hcat([s(input[:,i]) for s in m.solutions]...) for i=1:B]) #output is (Na, Nt)
    Flux.stack(qs .* Flux.unstack(w[1:Nt, :], 2), 2) .+ w[Nt+1:Nt+1, :] .* b
end

Flux.@functor A2TNetwork

Flux.trainable(m::A2TNetwork) = (m.base, m.attn)

function Base.iterate(m::A2TNetwork, i=1)
    i > length(m.base.layers) + length(m.attn.layers) && return nothing
    if i <= length(m.base.layers)
        return (m.base[i], i+1)
    elseif i <= length(m.base.layers) + length(m.attn.layers)
        return (m.attn[i - length(m.base.layers)], i+1)
    end
end

function Base.deepcopy(m::A2TNetwork)
  A2TNetwork(deepcopy(m.base), deepcopy(m.attn), deepcopy(m.solutions))
end


text:
module SoftConfidenceWeighted

using Devectorize

import Distributions: Normal, cdf
import SVMLightLoader: SVMLightFile

export init, fit!, predict, SCW1, SCW2

typealias AA AbstractArray

@enum SCWType SCW1 SCW2

type CDF
    phi
    psi
    zeta

    function CDF(ETA)
        phi = cdf(normal_distribution, ETA)
        psi = 1 + phi^2 / 2
        zeta = 1 + phi^2
        new(phi, psi, zeta)
    end
end

#calc cdf in a constructor
type SCW
    C::Float64
    cdf::CDF
    ndim::Int64
    weights::Array{Float64, 1}
    covariance::Array{Float64, 1}
    has_fitted::Bool

    function SCW(C, ETA)
        new(C, CDF(ETA), -1, [], [], false)
    end
end

function set_dimension!(scw::SCW, ndim::Int)
    assert(!scw.has_fitted)
    scw.ndim = ndim
    scw.weights = zeros(ndim)
    scw.covariance = ones(ndim)
    scw.has_fitted = true
    return scw
end

normal_distribution = Normal(0, 1)

function calc_margin{T<:AA,R<:Real}(scw::SCW, x::T, label::R)
    #Devectorize.jl requires assignment
    @devec t = label .* dot(scw.weights, x)
end

function calc_confidence{T<:AA}(scw::SCW, x::T)
    @devec t = dot(x, scw.covariance .* x)
end

function calc_alpha1{T<:AA,R<:Real}(scw::SCW, x::T, label::R)
    v = calc_confidence(scw, x)
    m = calc_margin(scw, x, label)
    cdf = scw.cdf
    (phi, psi, zeta) = (cdf.phi, cdf.psi, cdf.zeta)

    j = m^2 * phi^4 / 4
    k = v * zeta * phi^2
    t = (-m*psi + sqrt(j+k)) / (v*zeta)
    return min(scw.C, max(0, t))
end

function calc_alpha2{T<:AA,R<:Real}(scw::SCW, x::T, label::R)
    v = calc_confidence(scw, x)
    m = calc_margin(scw, x, label)
    cdf = scw.cdf
    (phi, psi, zeta) = (cdf.phi, cdf.psi, cdf.zeta)

    n = v+1 / scw.C
    a = (phi*m*v)^2
    b = 4*n*v * (n + v * phi^2)
    gamma = phi * sqrt(a+b)

    c = -(2 * m * n + m * v * phi^2)
    d = n^2 + n * v * phi^2
    t = (c+gamma) / 2d
    return max(0, t)
end

function init(;C = 1, ETA = 1, type_ = SCW1::SCWType)
    global calc_alpha
    if type_ == SCW1
        calc_alpha = calc_alpha1
    elseif type_ == SCW2
        calc_alpha = calc_alpha2
    else
        assert(true)
    end

    return SCW(C, ETA)
end

function loss{T<:AA,R<:Real}(scw::SCW, x::T, label::R)
    @devec t = label .* dot(scw.weights, x)
    if t >= 1
        return 0
    end
    return 1-t
end

function calc_beta{T<:AA,R<:Real}(scw::SCW, x::T, label::R)
    alpha = calc_alpha(scw, x, label)
    v = calc_confidence(scw, x)
    m = calc_margin(scw, x, label)
    cdf = scw.cdf
    (phi, psi, zeta) = (cdf.phi, cdf.psi, cdf.zeta)

    j = -alpha * v * phi
    k = sqrt((alpha*v*phi)^2 + 4v)
    u = (j+k)^2 / 4
    return (alpha * phi) / (sqrt(u) + v*alpha*phi)
end

function update_covariance!{S<:AA,T<:AA,R<:Real}(t::S, scw::SCW, x::T, label::R)
    beta = calc_beta(scw, x, label)
    c = scw.covariance

    # same as
    # scw.covariance -= beta * (c .* x) .* (c .* x)
    @devec t[:] = (c .* x) .* (c .* x)
    BLAS.axpy!(-beta, t, scw.covariance)
    return scw
end

function update_weights!{S<:AA,T<:AA,R<:Real}(t::S, scw::SCW, x::T, label::R)
    alpha = calc_alpha(scw, x, label)

    # same as
    # scw.weights += alpha * label * (scw.covariance .* x)
    @devec t[:] = scw.covariance .* x
    BLAS.axpy!(alpha * label, t, scw.weights)
    return scw
end

function update!{S<:AA,T<:AA,R<:Real}(t::S, scw::SCW, x::T, label::R)
    x = vec(full(x))
    if loss(scw, x, label) > 0
        update_weights!(t, scw, x, label)
        update_covariance!(t, scw, x, label)
    end
    return scw
end

function train!{T<:AA,R<:AA}(scw::SCW, X::T, labels::R)
    # preallocate for performance optimization
    t = Array(Float64, size(X, 1))
    for i in 1:size(X, 2)
        label = labels[i]

        if label != 1 && label != -1
            throw(ArgumentError("Label must be 1 or -1"))
        end

        update!(t, scw, slice(X, :, i), label)
    end
    return scw
end

function fit!{T, R}(scw::SCW, X::AA{T, 2}, labels::AA{R, 1})
    if size(X, 2) != length(labels)
        message = "X and y have incompatible shapes. " *
                  "size(X) = $(size(X)) size(y) = $(size(y))"
        throw(ArgumentError(message))
    end

    if !scw.has_fitted
        ndim = size(X, 1)
        set_dimension!(scw, ndim)
    end

    train!(scw, X, labels)
    return scw
end

function fit!{T}(scw::SCW, x::AA{T, 1}, label::Real)
    ndim = size(x, 1)
    if !scw.has_fitted
        set_dimension!(scw, ndim)
    end

    t = Array(Float64, ndim)
    update!(t, scw, x, label)
end

function fit!(scw::SCW, filename::AbstractString, ndim::Int64)
    if !scw.has_fitted
        set_dimension!(scw, ndim)
    end

    t = Array(Float64, ndim)
    for (vector, label) in SVMLightFile(filename, ndim)
        update!(t, scw, vector, label)
    end
    return scw
end

@deprecate fit fit!

function compute{T<:AA}(scw::SCW, x::T)
    x = vec(full(x))
    if dot(x, scw.weights) > 0
        return 1
    else
        return -1
    end
end

function predict{T<:AA}(scw::SCW, X::T)
    N = size(X, 2)
    labels = Array(Int, N)
    for i in 1:size(X, 2)
        labels[i] = compute(scw, slice(X, :, i))
    end
    return labels
end

function predict(scw::SCW, filename::AbstractString)
    labels = Int64[]
    for (x, _) in SVMLightFile(filename, scw.ndim)
        label = compute(scw, x)
        push!(labels, label)
    end
    return labels
end

end # module


text:
function generate_ccs9_to_icd9(download_cache = ensure_downloaded_files())
    ccs_filename = joinpath(download_cache, "ccs9.txt")
    ccs_filecontents = strip(read(ccs_filename, String)::String)
    elements = strip.(split(ccs_filecontents, "\n\n"))
    ccs9_to_icd9            = Dict{Class, Set{Class}}()
    ccs9_to_ccs_description = Dict{Class, String}()
    for element in elements
        if startswith(element, "Appendix A - ")
            continue
        end
        if startswith(element, "Revised ")
            continue
        end
        element2 = replace(element, '`' => ''')
        r = r"^(\d\d*?)\s\s*?([A-Za-z() :;,\-'\/]*?)\s*?\n([\S\s]*?)$"
        m = match(r, element2)
        (m === nothing) && @error("Could not parse element", element2)
        ccs_value       = strip(m[1])
        ccs_description = strip(m[2])
        ccs_node = construct_ccs9(ccs_value)
        icd_values = strip.(split(strip(replace(m[3], '\n' => ' ')), ' '))
        filter!(x -> !isempty(x), icd_values)
        icd_nodes = construct_icd9.(icd_values)
        icd_set = Set{Class}(icd_nodes)
        setindex_no_overwrite!(
            ccs9_to_icd9,
            ccs_node,
            icd_set,
        )
        setindex_no_overwrite!(
            ccs9_to_ccs_description,
            ccs_node,
            ccs_description,
        )
    end
    return (; ccs9_to_icd9, ccs9_to_ccs_description)
end

function generate_ccs10_to_icd10(download_cache = ensure_downloaded_files())
    ccs_filename = joinpath(download_cache, "ccs10.csv")
    ccs_file = CSV.File(ccs_filename)
    all_column_names_str = String.(ccs_file.names)
    unique!(all_column_names_str)
    sort!(all_column_names_str)
    ccs_column_names_str = filter(
        x -> startswith(x, "'CCSR CATEGORY ") && !occursin("DESCRIPTION", x),
        all_column_names_str,
    )
    unique!(ccs_column_names_str)
    sort!(ccs_column_names_str)
    remove_single_quotes = x -> begin
        r = r"^'(.*?)'$"
        m = match(r, x)
        (m === nothing) && @error("Could not parse", x)
        return convert(String, m[1])::String
    end
    ccscolname_to_ccsdesccolname = x -> begin
        return "'" * remove_single_quotes(x) * " DESCRIPTION'"
    end
    (length(ccs_column_names_str) < 1) && throw(ErrorException("Not enough CCS category columns"))

    ccs10_to_icd10           = Dict{Class, Set{Class}}()
    ccs10_to_ccs_description = Dict{Class, String}()
    for row in ccs_file
        icd_value = remove_single_quotes(row[Symbol("'ICD-10-CM CODE'")])
        icd_node  = construct_icd10(icd_value)
        for ccscolname in ccs_column_names_str
            ccs_value = row[Symbol(ccscolname)] |> remove_single_quotes |> strip
            if isempty(ccs_value)
                continue
            end
            ccs_node = construct_ccs10(ccs_value)
            if !haskey(ccs10_to_icd10, ccs_node)
                ccs10_to_icd10[ccs_node] = Set{Class}()
            end
            push!(ccs10_to_icd10[ccs_node], icd_node)

            ccsdesccolname = ccscolname_to_ccsdesccolname.(ccscolname)
            ccs_description = row[Symbol(ccsdesccolname)] |> strip
            ccs_description_string = convert(String, ccs_description)::String
            setindex_same_value!(
                ccs10_to_ccs_description,
                ccs_node,
                ccs_description_string,
            )
        end
    end

    return (; ccs10_to_icd10, ccs10_to_ccs_description)
end


text:
using OceanParameterizations
using WindMixing
using Plots

reconstruct_fluxes = true
println("Reconstruct fluxes? $(reconstruct_fluxes)")

enforce_surface_fluxes = true
println("Enforce surface fluxes? $(enforce_surface_fluxes)")

subsample_frequency = 4
println("Subsample frequency for training... $(subsample_frequency)")

file_labels = Dict(
    "free_convection" => "Free convection",
    "strong_wind" => "Strong wind",
    "strong_wind_no_coriolis" => "Strong wind, no rotation",
    "weak_wind_strong_cooling" => "Weak wind, strong cooling",
    "strong_wind_weak_cooling" => "Strong wind, weak cooling",
    "strong_wind_weak_heating" => "Strong wind, weak heating"
)

logγ_range = collect(-1.0:0.05:1.5)
n = length(logγ_range)

# Rows correspond to kernel functions, cols correspond to logγ values
errors_uw = zeros(4,n)
errors_vw = zeros(4,n)
errors_wT = zeros(4,n)

## Gaussian Process Regression
# Find the kernel that minimizes the prediction error on the training data
# * Sweeps over length-scale hyperparameter value in logγ_range

files =  ["free_convection", "strong_wind", "strong_wind_no_coriolis",
            "weak_wind_strong_cooling", "strong_wind_weak_cooling", "strong_wind_weak_heating"]
output_directory="GP/subsample_$(subsample_frequency)/reconstruct_$(reconstruct_fluxes)/enforce_surface_fluxes_$(enforce_surface_fluxes)"
# isdir(dirname(output_directory)) ||
mkpath(output_directory)

for i=1:length(files)
    # Train on all except file i
    train_files = files[1:end .!= i]
    𝒟train = data(train_files,
                        scale_type=ZeroMeanUnitVarianceScaling,
                        reconstruct_fluxes=reconstruct_fluxes,
                        subsample_frequency=subsample_frequency,
                        enforce_surface_fluxes=enforce_surface_fluxes)
    # Test on file i
    test_file = files[i]
    𝒟test = data(test_file,
                        override_scalings=𝒟train.scalings, # use the scalings from the training data
                        reconstruct_fluxes=reconstruct_fluxes,
                        subsample_frequency=subsample_frequency,
                        enforce_surface_fluxes=enforce_surface_fluxes)

    for k=1:4
        errors_uw[k,:] .+= error_per_gamma(𝒟train.uw, 𝒟test.uw, k; logγ_range=logγ_range)
        errors_vw[k,:] .+= error_per_gamma(𝒟train.vw, 𝒟test.vw, k; logγ_range=logγ_range)
        errors_wT[k,:] .+= error_per_gamma(𝒟train.wT, 𝒟test.wT, k; logγ_range=logγ_range)
    end
end

# total error across the 6 test files (where you train on all but that test file) at each log(γ) value
function hyperparameter_landscapes(errors)
    myplot(i) = plot(logγ_range, errors[i,:], linewidth=4, c=:purple, xlabel="log(γ)", ylabel="Error", legend=false, yscale=:log10)
    layout = @layout [a; b; c; d]
    p = plot(myplot(1), myplot(2), myplot(3), myplot(4), layout=layout)
end

# From top to bottom the rows in the plot correspond to:
# "Squared Exponential Kernel"
# "Matérn Kernel ν=1/2"
# "Matérn Kernel ν=3/2"
# "Matérn Kernel ν=5/2"
# From left to right the columns in the plot correspond to uw, vw, wT
p1 = hyperparameter_landscapes(errors_uw)
p2 = hyperparameter_landscapes(errors_vw)
p3 = hyperparameter_landscapes(errors_wT)
layout = @layout [a b c]
p = plot(p1, p2, p3, layout=layout)
savefig(p, output_directory*"/hyperparameter_landscapes_uw_vw_wT.pdf")

# Extract the optimal kernel function and log(γ) value from error matrix
k_logγ(errors) = (argmin(errors)[1], logγ_range[argmin(errors)[2]])
k_uw, logγ_uw = k_logγ(errors_uw)
k_vw, logγ_vw = k_logγ(errors_vw)
k_wT, logγ_wT = k_logγ(errors_wT)

# Create kernel objects
uw_kernel = get_kernel(k_uw, logγ_uw, 0.0, euclidean_distance)
vw_kernel = get_kernel(k_vw, logγ_vw, 0.0, euclidean_distance)
wT_kernel = get_kernel(k_wT, logγ_wT, 0.0, euclidean_distance)
println(uw_kernel)
println(vw_kernel)
println(wT_kernel)

i=1
train_files = files[1:end .!= i]

subsample_frequency
𝒟train = data(train_files,
                    scale_type=ZeroMeanUnitVarianceScaling,
                    reconstruct_fluxes=reconstruct_fluxes,
                    subsample_frequency=subsample_frequency,
                    enforce_surface_fluxes=enforce_surface_fluxes)
test_file = files[i]
𝒟test = data(test_file,
                    override_scalings=𝒟train.scalings, # use the scalings from the training data
                    reconstruct_fluxes=reconstruct_fluxes,
                    subsample_frequency=subsample_frequency,
                    enforce_surface_fluxes=enforce_surface_fluxes)

# Trained GP models
uw_GP_model = gp_model(𝒟train.uw, uw_kernel)
vw_GP_model = gp_model(𝒟train.vw, vw_kernel)
wT_GP_model = gp_model(𝒟train.wT, wT_kernel)

# GP predictions on test data
uw_GP = predict(𝒟test.uw, uw_GP_model)
vw_GP = predict(𝒟test.vw, vw_GP_model)
wT_GP = predict(𝒟test.wT, wT_GP_model)

# Compare GP predictions to truth
myanimate(xs, name) = animate_prediction(xs, name, 𝒟test, test_file;
                        filename=name*"_optimal_kernel_showing_free_convection", legend_labels=["GP(u,v,T)","Truth"], directory=output_directory)
myanimate(uw_GP, "uw")
myanimate(vw_GP, "vw")
myanimate(wT_GP, "wT")


text:
"""
    alternative_model_scVI(n_epochs, lr, n_genes, n_cells)
Trains a scVI model and extracts the contribution to the evidence
lower bound for each cell. If `reference=true` then scVI will be randomly
initialized and not trained to be used as a reference.
`n_epochs` number of epochs for training scVI when reference=false.
`n_genes` number of genes to read from dataset.
`n_cells` number of cells to subsample from dataset.
`reference` whether to initialize scVI randomly or train the model.
"""
function alternative_model_scVI(;
    data_path::String,
    n_epochs::Int64,
    lr::Float64,
    n_genes::Int64,
    n_cells::Int64,
    n_latent::Int64,
    n_hidden::Int64,
    n_layers::Int64,
    batch_size::Int64,
    pretrained_model::String="",
    subsampling::Bool=false,
    use_autotune::Bool=false,
    search_space)

    if isempty(pretrained_model) == false
        countmatrix = scvi_dataset.CsvDataset(data_path, save_path = "", new_n_genes = n_genes);

        if subsampling == true
            Random.seed!(111);
            subsample_scVI_cells!(countmatrix, n_cells);
        end;

        # Load pre-trained model
        trainer = torch.load(pretrained_model);

        full = trainer.create_posterior(trainer.model, countmatrix);
        latent, batch_indices, labels = full.sequential().get_latent(give_mean=true);

        elbo_per_cell_null = full.elbo_sample();
        elbo_null = full.elbo();

        random.seed(111);
        torch.manual_seed(111);
        trainer.show_progbar = false
        trainer.train(n_epochs=n_epochs,lr=lr);

        full = trainer.create_posterior(trainer.model, countmatrix);
        latent, batch_indices, labels = full.sequential().get_latent(give_mean=true);

        elbo_per_cell = full.elbo_sample();
        elbo = full.elbo();
    else

        if use_autotune == true
            countmatrix = scvi_dataset.CsvDataset(data_path, save_path = "", new_n_genes = n_genes);

            if subsampling == true
                Random.seed!(111);
                subsample_scVI_cells!(countmatrix, n_cells);
            end

            best_trainer, trials = scvi_inference_autotune.auto_tune_scvi_model(
                    gene_dataset=countmatrix,
                    space=search_space,
                    parallel=false,
                    pickle_result=false,
                    exp_key=string("autotune_out",n_cells),
                    max_evals=30,
            );

            # Set hyperparameters
            use_batches = false;
            use_cuda = true;

            # Train the model and output model likelihood every epoch
            vae = scvi_models.VAE(countmatrix.nb_genes,
                n_batch=countmatrix.n_batches * use_batches,
                n_latent=trials.best_trial["result"]["space"]["model_tunable_kwargs"]["n_latent"],
                n_hidden=trials.best_trial["result"]["space"]["model_tunable_kwargs"]["n_hidden"],
                n_layers=trials.best_trial["result"]["space"]["model_tunable_kwargs"]["n_layers"],
                reconstruction_loss=trials.best_trial["result"]["space"]["model_tunable_kwargs"]["reconstruction_loss"],
                dispersion = "gene-cell",
                latent_distribution = "ln",
                log_variational = false
            );

            trainer = scvi_inference.UnsupervisedTrainer(vae,
                countmatrix,
                use_cuda=use_cuda,
                frequency=1,
                batch_size=batch_size,
                n_epochs_kl_warmup= 30
            );

            full = trainer.create_posterior(trainer.model, countmatrix);
            latent, batch_indices, labels = full.sequential().get_latent(give_mean=true);

            elbo_per_cell_null = full.elbo_sample();
            elbo_null = full.elbo();

            full = best_trainer.create_posterior(best_trainer.model, countmatrix);
            latent, batch_indices, labels = full.sequential().get_latent(give_mean=true);

            elbo_per_cell = full.elbo_sample();
            elbo = full.elbo();
        else
            countmatrix = scvi_dataset.CsvDataset(data_path, save_path = "", new_n_genes = n_genes);

            if subsampling == true
                Random.seed!(111);
                subsample_scVI_cells!(countmatrix, n_cells);
            end

            # Set hyperparameters
            use_batches = false;
            use_cuda = true;

            # Set seeds
            random.seed(111);
            torch.manual_seed(111);

            # Train the model and output model likelihood every epoch
            vae = scvi_models.VAE(countmatrix.nb_genes,
                n_batch=countmatrix.n_batches * use_batches,
                n_latent=n_latent,
                n_hidden=n_hidden,
                n_layers=n_layers,
                dispersion = "gene-cell",
                latent_distribution = "ln",
                log_variational = false
            );

            trainer = scvi_inference.UnsupervisedTrainer(vae,
                countmatrix,
                train_size=0.7,
                use_cuda=use_cuda,
                frequency=1,
                batch_size=batch_size,
                n_epochs_kl_warmup= 30
            );

            full = trainer.create_posterior(trainer.model, countmatrix);
            latent, batch_indices, labels = full.sequential().get_latent(give_mean=true);

            elbo_per_cell_null = full.elbo_sample();
            elbo_null = full.elbo();

            trainer.show_progbar = false
            trainer.train(n_epochs=n_epochs,lr=lr);

            full = trainer.create_posterior(trainer.model, countmatrix);
            latent, batch_indices, labels = full.sequential().get_latent(give_mean=true);

            elbo_per_cell = full.elbo_sample();
            elbo = full.elbo();

            best = [];
            best_trainer = trainer;
        end;
    end;

    if use_autotune == true
        (elbo_per_cell=elbo_per_cell, elbo=elbo, elbo_per_cell_null=elbo_per_cell_null, elbo_null=elbo_null , trainer=best_trainer, model=best_trainer.model, countmatrix=countmatrix, autotune_out=trials.best_trial);
    else
        (elbo_per_cell=elbo_per_cell, elbo=elbo, elbo_per_cell_null=elbo_per_cell_null, elbo_null=elbo_null , trainer=best_trainer, countmatrix=countmatrix);
    end;
end

"""
    null_model_scVI(n_epochs, lr, n_genes, n_cells)
Trains a scVI model and extracts the contribution to the evidence
lower bound for each cell by using the jackstraw approach. A prespecified
number of cells `s` is extracted and the columns are sampled with replacement,
such that for these cells a null distribution is generated.
`n_epochs` number of epochs for training scVI when reference=false.
`n_genes` number of genes to read from dataset.
`n_cells` number of cells to subsample from dataset.
`s` number of samples to use for jackstraw sampling.
`B` number of jackstraw iterations.
`reference` whether to initialize scVI randomly or train the model.
"""
function null_model_scVI(;
    data_path::String,
    n_epochs::Int64,
    lr::Float64,
    n_genes::Int64,
    n_cells::Int64,
    n_latent::Int64,
    n_hidden::Int64,
    n_layers::Int64,
    batch_size::Int64,
    s::Int64,
    B::Int64,
    seed::Int64,
    pretrained_model::String="",
    subsampling::Bool=false)

    #using StatsBase
    #using PyCall
    #using Random
    #using Distributed

    # Load python packages
    #random = pyimport("random")
    #os = pyimport("os")
    #np = pyimport("numpy")
    #pd = pyimport("pandas")
    #scvi = pyimport("scvi")
    #scvi_dataset = pyimport("scvi.dataset")
    #scvi_models = pyimport("scvi.models")
    #scvi_inference = pyimport("scvi.inference")
    #torch = pyimport("torch")

    # TODO: remove random seed and use repetitions instead.
    #Random.seed!(seed);
    elbo_per_cell_append = Array{Float64,1}();
    arr_sample_append = Array{Int64,1}();
    elbo_per_cell_append_null = Array{Float64,1}();
    for i in 1:B

        if isempty(pretrained_model) == false

            # Read data
            countmatrix = scvi_dataset.CsvDataset(data_path, save_path = "", new_n_genes = n_genes);

            if subsampling == true
                Random.seed!(111);
                subsample_scVI_cells!(countmatrix, n_cells);
            end

            # Make copy of original countmatrix
            countmatrix_tmp = deepcopy(countmatrix.X);

            # Load pre-trained model
            trainer = torch.load(pretrained_model);

            # Subsample s jackstraw sample
            s_cells = s;
            arr = collect(1:size(countmatrix.X,1));
            arr_sample = StatsBase.sample(arr, s_cells; replace=false);

            # Sample variables of s jackstraw observations with replacement
            jackstraw_variables = Array{Float64,2}(undef, size(countmatrix.X[arr_sample,:]));
            for y = 1:size(arr_sample,1)
                jackstraw_variables[y,:] .= StatsBase.sample(countmatrix.X[arr_sample[y],:], size(countmatrix.X,2), replace=true);
            end
            countmatrix_tmp[arr_sample,:] .= jackstraw_variables;

            # Override original countmatrix
            countmatrix.X = countmatrix_tmp;

            # Override countmatrix in pre-trained model
            trainer.gene_dataset = countmatrix;

            full = trainer.create_posterior(trainer.model, countmatrix);

            elbo_per_cell_null = full.elbo_sample();

            elbo_per_cell_tmp_null = elbo_per_cell_null[arr_sample];
            append!(elbo_per_cell_append_null, elbo_per_cell_tmp_null);

            # Training based on pre-trained model
            random.seed(111);
            torch.manual_seed(111);
            trainer.show_progbar = false
            trainer.train(n_epochs=n_epochs,lr=lr);

            full = trainer.create_posterior(trainer.model, countmatrix);
            global posterior_object = full;

            elbo_per_cell = full.elbo_sample();
            elbo_per_cell_tmp = elbo_per_cell[arr_sample];
            append!(elbo_per_cell_append, elbo_per_cell_tmp);

            append!(arr_sample_append, arr_sample);

            vae = nothing;
            trainer = nothing;
            full = nothing;
        else

            # Read data
            countmatrix = scvi_dataset.CsvDataset(data_path, save_path = "", new_n_genes = n_genes);

            if subsampling == true
                Random.seed!(111);
                subsample_scVI_cells!(countmatrix, n_cells);
            end

            # Make copy of original countmatrix
            countmatrix_tmp = deepcopy(countmatrix.X);

            # Subsample s jackstraw sample
            s_cells = s;
            arr = collect(1:size(countmatrix.X,1));
            arr_sample = StatsBase.sample(arr, s_cells; replace=false);

            # Sample variables of s jackstraw observations with replacement
            jackstraw_variables = Array{Float64,2}(undef, size(countmatrix.X[arr_sample,:]));
            for y = 1:size(arr_sample,1)
                jackstraw_variables[y,:] .= StatsBase.sample(countmatrix.X[arr_sample[y],:], size(countmatrix.X,2), replace=true);
            end
            countmatrix_tmp[arr_sample,:] .= jackstraw_variables;

            # Override original countmatrix
            countmatrix.X = countmatrix_tmp;

            # Set hyperparameters
            use_batches = false;
            use_cuda = true;

            # Set seeds
            random.seed(111);
            torch.manual_seed(111);

            # Train the model and output model likelihood every epoch
            vae = scvi_models.VAE(countmatrix.nb_genes,
                n_batch=countmatrix.n_batches * use_batches,
                n_latent=n_latent,
                n_hidden=n_hidden,
                n_layers=n_layers,
                dispersion = "gene-cell",
                latent_distribution = "ln",
                log_variational = false
            );

            trainer = scvi_inference.UnsupervisedTrainer(vae,
                countmatrix,
                train_size=0.7,
                use_cuda=use_cuda,
                frequency=1,
                batch_size=batch_size,
                n_epochs_kl_warmup= 30
            );

            full = trainer.create_posterior(trainer.model, countmatrix);

            elbo_per_cell_null = full.elbo_sample();

            elbo_per_cell_tmp_null = elbo_per_cell_null[arr_sample];
            append!(elbo_per_cell_append_null, elbo_per_cell_tmp_null);

            trainer.show_progbar = false
            trainer.train(n_epochs=n_epochs,lr=lr);

            full = trainer.create_posterior(trainer.model, countmatrix);
            global posterior_object = full;

            elbo_per_cell = full.elbo_sample();
            elbo_per_cell_tmp = elbo_per_cell[arr_sample];
            append!(elbo_per_cell_append, elbo_per_cell_tmp);

            append!(arr_sample_append, arr_sample);

            vae = nothing;
            trainer = nothing;
            full = nothing;
        end;
    end;
    (elbo_per_cell_append=elbo_per_cell_append, jackstraw_indices=arr_sample_append, elbo_per_cell_append_null=elbo_per_cell_append_null, posterior_object = posterior_object);
end

function subsample_scVI_cells!(data ,n_cells::Int64)

    # Random cell indices
    arrayy = collect(1:size(data.X,1));
    array_sample = StatsBase.sample(arrayy, n_cells; replace=false);
    array_sample .= array_sample .- 1;

    # Override original countmatrix
    data.update_cells(array_sample);
end

function subsample_scVI_cells_idxs(data ,n_cells::Int64)

    # Random cell indices
    arrayy = collect(1:size(data.X,1));
    array_sample = StatsBase.sample(arrayy, n_cells; replace=false);
    array_sample .= array_sample .- 1;

    # Override original countmatrix
    data.update_cells(array_sample);

    array_sample;
end

"""
    jackstraw_scVI(n_epochs, lr, n_genes, n_cells)
`n_epochs` number of epochs for training scVI when reference=false.
`n_genes` number of genes to read from dataset.
`n_cells` number of cells to subsample from dataset.
`s` number of samples to use for jackstraw sampling.
`B` number of jackstraw iterations.
`reference` whether to initialize scVI randomly or train the model.
"""
function jackstraw_scVI(;
    data_path::String,
    n_epochs::Int64,
    lr::Float64,
    n_genes::Int64,
    n_cells::Int64,
    n_latent::Int64,
    n_hidden::Int64,
    n_layers::Int64,
    batch_size::Int64,
    s::Int64,
    B::Int64,
    seed::Array{Int64,1},
    pretrained_model::String="",
    subsampling::Bool=false,
    use_autotune::Bool=false,
    search_space)

    # Trained models
    alternative_model_scVI_out = alternative_model_scVI(data_path=data_path,
                                                        n_epochs=n_epochs,
                                                        lr=lr,
                                                        n_genes=n_genes,
                                                        n_cells=n_cells,
                                                        n_latent = n_latent,
                                                        n_hidden = n_hidden,
                                                        n_layers = n_layers,
                                                        batch_size = batch_size,
                                                        pretrained_model=pretrained_model,
                                                        subsampling=subsampling,
                                                        use_autotune=use_autotune,
                                                        search_space=search_space
    );

    vae = nothing;
    trainer = nothing;
    full = nothing;

    if isempty(pretrained_model) == false
        null_model_scVI_out = vcat(pmap(params -> null_model_scVI(data_path=data_path,
                                                                n_epochs=params[1],
                                                                lr=params[2],
                                                                n_genes=params[3],
                                                                n_cells=params[4],
                                                                n_latent = params[5],
                                                                n_hidden = params[6],
                                                                n_layers = params[7],
                                                                batch_size = params[8],
                                                                s=params[9],
                                                                B=params[10],
                                                                seed=params[11],
                                                                pretrained_model=pretrained_model,
                                                                subsampling=subsampling),
                Iterators.product(n_epochs, lr, n_genes, n_cells, n_latent, n_hidden, n_layers, batch_size, s, B, seed))...
        );
    else
        if use_autotune == true
            null_model_scVI_out = vcat(pmap(params -> null_model_scVI(data_path=data_path,
                                                                    n_epochs=alternative_model_scVI_out.autotune_out["result"]["best_epoch"],
                                                                    lr=alternative_model_scVI_out.autotune_out["result"]["space"]["train_func_tunable_kwargs"]["lr"],
                                                                    n_genes=params[1],
                                                                    n_cells=params[2],
                                                                    n_latent=alternative_model_scVI_out.autotune_out["result"]["space"]["model_tunable_kwargs"]["n_latent"],
                                                                    n_hidden=alternative_model_scVI_out.autotune_out["result"]["space"]["model_tunable_kwargs"]["n_hidden"],
                                                                    n_layers=alternative_model_scVI_out.autotune_out["result"]["space"]["model_tunable_kwargs"]["n_layers"],
                                                                    batch_size = params[3],
                                                                    s=params[4],
                                                                    B=params[5],
                                                                    seed=params[6],
                                                                    pretrained_model=pretrained_model,
                                                                    subsampling=subsampling),
                Iterators.product(n_genes, n_cells, batch_size, s, B, seed))...
            );
        else
            null_model_scVI_out = vcat(pmap(params -> null_model_scVI(data_path=data_path,
                                                                    n_epochs=params[1],
                                                                    lr=params[2],
                                                                    n_genes=params[3],
                                                                    n_cells=params[4],
                                                                    n_latent = params[5],
                                                                    n_hidden = params[6],
                                                                    n_layers = params[7],
                                                                    batch_size = params[8],
                                                                    s=params[9],
                                                                    B=params[10],
                                                                    seed=params[11],
                                                                    pretrained_model=pretrained_model,
                                                                    subsampling=subsampling),
                    Iterators.product(n_epochs, lr, n_genes, n_cells, n_latent, n_hidden, n_layers, batch_size, s, B, seed))...
            );
        end;
    end;

    interrupt();
    rmprocs(workers());

    (alternative_model_scVI_out=alternative_model_scVI_out,
    null_model_scVI_out=null_model_scVI_out);
end


text:
#=
   SEND+MORE=MONEY any base in Julia + ConstraintSolver.jl

   Number of solutions for base 2:30
   [0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 120, 136, 153, 171, 190, 210, 231]

   Model created by Hakan Kjellerstrand, hakank@gmail.com
   See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function send_more_money_any_base(base=10,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,
                                            "all_solutions"=> all_solutions,
                                            # "all_optimal_solutions"=>all_solutions,
                                            "logging"=>[],

                                            "traverse_strategy"=>:BFS,
                                            # "traverse_strategy"=>:DFS,
                                            # "traverse_strategy"=>:DBFS,

                                            # "branch_split"=>:Smallest,
                                            # "branch_split"=>:Biggest,
                                            "branch_split"=>:InHalf,

                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                            # "branch_strategy" => :IMPS, # default
                                            "branch_strategy" => :ABS, # Activity Based Search
                                            # "activity.decay" => 0.999, # default 0.999
                                            # "activity.max_probes" => 100, # default, 10
                                            # "activity.max_confidence_deviation" => 20, # default 20

                                            # "simplify"=>false,
                                            # "simplify"=>true, # default

                                            "time_limit"=>6,

                                            # "backtrack" => false, # default true
                                            # "backtrack_sorting" => false, # default true

                                            # "lp_optimizer" => cbc_optimizer,
                                            # "lp_optimizer" => glpk_optimizer,
                                            # "lp_optimizer" => ipopt_optimizer,
                                            ))

    println("base: $base")
    @variable(model, 0 <= x[1:8] <= base-1, Int)
    s,e,n,d,m,o,r,y = x

    @constraint(model, x in CS.AllDifferentSet())
    @constraint(model, s != 0)
    @constraint(model, m != 0)

    @constraint(model,
                                 (base^3)*s + (base^2)*e + base*n + d
                               + (base^3)*m + (base^2)*o + base*r + e
                 == (base^4)*m + (base^3)*o + (base^2)*n + base*e + y
        )

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")
                println()
            end
        end
        println("base:$base num_sols:$num_sols\n")
    else
        println("status:$status")
    end

    return status,num_sols
end

num_sols_a = []
for base in 2:30
    local status
    @time status,num_sols = send_more_money_any_base(base,true,true)
    if status == MOI.TIME_LIMIT
        break
    else
        push!(num_sols_a,num_sols)
    end

end

println("Num solutions: $num_sols_a")

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module MatchExt
import MLStyle: Feature, (⇒), Fun
import MLStyle.Match: pattern_match, PatternDef

using MLStyle.Private
using MLStyle.Err
using MLStyle

export (..), enum_next
function enum_next(x :: Number)
    x + 1
end

function (..)(predicate, start :: T, _end :: T) where T
    while start < _end
        if predicate(start)
            return true
        end
        start = enum_next(start)
    end
    false
end

pattern_match(num :: Number, guard, tag, mod :: Module) =
    if guard === nothing
        :($tag === $num)
    else
        :($tag === $num && $guard)
    end

pattern_match(str :: AbstractString, guard, tag, mod :: Module) =
    if guard === nothing
        :($tag == $str)
    else
        :($tag == $str && $guard)
    end

pattern_match(quoted :: QuoteNode, guard, tag, mod :: Module) =
    if guard === nothing
        :($tag == $quoted)
    else
        :($tag == $quoted && $guard)
    end

pattern_match(sym :: Symbol, guard, tag, mod :: Module) =

        if sym === :_
            if nothing === guard
                quote true end
            else
                guard
            end
        else
            let ret =
                quote
                    $sym = $tag
                    true
                end
                if guard === nothing
                    ret
                else
                    :($ret && $guard)
                end
            end
        end
# """
# @match :(x + 1) begin
#     :(x + 1) => 1
# end
# """
PatternDef.Meta(expr :: Expr -> expr.head === :quote) do expr, guard, tag, mod
    if guard === nothing
        :($tag == $expr)
    else
        :($tag == $expr && $guard)
    end
end

# """
# like ^ in Erlang/Elixir
# """
PatternDef.Meta(expr :: Expr -> expr.head == :&)     do expr, guard, tag, mod
    value = expr.args[1]
    if guard === nothing
        :($tag === $value)
    else
        :($tag === $value && $guard)
    end
end

# """
# tuple pattern
# """
PatternDef.Meta(expr :: Expr -> expr.head == :tuple) do expr, guard, tag, mod
    args = expr.args
    len = length(args)

    if len !== 0
        check_len = :($length($tag) === $len)
        matching = map(zip(1:len, args)) do (idx, arg)
            pattern_match(arg, nothing, :($tag[$idx]), mod)
        end |>
        function (last)
            reduce((a, b) -> Expr(:&&, a, b), last)
        end |>
        function(last)
            :($check_len && $last)
        end
    else
        :($isempty($tag))
    end |>
    function (last)
        if guard !== nothing
            :($last && $guard)
        else
            last
        end
    end
end

# """
# @match expr {
#     :: Ty => # do something if isa(expr, Ty)
#     x :: Ty => # do something if isa(expr, Ty) and perform capturing.
# }
# """
PatternDef.Meta((expr :: Expr -> expr.head == :(::))) do expr, guard, tag, mod

    args = expr.args
    len = length(args)

    annotation_processing(arg) =
        if Feature.is_activated(:TypeLevel, mod)
            pattern_match(arg, nothing, :($typeof($tag)), mod)
        else
            :($isa($tag, $(arg)))
        end

    check_ty =
        if len === 1
            annotation_processing(args[1])
        else
            @assert len === 2
            pat, ty = args
            let guard = pattern_match(pat, nothing, tag, mod)
                annotation_processing(ty) |>
                function (last)
                    if guard === nothing
                        last
                    else
                        :($last && $guard)
                    end
                end
            end
        end
    if guard === nothing
        check_ty
    else
        :($check_ty && $guard)
    end
end

# """
#  [a, b, c, d..., e, f, g] => ...
# """
PatternDef.Meta(expr::Expr -> expr.head == :vect) do expr, guard, tag, mod
    args = expr.args
    if length(args) == 0
        return :($isempty($tag)) |>
                last ->
                if guard === nothing
                    last
                else
                    :($last && $guard)
                end
    end

    atleast_element_count = 0
    unpack_begin = nothing
    unpack_end = 0
    unpack = []
    foreach(args) do arg
        if isa(arg, Expr) && arg.head === :...
            if unpack_begin === nothing
                unpack_begin = atleast_element_count + 1
             else
                SyntaxError("Vector unpack can only perform sequential unpack at most once.") |> throw
            end
            push!(unpack, arg.args[1])
        else
            atleast_element_count = atleast_element_count + 1
            if unpack_begin !== nothing
                push!(unpack, pattern_match(arg, nothing, :($tag[end-$unpack_end]), mod))
                unpack_end = unpack_end + 1
            else
                push!(unpack, pattern_match(arg, nothing, :($tag[$atleast_element_count]), mod))
            end
        end
    end

    if unpack_begin !== nothing
        arg = unpack[unpack_begin]
        unpack[unpack_begin] = pattern_match(arg, nothing, :($view($tag, $unpack_begin:($length($tag) - $unpack_end))), mod)
        reduce(unpack, init=:($length($tag) >= $atleast_element_count)) do last, each
            :($last && $each)
        end
    else
        reduce(unpack, init=:($length($tag) == $atleast_element_count)) do last, each
            :($last && $each)
        end
    end
end
# """
# Feature.@activate TypeLevel
# @match Int ⇒ Int begin
#     ::(T ⇒ G) => (T, G)
# end # => (Int, Int)
# """
PatternDef.Meta(expr :: Expr -> expr.head == :curly) do expr, guard, tag, mod

    if !Feature.is_activated(:TypeLevel, mod)
        :($tag isa $expr)
    else
        @match expr.args begin
            [head, tail...] =>
            begin
                len = length(tail)

                pat1 = pattern_match(head, nothing, :($tag.name.wrapper), mod)

                patn = map(enumerate(tail)) do (idx, each)
                    pattern_match(each, nothing, :($tag.parameters[$idx]), mod)
                end |> last -> reduce(ast_and, last, init = :($length($tag.parameters) == $len))

                :($pat1 && $patn)
            end
            _ => SyntaxError("$(string(expr))") |> throw
        end
    end
end

# """
# @match 1 {
#     1..10   in x => x
#     11..20  in x => x * 10
# }
# """
PatternDef.App(..) do args, guard, tag, mod

    start, _end = args
    quote
        let value = $tag
        ($..)($start, $_end) do it
            it === value
        end
        end
    end
end

PatternDef.App(⇒) do args, guard, tag, mod
    from, to = args
    pat1 = pattern_match(from, nothing, :($tag.parameters[1]), mod)
    pat2 = pattern_match(to,   nothing, :($tag.parameters[2]), mod)
    :($tag <: $Fun && $pat1 && $pat2)
end

end


text:
include("basics.jl")
include("productbench.jl")