text:
using Clusters
using DataFrames
using Base.Test
using CovarianceMatrices
using CSV
const sim = 300

## Two test.

#=
 1. - Make sure that the estimates obtained from estimatemodel
      coincides with those of STATA
=#

df = CSV.read("test.txt")
ng = vec(readcsv("states.txt", header = true)[1])
ng = floor.(Int64, ng./5000)

# G::Int64
# ng::Vector{Int64}
# ## Regressors std. dev.
# σ_z::Float64
# σ_ξ::Float64
# ## Errors std. dev
# σ_ϵ::Float64
# σ_α::Float64
# ## Third/Fourth design parms
# p::Float64
# γ::Float64
# δ::Float64
# σ_η::Float64
# μ_q::Float64
# ## Null value
# β₀::Float64
# design::Int

opt = LinearRegressionClusterOpt(length(ng), ng,
                                  1,  ## σ_z
                                1/9,  ## σ_ξ
                                  1,  ## σ_ϵ
                                1/9,  ## σ_α
                                0.0,  ## p
                                0.0,  ## γ
                                0.0,  ## δ
                                0.0,  ## σ_η
                                0.0,  ## μ_q
                                0.0,  ## β₀
                                  1)
m = initialize(LinearRegressionCluster, opt)
simulate!(m)
## Inject data from test.txt
copy!(m.y, df[:y])
copy!(m.X, df[:x])
out = estimatemodel(m)
## Make sure everything is fine
@testset "Set 1" begin
      @test m.cl == df[:cl]
      @test m.wts == df[:wts]
      @test all(m.y .== df[:y])
      @test all(m.X .== df[:x])
      @testset "Coefficients and Variances" begin
            @test out[1] ≈ -.0421262 atol = 0.0001  ## Unweighted OLS
            @test out[2] ≈ 0.0275297 atol = 0.0001  ## iid std
            @test out[3] ≈ 0.0276643 atol = 0.0001  ## HC std
            @test out[4] ≈ 0.0296611 atol = 0.0001  ## CRHC1 std
            @test out[5] ≈ 0.0296611 atol = 0.01    ## CRHC2 std
            @test out[6] ≈ 0.0296611 atol = 0.01    ## CRHC3 std

            @test out[17] ≈  -.1051925 atol = 0.0001 ## Weighted OLS
            @test out[18] ≈  0.0269856 atol = 0.0001 ## iid std
            @test out[19] ≈  0.0427102 atol = 0.0001 ## HC std
            @test out[20] ≈  0.0487519 atol = 0.0001 ## CRHC1 std
      end
      @testset "Check correctness of parameters" begin
            σ_ϵ = m.opt.σ_z
            σ_α = m.opt.σ_ξ
            σ_z = m.opt.σ_z
            σ_ξ = m.opt.σ_ξ
            γ   = m.opt.γ
            p   = m.opt.p
            ## The parameters are std. dev
            @test σ_ϵ == 1
            @test σ_α == 1/9
            @test σ_z == 1
            @test σ_ξ == 1/9
            ## Check ICC of regressors and errors
            @test icc_x(m) == σ_ξ^2/(σ_z^2+σ_ξ^2)
            @test icc_u(m) == σ_α^2/(σ_ϵ^2+σ_α^2)
      end
end

#=
# 2. - Run a Small Montecarlo
=#
ng = repeat([30], inner = 30)
opt = opt = LinearRegressionClusterOpt(length(ng),
                                       ng,
                                       1.0, ## σ_z
                                       1.0, ## σ_ξ
                                       1.0, ## σ_ϵ
                                       1.0, ## σ_α
                                       0.0, ## p
                                       0.0, ## γ
                                       0.0, ## δ
                                       0.0, ## σ_η
                                       0.0, ## μ_q
                                       0.0, ## β₀
                                       1) ## Design

m = initialize(LinearRegressionCluster, opt)
srand(1234512345)
out = montecarlo(LinearRegressionCluster, opt)

@testset "Set 2" begin
      @test maximum(out[:theta_u] -  out[:theta_w]) ≈ 0 atol = 1e-10
      ## Confidence interval
      ci = ConfInterval(out[:theta_u], out[:V3_u], 1.96)
      @test average_length(ci) <= 0.4
      ## Construct WB confidence interval
      ciwb = ConfInterval(out[:theta_w], out[:V3_w], out[:qw0_025], out[:qw0_975])
      @test coverage(ciwb, 0) ≈ 0.93 atol = 0.01
      ciwb = ConfInterval(out[:theta_u], out[:V3_u], out[:qu0_025], out[:qu0_975])
      @test coverage(ciwb, 0) ≈ 0.93 atol = 0.01
      ciwb = ConfInterval(out[:theta_w], out[:V3_w], out[:qw0_050], out[:qw0_950])
      @test coverage(ciwb, 0) ≈ 0.90 atol = 0.04
      ciwb = ConfInterval(out[:theta_u], out[:V3_u], out[:qu0_050], out[:qu0_950])
      @test coverage(ciwb, 0) ≈ 0.90 atol = 0.04
end

function run_montecarlo(d, sim)
      dd = Array{Float64}(sim)
      od = similar(dd)
      N = length(d.opt.ng)
      ng = maximum(d.opt.ng)
      M = Array{Float64}(ng,ng);

      for i in 1:sim
            simulate!(d)
            fill!(M, 0.0)
            for j in d.bstarts
                  E = d.y[j]
                  M = M + E*E'./N
            end
            dd[i] = mean(diag(M))
            od[i] = M[1,2]
      end
      (dd, od)
end

ci_lower_lim(x, z = 2.5) = mean(x) - z*std(x)/sqrt(length(x))
ci_upper_lim(x, z = 2.5) = mean(x) + z*std(x)/sqrt(length(x))

theoretical_Λ_ii(d, ::Type{Val{1}}) = d.opt.σ_ϵ^2+d.opt.σ_α^2
theoretical_Λ_ik(d, ::Type{Val{1}}) = d.opt.σ_α^2

theoretical_Λ_ii(d, ::Type{Val{2}}) = theoretical_Λ_ii(d, Val{1})
theoretical_Λ_ik(d, ::Type{Val{2}}) = theoretical_Λ_ik(d, Val{1})

function theoretical_Λ_ik(d, ::Type{Val{3}})
      (σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, μ_q) = Clusters.getparms(d)
      σ_α^2+γ^2*(1-p)*p*σ_ξ^2
end

function theoretical_Λ_ii(d, ::Type{Val{3}})
      (σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, μ_q) = Clusters.getparms(d)
      σ_ϵ^2 + σ_α^2+γ^2*(1-p)*p*(1+σ_ξ^2)
end

function theoretical_Λ_ik(d, ::Type{Val{4}})
      (σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, μ_q) = Clusters.getparms(d)
      a = σ_η^2 + δ*σ_ξ^2*(2*δ + 3*μ_q) + δ^2*σ_z^2
      γ^2*σ_ξ^2*a
end

function theoretical_Λ_ii(d, ::Type{Val{4}})
      (σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, μ_q) = Clusters.getparms(d)
      a = 2*δ^2*σ_ξ^4
      b = σ_ξ^2*(σ_η^2 + μ_q^2 + δ^2*σ_z^2)
      c = σ_z^2*(σ_η^2 + μ_q^2)
      1 + γ^2*(a+b+c)
end

@testset "Design 1" begin
      ng = repeat([2], outer=2000)
      opt = LinearRegressionClusterOpt(length(ng), ng,
                                       1,    ## σ_z
                                       1/9,  ## σ_ξ
                                       1,    ## σ_ϵ
                                       1/9,  ## σ_α
                                       0.0,  ## p
                                       0.0,  ## γ
                                       0.0,  ## δ
                                       0.0,  ## σ_η
                                       0.0,  ## μ_q
                                       0.0,  ## β₀
                                        1)   ## Desing

      d = initialize(LinearRegressionCluster, opt)
      simulate!(d)
      srand(1234554321)
      Λ_ii, Λ_ik = run_montecarlo(d, sim)
      @test ci_lower_lim(Λ_ii) <= theoretical_Λ_ii(d, Val{1})
      @test ci_upper_lim(Λ_ii) >= theoretical_Λ_ii(d, Val{1})

      @test ci_lower_lim(Λ_ik) <= theoretical_Λ_ik(d, Val{1})
      @test ci_upper_lim(Λ_ik) >= theoretical_Λ_ik(d, Val{1})
end

@testset "Design 2" begin
      ng = repeat([2], outer=2000)

      opt = LinearRegressionClusterOpt(length(ng), ng,
                                       1,    ## σ_z
                                       1/9,  ## σ_ξ
                                       1,    ## σ_ϵ
                                       1/9,  ## σ_α
                                       0.0,  ## p
                                       0.0,  ## γ
                                       0.0,  ## δ
                                       0.0,  ## σ_η
                                       0.0,  ## μ_q
                                       0.0,  ## β₀
                                        2)   ## Desing

      d = initialize(LinearRegressionCluster, opt)
      simulate!(d)
      srand(1234554321)
      Λ_ii, Λ_ik = run_montecarlo(d, sim)
      @test ci_lower_lim(Λ_ii) <= theoretical_Λ_ii(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ii) >= theoretical_Λ_ii(d, Val{d.opt.design})

      @test ci_lower_lim(Λ_ik) <= theoretical_Λ_ik(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ik) >= theoretical_Λ_ik(d, Val{d.opt.design})

      opt = LinearRegressionClusterOpt(length(ng), ng,
                                       1,    ## σ_z
                                       1/9,  ## σ_ξ
                                       1,    ## σ_ϵ
                                       1/2,  ## σ_α
                                       0.0,  ## p
                                       0.0,  ## γ
                                       0.0,  ## δ
                                       0.0,  ## σ_η
                                       0.0,  ## μ_q
                                       0.0,  ## β₀
                                        2)   ## Desing

      srand(1234554321)
      Λ_ii, Λ_ik = run_montecarlo(d, sim)
      @test ci_lower_lim(Λ_ii) <= theoretical_Λ_ii(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ii) >= theoretical_Λ_ii(d, Val{d.opt.design})

      @test ci_lower_lim(Λ_ik) <= theoretical_Λ_ik(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ik) >= theoretical_Λ_ik(d, Val{d.opt.design})
end

@testset "Design 3" begin
      ng = repeat([2], outer=20000)
      ## ρᵤ = 0.2  ρₓ = 0.2
      opt = LinearRegressionClusterOpt(length(ng), ng,
                                       1,    ## σ_z
                                       0.5,  ## σ_ξ
                                       1,    ## σ_ϵ
                                       0.5,  ## σ_α
                                       0.5,  ## p
                                       0.1,  ## γ
                                       0.0,  ## δ
                                       0.0,  ## σ_η
                                       0.0,  ## μ_q
                                       0.0,  ## β₀
                                        3)   ## Desing

      d = initialize(LinearRegressionCluster, opt)
      simulate!(d)
      srand(1234554321)
      Λ_ii, Λ_ik = run_montecarlo(d, sim)
      @test ci_lower_lim(Λ_ii) <= theoretical_Λ_ii(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ii) >= theoretical_Λ_ii(d, Val{d.opt.design})
      @test ci_lower_lim(Λ_ik) <= theoretical_Λ_ik(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ik) >= theoretical_Λ_ik(d, Val{d.opt.design})
end

@testset "Design 4" begin
      ## This test that as n_g -> ∞
      ng = repeat([150], outer=15)

      opt = LinearRegressionClusterOpt(length(ng), ng,
                                       1, ## σ_z
                                     0.5, ## σ_ξ
                                       0, ## σ_ϵ
                                       0, ## σ_α
                                     0.0, ## p
                                    -1.0, ## γ
                       .9354143466934853, ## δ
                       .5361902647381804, ## σ_η
                                     -0.5, ## μ_q
                                     0.0, ## β₀
                                       4) ## Design

      d = initialize(LinearRegressionCluster, opt)
      simulate!(d)
      srand(1234554321)
      Λ_ii, Λ_ik = run_montecarlo(d, 300)
      @test ci_lower_lim(Λ_ii) <= theoretical_Λ_ii(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ii) >= theoretical_Λ_ii(d, Val{d.opt.design})
      @test ci_lower_lim(Λ_ik) <= theoretical_Λ_ik(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ik) >= theoretical_Λ_ik(d, Val{d.opt.design})

      ng = repeat([150], outer=10)
      opt = LinearRegressionClusterOpt(length(ng), ng,
                                        1, ## σ_z
                                      0.5, ## σ_ξ
                                        0, ## σ_ϵ
                                        0, ## σ_α
                                      0.0, ## p
                                     -1.0, ## γ
                        .9354143466934853, ## δ
                        .5361902647381804, ## σ_η
                                     -0.5, ## μ_q
                                      0.0, ## β₀
                                        4) ## Design

      d = initialize(LinearRegressionCluster, opt)
      simulate!(d)
      srand(1234554321)
      Λ_ii, Λ_ik = run_montecarlo(d, sim)
      @test ci_lower_lim(Λ_ii) <= theoretical_Λ_ii(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ii) >= theoretical_Λ_ii(d, Val{d.opt.design})
      @test ci_lower_lim(Λ_ik) <= theoretical_Λ_ik(d, Val{d.opt.design})
      @test ci_upper_lim(Λ_ik) >= theoretical_Λ_ik(d, Val{d.opt.design})
end

parm_design_1 = readcsv("parameters/parameters-design1.csv")
parm_design_2 = readcsv("parameters/parameters-design2.csv")
parm_design_3 = readcsv("parameters/parameters-design3.csv")
parm_design_4 = readcsv("parameters/parameters-design4.csv")

@testset "icc_x and icc_u" begin
      ng = repeat([150], outer=20)
      rhox = repeat([.2, .5, .9, .99], outer = 4)
      rhou = repeat([.2, .5, .9, .99], inner = 4)

      parm_design_1 = mapslices(u->[u], parm_design_1, 2)
      parm_design_2 = mapslices(u->[u], parm_design_2, 2)
      parm_design_3 = mapslices(u->[u], parm_design_3, 2)
      parm_design_4 = mapslices(u->[u], parm_design_4, 2)

      for jopt in [parm_design_1, parm_design_2, parm_design_3, parm_design_4]
            opt = map(u -> LinearRegressionClusterOpt(length(ng), ng, u...), jopt)
            @test all(map(u -> icc_x(u), opt) .≈ rhox)
            @test all(map(u -> icc_u(u), opt) .≈ rhou)
      end

end

parm_design_1 = convert(Array, CSV.read("parameters/parameters-design1-xu.csv", datarow = 1))

@testset "icc_x and icc_u" begin
      ng = repeat([150], outer=20)
      rhoxu = [.2, .5, .9, .99]

      parm_design_1 = mapslices(u->[u], parm_design_1, 2)
      parm_design_2 = mapslices(u->[u], parm_design_2, 2)
      parm_design_3 = mapslices(u->[u], parm_design_3, 2)
      parm_design_4 = mapslices(u->[u], parm_design_4, 2)

      for jopt in [parm_design_1, parm_design_2, parm_design_3, parm_design_4]
            opt = map(u -> LinearRegressionClusterOpt(length(ng), ng, u...), jopt)
            @test all(map(u -> Clusters.icc_xu(u), opt) .≈ rhoxu)
      end
end

using RCall
R"library('ICC')"

parm_design_1 = convert(Array, CSV.read("parameters/parameters-design1-xu.csv", datarow = 1))
parm_design_2 = convert(Array, CSV.read("parameters/parameters-design2-xu.csv", datarow = 1))
parm_design_3 = convert(Array, CSV.read("parameters/parameters-design3-xu.csv", datarow = 1))
parm_design_4 = convert(Array, CSV.read("parameters/parameters-design4-xu.csv", datarow = 1))

parm_design_1 = mapslices(u->[u], parm_design_1, 2)
parm_design_2 = mapslices(u->[u], parm_design_2, 2)
parm_design_3 = mapslices(u->[u], parm_design_3, 2)
parm_design_4 = mapslices(u->[u], parm_design_4, 2)

@testset "icc_x and icc_u" begin
      ng = repeat([20], outer=10000)
      rhoxu = [.2, .5, .9, .99]

      #for jopt in  [parm_design_1, parm_design_2, parm_design_3, parm_design_4]
      jopt = parm_design_4
            opt = map(u -> LinearRegressionClusterOpt(length(ng), ng, u...), jopt)
            m = map(u -> initialize(LinearRegressionCluster, u), opt)
            map(simulate!, m)
            est_rhoxu = map(u -> begin
                  xu = u.X.*(u.y)
                  cl = u.cl
                  est_rhoxu = R"ICCbare(cl, xu, data.frame(cl=$cl, xu=$(xu)))"
                  (est_rhoxu[1], var(u.X), var(u.y))
            end, m)

            @test est_rhoxu ≈ rhoxu atol = 0.1
      end
end


text:
# rotatez.jl

using SPECTrecon: rotl90!, rotr90!, rot180!, rot_f90!, rot_f90_adj!
using Test: @test, @testset, @test_throws

@testset "rotl90!" begin
    T = Float32
    N = 20
    A = rand(T, N, N)
    B = similar(A, N, N)
    rotl90!(B, A)
    @test isequal(B, rotl90(A))
end

@testset "rotr90!" begin
    T = Float32
    N = 20
    A = rand(T, N, N)
    B = rand(T, N, N)
    rotr90!(B, A)
    @test isequal(B, rotr90(A))
end

@testset "rot180!" begin
    T = Float32
    N = 20
    A = rand(T, N, N)
    B = rand(T, N, N)
    rot180!(B, A)
    @test isequal(B, rot180(A))
end

@testset "rot_f90!" begin
    T = Float32
    N = 20
    A = rand(T, N, N)
    B = rand(T, N, N)
    @test_throws String rot_f90!(A, B, -1)
    @test_throws String rot_f90_adj!(A, B, -1)
    @test_throws String rot_f90!(A, B, 4)
    @test_throws String rot_f90_adj!(A, B, 4)
end

ext: jl


text:
ENV_name = "owl"
ENV_root = "output_owl"
println("Environment: OWL")

#memory=4G
#gpu_mem=4000M
function jobengine_head(name::String, dir::String, Ncores::Integer, gpu::Bool; hours::Int64=48, architecture::String="", memory::String="", gpu_mem::String="")
    # gpu = true #Shortcut gpu setting to not make job engine hickup
    return """
    #!/bin/bash
    #\$ -S /bin/bash
    $(Ncores > 1 ? "#\$ -pe openmp_fast $Ncores" : "")
    #\$ -q *
    #\$ -e $dir/terminal-err
    #\$ -o $dir/terminal-out
    $(architecture != "" ? "#\$ -l cm=$(architecture)" : "")
    #\$ -N $(jobname(name))
    #\$ -M bardenn@gwdg.de
    #\$ -m n
    #\$ -l h_rt=$(hours):00:00
    #\$ -wd $dir
    #\$ -V
    #\$ -hold_jid $(jobname(name))
    $(memory != "" ? "#\$ -l h_data=$memory" : "")
    $(gpu ? "#\$ -l gpu=1" : "")
    $(gpu_mem != "" ? "#\$ -l gpu_mem=$(gpu_mem)" : "")
    """
end

"""Launches the script via qsub, in case dir doesnt exist, create it"""
function launch_job(dir::String, Ncores::Integer, gpu::Bool, julia_script::String, successive_jobs::Integer=1; hours::Int64=48, architecture::String="", memory::String="", fresh::Bool=false, gpu_mem::String="")
    githead = check_git_status()
    full_dir = environment_path(dir)
    head = jobengine_head(dir, full_dir, Ncores, gpu; hours=hours, architecture=architecture, memory=memory, gpu_mem=gpu_mem)

    if fresh run(`rm -rf $(full_dir)`) end

    run(`mkdir -p $(full_dir)`)
    open("$(full_dir)/job.sh", "w+") do file
        write(file,"""#git-SHA: $githead
        $head

        $(Ncores > 1 ? "julia -p8" : "julia") << EOF
        $julia_script
        flush(STDOUT) #make sure output is flushed to terminal files
        EOF
        """)
    end

    #Let's queue some jobs, hold_jid takes care of chaining
    for n = 1:successive_jobs
        submit_job("$(full_dir)/job.sh")
    end
end

"""Submits a job to the queue"""
function submit_job(path)
    run(`qsub $path`)
end


text:
@testset "718.maximum-length-of-repeated-subarray.jl" begin
    @test find_length718([1, 2, 3, 2, 1], [3, 2, 1, 4, 7]) == 3
    @test find_length718([1, 2, 3, 2, 1], [3, 2, 4, 1, 7]) == 2
end


text:
using FastMarching
using Base.Test
using LinearAlgebra

function eye(n)
  fill(0,(n,n))+I
end

function testeye(n::Integer,endpoint=[1.,1.],stepsize=0.1)
  speedmap = eye(n) .* 1000 .+ 0.001
  source = [float(size(speedmap,1)), float(size(speedmap,2))]
  distancemap = FastMarching.msfm(speedmap,source,true,true)
  endpoint = [1.,1.]
  shortestline = FastMarching.shortestpath(distancemap,endpoint,source,0.1)
end

function runtests()
  const endpoint = [1.,1.]
  for t in 2:2:10, stepsize in [0.1,0.01]
    println("size: $(lpad(t,2," ")) stepsize: $(rpad(string(stepsize),4,"0"))")
    shortestline = testeye(t,endpoint,stepsize)
    @test sum(shortestline[:,2]-shortestline[:,1]) == 0
    @test sqrt((shortestline[1,1]-1.)^2-(shortestline[1,2]-1.)^2) <= stepsize
  end

  try
    rm("multiplestarts.png")
  end
  include(joinpath(Pkg.dir("FastMarching"),"examples/multiplestart.jl"))
  @test isfile("multiplestart.png")

end

runtests()


text:
using JSON

include("./jupyter.jl")

install_jupyter_nbextension()
try
    install_jupyter_labextension()
catch exc
    @warn(
        "Unable to install Jupyter Lab extension; make sure jupyter is in your "
            * "PATH and rebuild WebIO if you want to use Jupyter Lab.",
        exception=exc,
    )
end


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

The exported names are:

$(EXPORTS)
"""
module Polynomials

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

import Gridap.Fields: evaluate_field!
import Gridap.Fields: field_cache
import Gridap.Fields: evaluate_gradient!
import Gridap.Fields: gradient_cache
import Gridap.Fields: evaluate_hessian!
import Gridap.Fields: hessian_cache
import Gridap.Arrays: kernel_return_type

export MonomialBasis
export QGradMonomialBasis
export QCurlGradMonomialBasis
export change_basis
export get_exponents

export get_order
export get_orders
export get_value_type
export num_terms

include("MonomialBases.jl")

include("QGradMonomialBases.jl")

include("QCurlGradMonomialBases.jl")

include("ChangeBasis.jl")

end # module


text:
include("c3.jl")
include("boundarylayer.jl")
include("stomata.jl")
include("intercellularspace.jl")
include("irradiance.jl")
include("energybalance.jl")

@system GasExchange(
    WeatherStub, SoilStub,
    BoundaryLayer, Stomata, IntercellularSpace, Irradiance, EnergyBalance,
    C3
) begin
    weather ~ ::Weather(override)
    soil ~ ::Soil(override)

    #FIXME: confusion between PFD vs. PPFD
    PPFD: photosynthetic_photon_flux_density ~ track(u"μmol/m^2/s" #= Quanta =#, override)
    LAI: leaf_area_index ~ track(override)

    A_net_total(A_net, LAI): net_photosynthesis_total => A_net * LAI ~ track(u"μmol/m^2/s" #= CO2 =#)
    A_gross_total(A_gross, LAI): gross_photosynthesis_total => A_gross * LAI ~ track(u"μmol/m^2/s" #= CO2 =#)
    E_total(E, LAI): transpiration_total => E * LAI ~ track(u"mmol/m^2/s" #= H2O =#)

    #TODO: nitrogen response not implemented for C3 yet
    #N: nitrogen => 2.0 ~ preserve(parameter)
end


text:
function cem(system::MarketClearingProblem{Z, T},
    solver::JuMP.MOI.OptimizerWithAttributes,
    resultfile::String="") where {Z, T}

    lines = [line.name for line in system.lines] # Lines
    projects = [project.name for project in system.projects] # Projects
    invperiods = 1:length(system.inv_periods)     # Investment periods
    opperiods = 1:T
    rep_hour_weight = system.rep_hour_weight
    capital_cost_multiplier = system.capital_cost_multiplier

    end_of_day = collect(24:24:T)

    tech_type,                                # technology type of project
    project_type,                             # is the project of storage type or generator type
    marginal_energy_cost,                      # $/MW
    marginal_reserveup_cost,                   # $/MW
    marginal_reservedown_cost,                 # $/MW
    fixed_cost,                                # $/unit
    queue_cost,                                # $/unit/year
    expansion_cost,                            # $/unit
    project_discount_rate,                    # discount rate
    min_gen, max_gen,                           # MW/unit
    min_input, max_input,                       # MW/unit
    efficiency_in, efficiency_out,            # percentage
    min_storage, max_storage,                   # MWh/unit
    init_storage,                              # MWh/unit
    availability,                             # percentage
    derating_factor,                           # percentage
    ramp_limits,                               # MW/unit/hour
    max_reserveup, max_reservedown,             # MW/unit/hour
    existing_units,                            # existing units
    units_in_queue,                             # units in queue
    build_lead_time, remaining_buildtime,       # construction lag times
    max_new_options,                            # maximum units
    base_cost_units,                            # base cost units (for capital cost multiplier)
    capex_years,                               # for investment cost annualization
    life_time,                                 # total possible life_time
    remaining_life_time,                       # remaining life_time
    capacity_eligible,                        # project eligible for capacity market
    rec_eligible,                             # project eligible for REC market
    location =                                # project zone location
    make_parameter_vector.(
    Ref(system.projects), :name,
    [:tech_type, :project_type, :marginal_energy_cost, :marginal_reserveup_cost, :marginal_reservedown_cost,
    :fixed_cost, :queue_cost, :expansion_cost, :discount_rate,
    :min_gen, :max_gen, :min_input, :max_input,
    :efficiency_in, :efficiency_out, :min_storage, :max_storage, :init_storage,
    :availability, :derating_factor,
    :ramp_limits, :max_reserveup, :max_reservedown,
    :existing_units, :units_in_queue, :build_lead_time, :remaining_build_time,
    :max_new_options, :base_cost_units, :capex_years, :life_time, :remaining_life,
    :rec_eligible, :capacity_eligible, :zone])
    life_range = AxisArrays.AxisArray(convert.(Int64, ones(length(projects))), projects)

    tech_types = unique(tech_type)

    option_projects = String[]
    decided_projects = String[]

    option_projects_by_type = Dict(t => Vector{String}() for t in tech_types)

    for g in projects
        if sum(units_in_queue[g]) <= 0
            push!(option_projects, g)

            for t in tech_types
                if tech_type[g] == t
                    push!(option_projects_by_type[t], g)
                end
            end

            if remaining_life_time[g] >= length(invperiods)
                life_range[g] = length(invperiods) - 1
            else
                life_range[g] = remaining_life_time[g] - 1
            end
        else
            push!(decided_projects, g)
        end
    end

    social_discount_rate = maximum(project_discount_rate)   # Social discount rate is kept as the max of all projects' discount rates

    social_npv_array = [(1 / (1 + social_discount_rate)) ^ p for p in invperiods]

    α = AxisArrays.AxisArray(zeros(length(projects)), projects)
    adjusted_investment_cost = AxisArrays.AxisArray(zeros(length(projects), length(invperiods)), projects, invperiods)
    adjusted_queue_cost = AxisArrays.AxisArray(zeros(length(projects)), projects)

    annualized_cap_cost = AxisArrays.AxisArray(zeros(length(projects), length(invperiods)), projects, invperiods)

    for g in projects
        α[g] = project_discount_rate[g] / (1 - (1 + project_discount_rate[g]) ^ (- (capex_years[g])))

        for q in 1:length(queue_cost[g])
            adjusted_queue_cost[g] += queue_cost[g][q] * (1 + project_discount_rate[g]) ^ (build_lead_time[g] + (length(queue_cost[g]) - q))
        end
        for p in invperiods
            adjusted_investment_cost[g, p] = expansion_cost[g][p] * (1 + project_discount_rate[g]) ^ (build_lead_time[g] - 1)
            annualized_cap_cost[g, p] = α[g] * (adjusted_investment_cost[g, p])

        end
    end

    capcost_segmentsize, capcost_segmentgrad, capcost_price_points, capcost_numsegments = make_capital_cost_curve(option_projects_by_type,
                                                                                annualized_cap_cost,
                                                                                base_cost_units,
                                                                                max_new_options,
                                                                                capital_cost_multiplier)

    # Populate markets data
    zones = system.zones
    price_cap_e = AxisArrays.AxisArray(zeros(length(zones), length(invperiods)), zones, invperiods)
    demand_e = AxisArrays.AxisArray(zeros(length(zones), length(invperiods), T), zones, invperiods, opperiods)

    demand_ru = AxisArrays.AxisArray(zeros(length(zones), length(invperiods), T), zones, invperiods, opperiods)

    price_cap_rd = AxisArrays.AxisArray(zeros(length(zones), length(invperiods)), zones, invperiods)
    demand_rd = AxisArrays.AxisArray(zeros(length(zones), length(invperiods), T), zones, invperiods, opperiods)

    price_cap_rec = AxisArrays.AxisArray(zeros(length(invperiods)), invperiods)
    rec_requirement = AxisArrays.AxisArray(zeros(length(invperiods)), invperiods)

    for p in invperiods
        price_cap_rec[p] = getproperty(getproperty(system.inv_periods[p], :rec_market), :price_cap)
        rec_requirement[p] = getproperty(getproperty(system.inv_periods[p], :rec_market), :rec_req)
        for z in zones
            price_cap_e[z, p] = getproperty(getproperty(system.inv_periods[p], :energy_market), :price_cap)[z]
            demand_e[z, p, :] = getproperty(getproperty(system.inv_periods[p], :energy_market), :demand)[z, :]

            demand_ru[z, p, :] = getproperty(getproperty(system.inv_periods[p], :reserveup_market), :demand)[z, :]

            price_cap_rd[z, p] = getproperty(getproperty(system.inv_periods[p], :reservedown_market), :price_cap)[z]
            demand_rd[z, p, :] = getproperty(getproperty(system.inv_periods[p], :reservedown_market), :demand)[z, :]
        end
    end

    ru_segmentsize, ru_segmentgrad, ru_price_points, ru_numsegments = make_ORDC_vectors(getproperty.(system.inv_periods, :reserveup_market))

    cap_segmentsize, cap_segmentgrad, cap_price_points, cap_numsegments = make_capacity_demand_vectors(getproperty.(system.inv_periods, :capacity_market))

    capacity_mkt_projects = Vector{String}()
    rps_compliant_projects = Vector{String}()

    generator_projects = Vector{String}()
    storage_projects = Vector{String}()

    zone_projects = Dict(z => Vector{String}() for z in zones)
    zone_storage = Dict(z => Vector{String}() for z in zones)

    ramp_lim_projects = Vector{String}()

    for project in projects

        if ramp_limits[project] !== nothing
            push!(ramp_lim_projects, project)
        end

        for zone in zones
            if location[project] == zone
                push!(zone_projects[zone], project)

                if project_type[project] == "generator"
                    push!(generator_projects, project)
                elseif project_type[project] == "storage"
                    push!(storage_projects, project)
                    push!(zone_storage[zone], project)
                end
            end
        end

        if capacity_eligible[project]
            push!(capacity_mkt_projects, project) # Find Capacity market eligible projects
        end

        if rec_eligible[project]
            push!(rps_compliant_projects, project) # Find REC eligible projects
        end
    end

    # Populate line data
    from_zone,
    to_zone,
    linepowerlimit =
    make_parameter_vector.(
    Ref(system.lines), :name,
    [:from_zone, :to_zone, :active_power_limit]
    )

    lines_from_zone = Dict(z => Vector{String}() for z in zones)
    lines_to_zone = Dict(z => Vector{String}() for z in zones)

    for line in lines
        for zone in zones
            if from_zone[line] == zone
                push!(lines_from_zone[zone], line)
            elseif to_zone[line] == zone
                push!(lines_to_zone[zone], line)
            end
        end
    end

    m = JuMP.Model(solver)
    # INVESTMENT DECISION VARIABLES

    JuMP.@variable(m, n[g in projects, p in invperiods] >= 0)  # New units built
    JuMP.@variable(m, n_by_type[type in tech_types, p in invperiods, s in 1:capcost_numsegments[type, p]] >= 0) # New units by technology type for capital cost multiplier
    JuMP.@variable(m, r[g in projects, p in invperiods] >= 0)  # Units retired
    JuMP.@variable(m, c[g in projects, p in invperiods] >= 0)  # Unit cleared in capacity market
    JuMP.@variable(m, temp == 0) # For avoiding add_to_expression errors

    # ORDC DEMAND VARIABLE
    JuMP.@variable(m, d_ru[z in zones, p in invperiods, t in opperiods, s in 1:ru_numsegments[z, p]] >= 0)

    # CAPACITY MARKET ELASTIC DEMAND
    JuMP.@variable(m, d_cap[p in invperiods, s in 1:cap_numsegments[p]] >= 0)

    # OPERATION DECISION VARIABLES

    JuMP.@variable(m, p_e[g in projects, p in invperiods, t in opperiods] >= 0) # Unit energy production [MW]
    JuMP.@variable(m, p_in[g in storage_projects, p in invperiods, t in opperiods] >= 0) # Storage charging [MW]
    JuMP.@variable(m, p_ru[g in projects, p in invperiods, t in opperiods] >= 0) # reserve up provided [MW]
    JuMP.@variable(m, p_rd[g in projects, p in invperiods, t in opperiods] >= 0) # reserve down provided [MW]

    JuMP.@variable(m, p_in_ru[g in storage_projects, p in invperiods, t in opperiods] >= 0) # reserve up provided by storage charging [MW]
    JuMP.@variable(m, p_in_rd[g in storage_projects, p in invperiods, t in opperiods] >= 0) # reserve down provided by storage charging [MW]

    JuMP.@variable(m, p_out_ru[g in storage_projects, p in invperiods, t in opperiods] >= 0) # reserve up provided by storage discharging [MW]
    JuMP.@variable(m, p_out_rd[g in storage_projects, p in invperiods, t in opperiods] >= 0) # reserve down provided by storage discharging [MW]

    JuMP.@variable(m, flow[l in lines, p in invperiods, t in opperiods])     # Line flow

    JuMP.@variable(m, v_e[z in zones, p in invperiods, t in opperiods] >= 0) # Load shortfall [MW]
    JuMP.@variable(m, v_rd[z in zones, p in invperiods, t in opperiods] >= 0) # reserve down shortfall [MW]
    JuMP.@variable(m, v_rec[p in invperiods] >= 0) # REC shortfall [MWh]

    # Populate units dispatchable expression based on remaining build/construction time
    JuMP.@expression(m, unitsdispatchable[g in projects, p in invperiods], 0
                            + sum(n[g, i] for i in 1:(p - remaining_buildtime[g]))
                           - r[g, p])

    JuMP.@expression(m, unitsmustretire[g in projects, p in invperiods], 0 + temp)

    JuMP.@expression(m, yearly_cap_cost_decided[g in decided_projects, p in invperiods], 0 + temp)

    JuMP.@expression(m, yearly_cap_cost_options[type in tech_types, p in invperiods], 0 + temp ^ 2)

    JuMP.@expression(m, yearly_queue_cost[g in projects, p in invperiods], 0 + temp)

    # Storage level evolution
    JuMP.@expression(m, storage_level[g in storage_projects, p in invperiods, t in opperiods],
                    p_in[g, p, t] * efficiency_in[g] - p_e[g, p, t] / efficiency_out[g])

    for g in projects

        min_life_time = min(life_time[g], remaining_life_time[g])
        remaining_queuetime = max(0, remaining_buildtime[g] - build_lead_time[g])
        total_queuetime = length(queue_cost[g])

        # Units must retire after their life_time
        for p in invperiods[1:end - min_life_time]
            JuMP.add_to_expression!(unitsmustretire[g, p + min_life_time], unitsdispatchable[g, p] + r[g, p])
        end

        # Calculate yearly capital costs allowing different technology costs each year
        for p in invperiods
            start_year = max(1, p - remaining_buildtime[g] - min_life_time + 1)

                for i in start_year:(p - remaining_buildtime[g])

                    JuMP.add_to_expression!(yearly_queue_cost[g, p], α[g] * adjusted_queue_cost[g] * n[g, i])
                    if in(g, decided_projects)
                        JuMP.add_to_expression!(yearly_cap_cost_decided[g, p], n[g, i] * annualized_cap_cost[g, i + remaining_queuetime])
                    end

                end
            if in(g, storage_projects)
                for t in opperiods
                    if t == 1 || in(t, end_of_day + ones(length(end_of_day)))
                        JuMP.add_to_expression!(storage_level[g, p, t], init_storage[g] * max_storage[g] * unitsdispatchable[g, p])
                    else
                        JuMP.add_to_expression!(storage_level[g, p, t], storage_level[g, p, t - 1])
                    end
                end
            end

        end

    end

    for type in tech_types
        if length(option_projects_by_type[type]) >= 1
            g = option_projects_by_type[type][1]
            min_life_time = min(life_time[g], remaining_life_time[g])
            remaining_queuetime = max(0, remaining_buildtime[g] - build_lead_time[g])
            for p in invperiods
                start_year = max(1, p - remaining_buildtime[g] - min_life_time + 1)
                for i in start_year:(p - remaining_buildtime[g])
                    JuMP.add_to_expression!(yearly_cap_cost_options[type, p],
                        sum(n_by_type[type, i, s] * capcost_price_points[type, i + remaining_queuetime][s] for s in 1:capcost_numsegments[type, i])
                        + 1/2 * sum(capcost_segmentgrad[type, i + remaining_queuetime][s] * (n_by_type[type, i, s]^2) for s in 1:capcost_numsegments[type, i]))
                end
            end
        end
    end

    JuMP.@expression(m, annual_project_costs[g in projects, p in invperiods], 0
                          + yearly_queue_cost[g, p]                               # yearly queue_cost
                          + fixed_cost[g] * unitsdispatchable[g, p]                # fixed operating costs
                          + sum((p_e[g, p, t] * marginal_energy_cost[g]               # project operating costs
                              + p_ru[g, p, t] * marginal_reserveup_cost[g]
                              + p_rd[g, p, t] * marginal_reservedown_cost[g]) * rep_hour_weight[t] for t in opperiods))

    # OBJECTIVE FUNCTION

    JuMP.@objective(m, Max,

    # Capacity market welfare
    sum(sum(d_cap[p, s] * cap_price_points[p][s] for s in 1:cap_numsegments[p]) * social_npv_array[p] for p in invperiods)
    + 1/2 * sum(sum(cap_segmentgrad[p][s] * d_cap[p, s]^2 for s in 1:cap_numsegments[p]) * social_npv_array[p] for p in invperiods)

    # ORDC welfare
    + sum(sum(d_ru[z, p, t, s] * demand_ru[z, p, t] * rep_hour_weight[t] * ru_price_points[z, p][s] for z in zones, t in opperiods, s in 1:ru_numsegments[z, p])
                * social_npv_array[p] for p in invperiods)
    + 1/2 * sum(sum(ru_segmentgrad[z, p][s] * (d_ru[z, p, t, s] * demand_ru[z, p, t])^2 * rep_hour_weight[t] for z in zones, t in opperiods, s in 1:ru_numsegments[z, p])
                * social_npv_array[p] for p in invperiods)

    # REC market Alternative Compliance Payment
    - sum(v_rec[p] * price_cap_rec[p] * social_npv_array[p] for p in invperiods)

    - sum(yearly_cap_cost_decided[g, p] * social_npv_array[p] for g in decided_projects, p in invperiods)

    - sum(yearly_cap_cost_options[type, p] * social_npv_array[p] for type in tech_types, p in invperiods)

    # Annual costs incurred by the project
    - sum(annual_project_costs[g, p] * social_npv_array[p]
        for g in projects, p in invperiods)

    +  (# weigh operations according to actual duration
    # operations market shortfall penalties
             - sum((v_e[z, p, t] * price_cap_e[z, p] + v_rd[z, p, t] * price_cap_rd[z, p]) * rep_hour_weight[t] * social_npv_array[p]
                for z in zones, p in invperiods, t in opperiods)

                        )
        )

    # CONSTRAINTS

    JuMP.@constraint(m, maxoptions[g in projects],
        sum(n[g, p] for p in invperiods) - max_new_options[g] <= 0) # Maximum options purchased

    for p in invperiods

        for type in tech_types
            JuMP.@constraint(m, sum(n_by_type[type, p, s] for s in 1:capcost_numsegments[type, p]) == sum(n[g, p] for g in option_projects_by_type[type]))

            for s in 1:capcost_numsegments[type, p]
                JuMP.@constraint(m, n_by_type[type, p, s] <= capcost_segmentsize[type, p][s])
            end

        end
        # Planning Constraints
        for g in projects
            JuMP.@constraint(m, unitsmustretire[g, p] - r[g, p] <= 0) # Must retire units

            if p < length(invperiods)
                JuMP.@constraint(m, r[g, p] - r[g, p + 1] <= 0) # Units should not become dispatchable after being retired
            end

            if p >= remaining_life_time[g] + 1
                JuMP.@constraint(m, unitsdispatchable[g, p] == 0) # Units should not be dispatchable beyond remaining life_time
            end

            if in(g, capacity_mkt_projects)
                JuMP.@constraint(m, c[g, p] - unitsdispatchable[g, p] <= 0) # Maximum capacity cleared in the capacity market
            end

        end

        # Capacity Demand constraints
        for s in 1:cap_numsegments[p]
            JuMP.@constraint(m, d_cap[p, s] <= cap_segmentsize[p][s]) # Capacity market demand curve segment limit
        end

        for t in opperiods
            # Unit operating constraints (dispatch, reserves and ramping)
            for g in projects
                if in(g, generator_projects)
                    # Generator technical constraints
                    JuMP.@constraint(m, p_e[g, p, t] - p_rd[g, p, t] >= unitsdispatchable[g, p] * 0.0) # Minimum dispatch
                    JuMP.@constraint(m, p_e[g, p, t] + p_ru[g, p, t] <= unitsdispatchable[g, p] * max_gen[g] * availability[g][p, t]) # Maximum dispatch
                    JuMP.@constraint(m, p_rd[g, p, t] <= unitsdispatchable[g, p] * max_reservedown[g]) # Maximum reserve down
                    JuMP.@constraint(m, p_ru[g, p, t] <= unitsdispatchable[g, p] * max_reserveup[g]) # Maximum reserve up

                elseif in(g, storage_projects)
                    # Storage technical constraints
                    JuMP.@constraint(m, p_e[g, p, t] - p_out_rd[g, p, t] >= unitsdispatchable[g, p] * 0.0) # Minimum output dispatch
                    JuMP.@constraint(m, p_e[g, p, t] + p_out_ru[g, p, t] <= unitsdispatchable[g, p] * max_gen[g]) # Maximum output dispatch

                    JuMP.@constraint(m, p_in[g, p, t] - p_in_ru[g, p, t] >= unitsdispatchable[g, p] * 0.0) # Minimum input dispatch
                    JuMP.@constraint(m, p_in[g, p, t] + p_in_rd[g, p, t] <= unitsdispatchable[g, p] * max_input[g]) # Maximum input dispatch

                    JuMP.@constraint(m, p_ru[g, p, t] <= p_out_ru[g, p, t] + p_in_ru[g, p, t]) # Total storage reserve up
                    JuMP.@constraint(m, p_rd[g, p, t] <= p_out_rd[g, p, t] + p_in_rd[g, p, t]) # Total storage reserve down

                    JuMP.@constraint(m, storage_level[g, p, t] >= unitsdispatchable[g, p] * min_storage[g]) # Minimum storage level
                    JuMP.@constraint(m, storage_level[g, p, t] <= unitsdispatchable[g, p] * max_storage[g]) # Maximum storage level

                    if in(t, end_of_day)
                        JuMP.@constraint(m, storage_level[g, p, t] == init_storage[g] * max_storage[g] * unitsdispatchable[g, p])
                    end
                end
            end

            # Line flow constraints
            for l in lines
                JuMP.@constraint(m, flow[l, p, t] >= -linepowerlimit[l])
                JuMP.@constraint(m, flow[l, p, t] <= linepowerlimit[l])
            end

            for z in zones
                for s in 1:ru_numsegments[z, p]
                    JuMP.@constraint(m, d_ru[z, p, t, s] <= ru_segmentsize[z, p][s]) #ORDC segment limit
                end
            end
        end
    end

    # Market Clearing Constraints:

    # Energy market
    JuMP.@constraint(m, energy_market[z in zones, p in invperiods, t in opperiods], # Power balance
        sum(p_e[g ,p, t] for g in zone_projects[z])
        + sum(flow[l, p, t] for l in lines_to_zone[z])
        - sum(flow[l, p, t] for l in lines_from_zone[z])
        + v_e[z, p, t] >= demand_e[z, p, t] + sum(p_in[g, p, t] for g in zone_storage[z]))

    # reserve up market
    JuMP.@constraint(m, reserve_up_market[z in zones, p in invperiods, t in opperiods], # reserve up balance
        sum(p_ru[g ,p, t] for g in zone_projects[z]) == sum(d_ru[z, p, t, s] for s in 1:ru_numsegments[z, p]) * demand_ru[z, p, t] )

    # reserve down market
    JuMP.@constraint(m, reserve_down_market[z in zones, p in invperiods, t in opperiods], # reserve down balance
        sum(p_rd[g ,p, t] for g in zone_projects[z]) + v_rd[z, p, t]  == demand_rd[z, p, t])

    # Capacity market
    JuMP.@constraint(m, capacity_market[p in invperiods],
        sum(c[g, p] * max_gen[g] * derating_factor[g] for g in capacity_mkt_projects) == sum(d_cap[p, s] for s in 1:cap_numsegments[p]))

    # REC market
    JuMP.@constraint(m, rps_compliance[p in invperiods],
        sum(p_e[g, p, t] * rep_hour_weight[t] for g in rps_compliant_projects, t in opperiods) + v_rec[p]
        >= sum(demand_e[z, p, t] * rep_hour_weight[t] for z in zones, t in opperiods) * rec_requirement[p])

    println("Price Projection:")
    @time JuMP.optimize!(m)
    println(JuMP.termination_status(m))
    println(JuMP.objective_value(m))

    #Post processing--------------------------------------------------------------------------------------------------------
    capacity_price = JuMP.dual.(capacity_market)
    energy_price = JuMP.dual.(energy_market)
    reserve_up_price = JuMP.dual.(reserve_up_market)
    reserve_down_price = JuMP.dual.(reserve_down_market)
    REC_price = JuMP.dual.(rps_compliance)

    capacity_factor = Dict([g => zeros(length(invperiods), length(opperiods)) for g in projects])
    capacity_accepted_perc = Dict([g => zeros(length(invperiods)) for g in projects])

    # Number of new options invested in the first year of the CEM.
    new_options = Dict([g => 0 for g in option_projects])

    ϵ = 2e-1

    for g in projects
        for p in invperiods
            if value.(unitsdispatchable[g, p]) < 0.0001
                capacity_accepted_perc[g][p] = 0.0
            else
                capacity_accepted_perc[g][p] = value.(c[g, p]) / value.(unitsdispatchable[g, p])
            end
            for t in opperiods
                if value.(unitsdispatchable[g, p]) < 0.0001
                    capacity_factor[g][p, t] = 0.0
                else
                    capacity_factor[g][p, t] = value.(p_e[g, p, t]) /(value.(unitsdispatchable[g, p]) * max_gen[g])
                end
            end

        end

        if in(g, option_projects)
            for i in 0:base_cost_units[g]:max_new_options[g] + base_cost_units[g]
                if value.(n[g, 1]) > i + ϵ
                    new_options[g] += 1
                end
            end
        end
    end

    nominal_capacity_price = AxisArrays.AxisArray(zeros(length(capacity_price)), invperiods)
    nominal_REC_price = AxisArrays.AxisArray(zeros(length(REC_price)), invperiods)

    nominal_energy_price = AxisArrays.AxisArray(zeros(size(energy_price)), zones, invperiods, opperiods)
    nominal_reserve_up_price = AxisArrays.AxisArray(zeros(size(reserve_up_price)), zones, invperiods, opperiods)
    nominal_reserve_down_price = AxisArrays.AxisArray(zeros(size(reserve_down_price)), zones, invperiods, opperiods)

    for p in invperiods
        nominal_capacity_price[p] = capacity_price[p] / social_npv_array[p]
        nominal_REC_price[p] = REC_price[p] / social_npv_array[p]
        for t in opperiods
            nominal_energy_price[:, p, t] = (energy_price[:, p, t] / rep_hour_weight[t]) / social_npv_array[p]
            nominal_reserve_up_price[:, p, t] = (reserve_up_price[:, p, t] / rep_hour_weight[t]) / social_npv_array[p]
            nominal_reserve_down_price[:, p, t] = (reserve_down_price[:, p, t] / rep_hour_weight[t]) / social_npv_array[p]
        end
    end

    for z in zones
        for p in invperiods
            for t in opperiods
                if nominal_energy_price[z, p, t] <= 1e-5
                    nominal_energy_price[z, p, t] = 1e-5 # replace zeros with 1e-5, to avoid no energy production by renewables when price is 0.
                end
            end
        end
    end

    return nominal_capacity_price,
           nominal_energy_price,
           nominal_reserve_up_price,
           nominal_reserve_down_price,
           nominal_REC_price,
           capacity_factor,
           capacity_accepted_perc,
           new_options;

end


text:
using ToyPkg
using Test

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


text:
using BSON: @save, @load

"""
  train_model()

This function takes no arguments. All the variables used in here
are to be `globally initialized`.

This function trains the `model` using the `data_dict["train"]`
dataset. I call to the function essentially signifies `1 epoch`.

Apart from these variable it expects `save_interval` and
`overwrite_save_interval` to be set prior to the call. THese essentially
help to automate the checkpointing the model.

At the end of the epoch the total loss and accuracy is printed but
if more detailed printing is needed then simply set `verbose` o `1`.

Also the `cost_metric` and `accuracy_metric` dictionaries are filled
up at the end of 1 epoch.
"""

function train_model()
  global model
  info("Starting to train model")
  start_time = time()
  i = 1
  local costs = []
  local accs = []
  for d in data_dict["train"]
    forward = model(d[1] |> gpu)
    res = d[2] |> gpu
    l = loss(res, forward, "train")
    a = accuracy(res, forward)
    push!(accs, Tracker.data(a))
    if verbose == 1
      @show l
      @show a
    end
    if isnan(l) || isinf(l)
      continue # Donot backprop through nan or inf loss instead neglect it
      # This happens since the network is absolutely sertain about its prediction
      # So the loss becomes NaN as it doesn't handle this case
      # error("Model Weights have diverged. Change the Learning Rate")
    end
    push!(costs, Tracker.data(l))
    Flux.back!(l)
    opt()
    end_time = time()
    if(end_time - start_time > save_interval)
      model = model |> cpu
      i == overwrite_save_interval ? i = 1 : i += 1
      save_path = model_save_path * "_i.bson"
      @save save_path model
      info("Model Saved at $(save_path)")
      model = model |> gpu
      start_time = time()
    end
  end
  push!(cost_metric["train"], sum(costs)/length(data_dict["train"]))
  push!(accuracy_metric["train"], sum(accs)/length(data_dict["train"]))
  info("Training Loss is $(cost_metric["train"][end]) Training Accuracy is $(accuracy_metric["train"][end])")
end

"""
  validate_model()

This function takes no arguments. All the variables used in here
are to be `globally initialized`.

This function tests the `model` using the `data_dict["valid"]`
dataset.

Also the `cost_metric` and `accuracy_metric` dictionaries are filled
up at the end. Additionally it displays the validation loss and
accuracy.
"""

function validate_model()
  global model
  info("Validating Model")
  Flux.testmode!(model)
  local costs = []
  local accs = []
  for d in data_dict["valid"]
    forward = model(d[1] |> gpu)
    res = d[2] |> gpu
    l = loss(res, forward, "train")
    a = accuracy(res, forward)
    push!(accs, Tracker.data(a))
    if isnan(l) || isinf(l)
      continue # Donot backprop through nan or inf loss instead neglect it
    end
    push!(costs, Tracker.data(l))
  end
  push!(cost_metric["valid"], sum(costs)/length(data_dict["valid"]))
  push!(accuracy_metric["valid"], sum(accs)/length(data_dict["valid"]))
  Flux.testmode!(model, false)
  info("Validation Loss is $(cost_metric["valid"][end]) Validation Accuracy is $(accuracy_metric["valid"][end])")
end


text:
import Base: display, length

struct ChequeredBlocking <: BlockingSchedule end

struct Block{S1,S2,S3,S4}
    Ls::Vector{S1}         # observation operators for both sets of blocks
    vs::Vector{S2}         # copied over observations of the process
    Σs::Vector{S3}         # covariance matrix of the noise for both sets of blocks
    aux_flags::Vector{S4}
    # two sets of knots
    knots::Vector{Int64}
    # two sets of blocks, where each block are indicies of intervals that make up a block
    blocks::Vector{Vector{Int64}}
    # info about the points at which to switch between the systems of ODEs
    change_pts::Vector{ODEChangePt}

    function Block(Ls::Vector{S1}, vs::Vector{S2}, Σs::Vector{S3},
                   aux_flags::Vector{S4}, knots::Vector{Int64},
                   blocks::Vector{Vector{Int64}},
                   change_pts::Vector{ODEChangePt}) where {S1,S2,S3,S4}
        new{S1,S2,S3,S4}(Ls, vs, Σs, aux_flags, knots, blocks, change_pts)
    end
end

length(bl::Block) = length(bl.blocks)
length(::NoBlocking) = 1

"""
    find_end_pts(𝔅::ChequeredBlocking, XX, idx)

Determine the observations for the update of the `idx`-th set of blocks. In
particular, on each block with interval indices [a₁,…,aₙ], observations vᵢ with
i∈{1,…,n-1} are made, whereas the last obesrvation is an exactly observed
process Xₜ at the terminal time t of the aₙ-th interval.
"""
function find_end_pts(bl::Block, XX)
    [( k in bl.knots ? X.yy[end] : bl.vs[k]) for (k,X) in enumerate(XX)]
end

"""
    display(𝔅::NoBlocking)

Nothing particular to display
"""
function display(bl::NoBlocking)
    print("No blocking...\n")
end

"""
    display(𝔅::ChequeredBlocking)

Display the pattern of blocks
"""
function display(bl::Block)
    function print_blocks(knots, idxLast, m)
        M = length(knots)
        get_knot(knots, i) = (M >= i > 0 ? knots[i] : idxLast * (i>0))
        function print_range(from, to)
            for i in from:to
                print("|", get_knot(knots,i-1), "|----",
                      get_knot(knots,i)-get_knot(knots,i-1), "----")
            end
        end

        if m > div(M, 2)
            print_range(1, M+1)
            print("|", idxLast, "|")
        else
            print_range(1,m)
            print("|", get_knot(knots, m) ,"|   ...   ")
            print_range(M+2-m,M+1)
            print("|", idxLast, "|")
        end
        print("  (number of blocks: ", M+1,")")
    end
    print("Chequered Blocking scheme\n",
          "-------------------------\n",
          "Format:\n",
          "block sizes in A: ")
    print_blocks(𝔅.knots[1], length(𝔅.vs), 3)
    print("\nblock sizes in B: ")
    print_blocks(𝔅.knots[2], length(𝔅.vs), 3)
    print("\n")
end

function _print_info(info)
    m = length(info)
    for i in 1:m
        print("b", i, ": ", round(info[i], digits=3), " | ")
    end
    print("\n- - - - - - - - - - - - - -\n")
end

function create_blocks(::ChequeredBlocking, P, params)
    knots, ϵ, change_pt = params
    find_knots(mod, rem) = [k for (i,k) in enumerate(knots) if i % mod == rem]
    knotsA, knotsB = find_knots(2, 1), find_knots(2, 0)

    m, d = size(P[end].L)

    find_L(knots) = [( k in knots ? SMatrix{d,d}(1.0*I) : p.L)
                     for (k,p) in enumerate(P)]
    LsA, LsB = find_L(knotsA), find_L(knotsB)

    vs = [p.v for p in P]

    find_Σ(knots) = [(k in knots ? SMatrix{d,d}(ϵ*I) : p.Σ)
                     for (k,p) in enumerate(P)]
    ΣsA, ΣsB = find_Σ(knotsA), find_Σ(knotsB)

    find_ch_pt(knots) = [(k in knots ? deepcopy(change_pt) : p.change_pt)
                         for (k,p) in enumerate(P)]
    chpA, chpB = find_ch_pt(knotsA), find_ch_pt(knotsB)

    find_aux_flag(knots) = [( k in knots ? Nothing : get_aux_flag(p.Pt))
                            for (k,p) in enumerate(P)]
    auxA, auxB = find_aux_flag(knotsA), find_aux_flag(knotsB)

    """
        knots_to_blocks(knots, idxLast, i)

    Given a list of `knots` fetch the indices of the intervals that together
    make up the `i`-th block. `idxLast` is the index of the last interval.
    """
    function knots_to_blocks(knots, idx_last, i)
        M = length(knots)
        @assert M > 0
        if M >= i > 1
            return (knots[i-1]+1):knots[i]
        elseif i == 1
            return 1:knots[1]
        else
            return (knots[M]+1):idx_last
        end
    end
    blocks = ([collect(knots_to_blocks(knotsA, length(P), i)) for i in 1:length(knotsA)+1],
              [collect(knots_to_blocks(knotsB, length(P), i)) for i in 1:length(knotsB)+1])

    (Block(LsA, vs, ΣsA, auxA, knotsA, blocks[1], chpA),
     Block(LsB, vs, ΣsB, auxB, knotsB, blocks[2], chpB))
end

create_blocks(blocking::NoBlocking, P, params) = blocking


text:
#=
    Test all my constraint models in ConstraintSolver.jl

    This program test all

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

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
using HTTP

# Skip testing these files
skip = Dict(
    "test_all_models.jl" => true, # Well, we run it just one, right. :-)
    "constraints_utils.jl" => true,
    "minesweeper_problems.jl" => true,
    "sudoku_problems.jl" => true,
    "survo_puzzle_problems.jl" => true,
    "quasigroup_completion_problems.jl" => true,
    "assignment_instances.jl" => true,

)

#
# Read all the public models from my http://hakank.org/julia/constraints/index.html
#
function read_from_hakank_org()
    url = "http://hakank.org/julia/constraints/index.html"
    res = HTTP.get(url)
    body = String(res.body)
    regex = r"<a href=\"([^\s]+?.(?:jl))\""sm
    all = []
    for t in [t[1] for t in eachmatch(regex,body)]
        if get(skip,t,false)
            println("Skipping $t")
            continue
        end
        push!(all, t)
    end
    return all
end

# Test all my public ConstraintSolver.jl models
function test_all()

    files = read_from_hakank_org()

    total_time = 0
    for file in files
        println("\nRunning $file")
        @time include(file)

    end

end

test_all()
# read_from_hakank_org()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#=

File defines the abstract types used within the package.

=#

abstract type AbstractProblem{T} end

abstract type AbstractPoll end

abstract type AbstractSearch end

"""
    abstract type AbstractMesh end

Parent type of any struct implementing the construction of a mesh. To maintain
compatibility with other aspects of the package, the naming convention for
variables within structs must be followed. These respect the notation used
within Audet, Le Digabel & Tribes 2019.
"""
abstract type AbstractMesh end

abstract type AbstractConstraint end

"""
    abstract type AbstractCache end

Parent type of any struct implementing the cache.
"""
abstract type AbstractCache end

"""
    abstract type AbstractStoppingCondition end

Parent type of any struct implementing a stopping condition.
"""
abstract type AbstractStoppingCondition end

@enum ProblemSense Min Max

@enum OptimizationStatus Unoptimized MeshPrecisionLimit PollPrecisionLimit IterationLimit FunctionEvaluationLimit RuntimeLimit KeyInterrupt HypervolumeLimit OtherStoppingCondition

"""
    Config{FT}(N::Int,
           poll::AbstractPoll,
           search::AbstractSearch,
           mesh::AbstractMesh=Mesh{FT}(N);
           opportunistic::Bool=false,
           cost_digits::Int=32,
           kwargs...
           ) where {FT<:AbstractFloat}

Encapsulates configuration options for the solver, generally shouldn't
be user-edited.

Generally these are set at the start (automatically or via setter functions)
and don't change.
"""
mutable struct Config{FT<:AbstractFloat, MT<:AbstractMesh, ST<:AbstractSearch, PT<:AbstractPoll}

    poll::PT
    search::ST

    mesh::MT

    num_threads::Int
    max_simultaneous_evaluations::Int
    parallel_lock::Threads.AbstractLock
    opportunistic::Bool

    cost_digits::Int

    function Config{FT}(N::Int,
                        poll::AbstractPoll,
                        search::AbstractSearch,
                        mesh::AbstractMesh=Mesh{FT}(N);
                        opportunistic::Bool=false,
                        cost_digits::Int=32,
                        kwargs...
                       ) where {FT<:AbstractFloat}
        c = new{FT, typeof(mesh), typeof(search), typeof(poll)}()

        c.poll = poll
        c.search = search

        c.mesh = mesh

        c.num_threads = Threads.nthreads()
        c.max_simultaneous_evaluations = 1
        c.opportunistic = opportunistic

        c.cost_digits = cost_digits

        return c
    end
end

"""
    Status{T}() where T

Holds the status information of the solver.
"""
mutable struct Status{T}
    function_evaluations::Int64
    cache_hits::Int64
    iteration::Int64
    directions::Union{Vector{Vector{T}}, Nothing}
    success_direction::Union{Vector{T},Nothing}
    # paretoCoverage::Float64

    optimization_status::OptimizationStatus
    optimization_status_string::String

    #= Time Running Totals =#
    runtime_total::Float64
    search_time_total::Float64
    poll_time_total::Float64
    blackbox_time_total::Float64

    #= Start/End Time =#
    start_time::Float64
    end_time::Float64

    function Status{T}() where T
        s = new()

        s.function_evaluations = 0
        s.iteration = 0
        s.cache_hits = 0
        s.directions = nothing
        s.success_direction = nothing
        s.optimization_status_string = "Unoptimized"
        s.optimization_status = Unoptimized
        # s.paretoCoverage=0.0

        s.runtime_total = 0.0
        s.search_time_total = 0.0
        s.poll_time_total = 0.0
        s.blackbox_time_total = 0.0

        return s
    end
end


text:
module JUtils

using Colors
using Animations
using Javis
using Random
using FileIO
using StatsBase
using Distributions
using ROCAnalysis
using LaTeXStrings

# import Luxor

const φ = Base.MathConstants.φ

include("Drawings.jl")
include("Objects.jl")
include("ActingObjects.jl")

export load_and_scale, draw_thickbox
export moving_object, draw_random_point
export draw_text, make_text, moving_text
export draw_continuous_lines!
export draw_lamp, make_lamp, draw_lens, make_lens
export draw_doc, make_doc, draw_banner, make_banner
export draw_roc

include("Videos.jl")
export CLT

include("Util.jl")
export sensitivity, specificity

include("NER.jl")
export make_ner

include("prelearn.jl")
export make_prelearn

end


text:
using ModelingToolkit, StaticArrays, LinearAlgebra, SparseArrays
using DiffEqBase
using Test

canonequal(a, b) = isequal(simplify(a), simplify(b))

# Calculus
@parameters t σ ρ β
@variables x y z

eqs = [σ*(y-x),
       x*(ρ-z)-y,
       x*y - β*z]

simpexpr = [
   :(σ * (y - x))
   :(x * (ρ - z) - y)
   :(x * y - β * z)
   ]

for i in 1:3
   @test ModelingToolkit.simplified_expr.(eqs)[i] == simpexpr[i]
   @test ModelingToolkit.simplified_expr.(eqs)[i] == simpexpr[i]
end

∂ = ModelingToolkit.jacobian(eqs,[x,y,z])
for i in 1:3
    ∇ = ModelingToolkit.gradient(eqs[i],[x,y,z])
    @test canonequal(∂[i,:],∇)
end

@test all(canonequal.(ModelingToolkit.gradient(eqs[1],[x,y,z]),[σ * -1,σ,0]))
@test all(canonequal.(ModelingToolkit.hessian(eqs[1],[x,y,z]),0))

Joop,Jiip = eval.(ModelingToolkit.build_function(∂,[x,y,z],[σ,ρ,β],t))
J = Joop([1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test J isa Matrix
J2 = copy(J)
Jiip(J2,[1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test J2 == J

Joop,Jiip = eval.(ModelingToolkit.build_function(vcat(∂,∂),[x,y,z],[σ,ρ,β],t))
J = Joop([1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test J isa Matrix
J2 = copy(J)
Jiip(J2,[1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test J2 == J

Joop,Jiip = eval.(ModelingToolkit.build_function(hcat(∂,∂),[x,y,z],[σ,ρ,β],t))
J = Joop([1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test J isa Matrix
J2 = copy(J)
Jiip(J2,[1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test J2 == J

∂3 = cat(∂,∂,dims=3)
Joop,Jiip = eval.(ModelingToolkit.build_function(∂3,[x,y,z],[σ,ρ,β],t))
J = Joop([1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test size(J) == (3,3,2)
J2 = copy(J)
Jiip(J2,[1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test J2 == J

s∂ = sparse(∂)
@test nnz(s∂) == 8
Joop,Jiip = eval.(ModelingToolkit.build_function(s∂,[x,y,z],[σ,ρ,β],t,linenumbers=true))
J = Joop([1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test length(nonzeros(s∂)) == 8
J2 = copy(J)
Jiip(J2,[1.0,2.0,3.0],[1.0,2.0,3.0],1.0)
@test J2 == J

# Function building

@parameters σ() ρ() β()
@variables x y z
eqs = [σ*(y-x),
       x*(ρ-z)-y,
       x*y - β*z]
f1,f2 = ModelingToolkit.build_function(eqs,[x,y,z],[σ,ρ,β])
f = eval(f1)
out = [1.0,2,3]
o1 = f([1.0,2,3],[1.0,2,3])
f = eval(f2)
f(out,[1.0,2,3],[1.0,2,3])
@test all(o1 .== out)

function test_worldage()
   @parameters σ() ρ() β()
   @variables x y z
   eqs = [σ*(y-x),
          x*(ρ-z)-y,
          x*y - β*z]
   f, f_iip = ModelingToolkit.build_function(eqs,[x,y,z],[σ,ρ,β];expression=Val{false})
   out = [1.0,2,3]
   o1 = f([1.0,2,3],[1.0,2,3])
   f_iip(out,[1.0,2,3],[1.0,2,3])
end
test_worldage()

## No parameters
@variables x y z
eqs = [(y-x)^2,
       x*(x-z)-y,
       x*y - y*z]
f1,f2 = ModelingToolkit.build_function(eqs,[x,y,z])
f = eval(f1)
out = zeros(3)
o1 = f([1.0,2,3])
f = eval(f2)
f(out,[1.0,2,3])
@test all(out .== o1)

function test_worldage()
   @variables x y z
   eqs = [(y-x)^2,
          x*(x-z)-y,
          x*y - y*z]
   f, f_iip = ModelingToolkit.build_function(eqs,[x,y,z];expression=Val{false})
   out = zeros(3)
   o1 = f([1.0,2,3])
   f_iip(out,[1.0,2,3])
end
test_worldage()

@test_nowarn muladd(x, y, ModelingToolkit.Constant(0))
@test promote(x, ModelingToolkit.Constant(0)) == (x, identity(0))
@test_nowarn [x, y, z]'


text:
@doc """
An ORDFile object represents the information in an ORD format file like that
found on the VoteView website. It contains information about legislators and
their roll-call votes. The available information includes:

* congresses: To which Congress does the data belong?
* icpsr_ids: What is the ICPSR ID for the legislator?
* icpsr_states: What is the ICPSR state for the legislator?
* cds: What is the congressional district for the legislator?
* kh_states: What is the K/H state for the legislator?
* parties: What is the party for the legislator?
* unknowns: UNKNOWN PURPOSE COLUMN IN ORD FILE
* legislators: What is the name of the legislator?
* roll_call: The full roll-call from the ORD file stored as a RollCall object.

""" ->
immutable ORDFile{T1 <: Tuple, T2 <: Tuple, T3 <: Tuple, T4 <: Tuple}
    congresses::Vector{Int}
    icpsr_ids::Vector{Int}
    icpsr_states::Vector{Int}
    cds::Vector{Int}
    kh_states::Vector{UTF8String}
    parties::Vector{Int}
    unknowns::Vector{UTF8String}
    legislators::Vector{UTF8String}
    roll_call::RollCall{T1, T2, T3, T4}

    function ORDFile(
        congresses::Vector{Int},
        icpsr_ids::Vector{Int},
        icpsr_states::Vector{Int},
        cds::Vector{Int},
        kh_states::Vector{UTF8String},
        parties::Vector{Int},
        unknown::Vector{UTF8String},
        legislators::Vector{UTF8String},
        roll_call::RollCall{T1, T2, T3, T4},
    )
        n_legislators, n_bills = size(roll_call.votes)
        @assert length(congresses) == n_legislators
        @assert length(icpsr_ids) == n_legislators
        @assert length(icpsr_states) == n_legislators
        @assert length(cds) == n_legislators
        @assert length(kh_states) == n_legislators
        @assert length(parties) == n_legislators
        @assert length(unknown) == n_legislators
        @assert length(legislators) == n_legislators
        new(
            congresses,
            icpsr_ids,
            icpsr_states,
            cds,
            kh_states,
            parties,
            unknown,
            legislators,
            roll_call,
        )
    end
end

# Outer constructor necessary because of the inner constructor used for a
# parametric type.
function ORDFile{T1 <: Tuple, T2 <: Tuple, T3 <: Tuple, T4 <: Tuple}(
    congresses::Vector{Int},
    icpsr_ids::Vector{Int},
    icpsr_states::Vector{Int},
    cds::Vector{Int},
    kh_states::Vector{UTF8String},
    parties::Vector{Int},
    unknown::Vector{UTF8String},
    legislators::Vector{UTF8String},
    roll_call::RollCall{T1, T2, T3, T4},
)
    ORDFile{T1, T2, T3, T4}(
        congresses,
        icpsr_ids,
        icpsr_states,
        cds,
        kh_states,
        parties,
        unknown,
        legislators,
        roll_call,
    )
end

# Getter methods
congresses(ord_file::ORDFile) = ord_file.congresses
icpsr_ids(ord_file::ORDFile) = ord_file.icpsr_ids
icpsr_states(ord_file::ORDFile) = ord_file.icpsr_states
cds(ord_file::ORDFile) = ord_file.cds
kh_states(ord_file::ORDFile) = ord_file.kh_states
parties(ord_file::ORDFile) = ord_file.parties
unknowns(ord_file::ORDFile) = ord_file.unknowns
legislators(ord_file::ORDFile) = ord_file.legislators
roll_call(ord_file::ORDFile) = ord_file.roll_call

# TODO: Make this a very low-memory streaming reader.
# NOTE: The use of numeric indexing into substrings only works because the data
#       is assumed to always be ASCII.
@doc """
Use this function to parse an ORDFile from an arbitrary IO object. Assumes that
the votes should be encoding using the default VoteView encoding for 'yea',
'nay', 'missing' and 'inactive' votes. If another encoding should be used,
this must be passed as a third argument.
""" ->
function Base.read(
    ::Type{ORDFile},
    io::IO,
    encoding::RollCallEncoding = RollCallEncoding(),
)
    congresses = Int[]
    icpsr_ids = Int[]
    icpsr_states = Int[]
    cds = Int[]
    kh_states = UTF8String[]
    parties = Int[]
    unknown = UTF8String[]
    legislators = UTF8String[]
    votes = Int[]

    n_legislators = 0
    n_bills = 0

    while !eof(io)
        line = chomp(readline(io))

        n_legislators += 1
        n_bills = length(line) - 36

        # Parse the Congressional session.
        push!(congresses, parse(Int, line[1:3]))

        # Parse the ICPSR legislator code.
        push!(icpsr_ids, parse(Int, line[4:8]))

        # Parse the ICPSR state code.
        push!(icpsr_states, parse(Int, line[9:10]))

        # Parse the congressional district.
        push!(cds, parse(Int, line[11:12]))

        # Parse the Keith/Howard state code.
        push!(kh_states, rstrip(line[13:20]))

        # Parse the party code.
        push!(parties, parse(Int, line[21:23]))

        # TODO: Figure out what this is.
        push!(unknown, line[24:25])

        # Parse the legislator name.
        push!(legislators, rstrip(line[26:36]))

        # Parse the votes for this legislator.
        for j in 1:n_bills
            push!(votes, parse(Int, line[36 + j]))
        end
    end

    return ORDFile(
        congresses,
        icpsr_ids,
        icpsr_states,
        cds,
        kh_states,
        parties,
        unknown,
        legislators,
        RollCall(
            transpose(reshape(votes, n_bills, n_legislators)),
            encoding,
        ),
    )
end

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
function construct_device!(ps_m::CanonicalModel,
                           device::Type{L},
                           device_formulation::Type{D},
                           system_formulation::Type{S},
                           sys::PSY.System,
                           time_range::UnitRange{Int64};
                           kwargs...) where {L <: PSY.ElectricLoad,
                                             D <: PSI.AbstractControllablePowerLoadForm,
                                             S <: PM.AbstractPowerFormulation}

    parameters = get(kwargs, :parameters, true)

    fixed_resources = [fs for fs in sys.loads if isa(fs,PSY.PowerLoad)]

    controllable_resources = [fs for fs in sys.loads if !isa(fs,PSY.PowerLoad)]

        if !isempty(controllable_resources)

            #Variables
            activepower_variables(ps_m, controllable_resources, time_range);

            reactivepower_variables(ps_m, controllable_resources, time_range);

            #Constraints
            activepower_constraints(ps_m, controllable_resources, device_formulation, system_formulation, time_range, parameters)

            reactivepower_constraints(ps_m, controllable_resources, device_formulation, system_formulation, time_range)

            #Cost Function
            cost_function(ps_m, controllable_resources, device_formulation, system_formulation)

        else
            @warn("The Data Doesn't Contain Controllable Loads, Consider Changing the Device Formulation to StaticPowerLoad")

        end

        #add to expression

        if !isempty(fixed_resources)
            nodal_expression(ps_m, fixed_resources, system_formulation, time_range, parameters)
        end

        return

end

function construct_device!(ps_m::CanonicalModel,
                           device::Type{L},
                           device_formulation::Type{D},
                           system_formulation::Type{S},
                           sys::PSY.System,
                           time_range::UnitRange{Int64};
                           kwargs...) where {L <: PSY.ElectricLoad,
                                             D <: PSI.AbstractControllablePowerLoadForm,
                                             S <: PM.AbstractActivePowerFormulation}

    parameters = get(kwargs, :parameters, true)

    fixed_resources = [fs for fs in sys.loads if isa(fs,PSY.PowerLoad)]

    controllable_resources = [fs for fs in sys.loads if !isa(fs,PSY.PowerLoad)]

    if !isempty(controllable_resources)

        #Variables
        activepower_variables(ps_m, controllable_resources, time_range);

        #Constraints
        activepower_constraints(ps_m, controllable_resources, device_formulation, system_formulation, time_range, parameters)

        #Cost Function
        cost_function(ps_m, controllable_resources, device_formulation, system_formulation)

    else
        @warn("The Data Doesn't Contain Controllable Loads, Consider Changing the Device Formulation to StaticPowerLoad")
    end

    #add to expression

    if !isempty(fixed_resources)
        nodal_expression(ps_m, fixed_resources, system_formulation, time_range, parameters)
    end

    return

end

function construct_device!(ps_m::CanonicalModel,
                           device::Type{L},
                           device_formulation::Type{PSI.StaticPowerLoad},
                           system_formulation::Type{S},
                           sys::PSY.System,
                           time_range::UnitRange{Int64};
                           kwargs...) where {L <: PSY.ElectricLoad,
                                             S <: PM.AbstractPowerFormulation}

    parameters = get(kwargs, :parameters, true)

    nodal_expression(ps_m, sys.loads, system_formulation, time_range, parameters)

    return

end


text:
# layout-related functions

ContViews{T,N} =  Union{ContiguousView{T,N},UnsafeContiguousView{T,N}}
NonContViews{T,N,M} =  Union{StridedView{T,N,M}, UnsafeStridedView{T,N,M}}

ContiguousArray{T,N} =  Union{Array{T,N}, ContiguousView{T,N}}
ContiguousVector{T} =  ContiguousArray{T,1}
ContiguousMatrix{T} =  ContiguousArray{T,2}

import Base: IndexStyle, IndexLinear, IndexCartesian
IndexStyle(a::ContViews) = IndexLinear()
IndexStyle(a::NonContViews) = IndexCartesian()
IndexStyle(a::NonContViews{T,1}) where {T} = IndexLinear()

## astrides method

astrides(a::Array{T,N}) where {T,N} = strides(a)

astrides(a::ContViews{T,1}) where {T} = (1,)
astrides(a::ContViews{T,2}) where {T} = (1, a.shp[1])
astrides(a::ContViews{T,3}) where {T} = ((d1, d2, d3) = size(a); (1, d1, d1 * d2))
astrides(a::ContViews{T,4}) where {T} = ((d1, d2, d3, d4) = size(a); (1, d1, d1 * d2, d1 * d2 * d3))

function astrides(a::ContViews{T,5}) where T
    (d1, d2, d3, d4, d5) = a.shp
    s3 = d1 * d2
    s4 = s3 * d3
    s5 = s4 * d4
    (1, d1, s3, s4, s5)
end

astrides(a::NonContViews) = a.astrides

## astride method

astride(a::Array{T,N}, d::Integer) where {T,N} = stride(a, d)

astride(a::ContViews{T,1}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 : length(a))::Int

astride(a::ContViews{T,2}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] : length(a))::Int

astride(a::ContViews{T,3}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] : length(a))::Int

 astride(a::ContViews{T,4}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] :
     d == 4 ? a.shp[1] * a.shp[2] * a.shp[3] : length(a))::Int

astride(a::ContViews{T,5}, d::Integer) where {T} =
    (d > 0 || error("dimension out of range.");
     d == 1 ? 1 :
     d == 2 ? a.shp[1] :
     d == 3 ? a.shp[1] * a.shp[2] :
     d == 4 ? a.shp[1] * a.shp[2] * a.shp[3] :
     d == 5 ? a.shp[1] * a.shp[2] * a.shp[3] * a.shp[4] : length(a))::Int

astride(a::NonContViews{T,N}, d::Integer) where {T,N} =
    (d > 0 || error("dimension out of range.");
     d <= N ? a.astrides[d] : a.astrides[N])

 # test contiguousness

iscontiguous(a::ContViews) = true

iscontiguous(a::NonContViews{T,N}) where {T,N} = _iscontiguous(size(a), astrides(a))

_iscontiguous(shp::NTuple{0,Int}, astrides::NTuple{0,Int}) = true
_iscontiguous(shp::NTuple{1,Int}, astrides::NTuple{1,Int}) = (astrides[1] == 1)
_iscontiguous(shp::NTuple{2,Int}, astrides::NTuple{2,Int}) = (astrides[1] == 1 && astrides[2] == shp[1])
_iscontiguous(shp::NTuple{3,Int}, astrides::NTuple{3,Int}) =
    (astrides[1] == 1 &&
     astrides[2] == shp[1] &&
     astrides[3] == shp[1] * shp[2])

function _iscontiguous(shp::NTuple{4,Int}, astrides::NTuple{4,Int})
    s2 = shp[1]
    s3 = s2 * shp[2]
    s4 = s3 * shp[3]
    astrides[1] == 1 && astrides[2] == s2 && astrides[3] == s3 && astrides[4] == s4
end

function _iscontiguous(shp::NTuple{5,Int}, astrides::NTuple{5,Int})
    s2 = shp[1]
    s3 = s2 * shp[2]
    s4 = s3 * shp[3]
    s5 = s4 * shp[4]
    astrides[1] == 1 && astrides[2] == s2 && astrides[3] == s3 && astrides[4] == s4 && astrides[5] == s5
end


text:
#=

  Euler #41 in Julia.

  Problem 41:
  """
  We shall say that an n-digit number is pandigital if it makes use of all
  the digits 1 to n exactly once. For example, 2143 is a 4-digit pandigital
  and is also prime.

  What is the largest n-digit pandigital prime that exists?
  """

  This Julia model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/Julia_progs/

=#

include("Euler.jl")

# 0.03632750s
function euler41a()
    # Simplification:
    # n=9 is not possible since 1+2+3+4+5+6+7+8+9=45 is divisible by 3
    #  n=8 is not possible since 1+2+3+4+5+6+7+8=36 is divisible by 3
    n = 7
    m = 0
    while m == 0 && n >= 4
        p = Array(1:n)
        for pp in sort(all_permutations(p),rev=true)
            v = parse(Int,join(pp,""))
            if isPrime(v)
                return v
            end
        end
        n += 1
    end

    return nothing
end

run_euler(euler41a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using CUDAapi
using CUDAdrv
using LLVM

## auxiliary routines

function build_error(reason)
    println("""
        $reason.

        This is not a fatal error, but GPU functionality will be unavailable.
        If you expected this to work, please open a thread on
        https://discourse.julialang.org/c/domain/gpu""")
    exit(1)
end

function llvm_support(version)
    @debug("Using LLVM $version")

    InitializeAllTargets()
    haskey(targets(), "nvptx") ||
        build_error("""
            Your LLVM does not support the NVPTX back-end.

            This is very strange; both the official binaries
            and an unmodified build should contain this back-end.""")

    target_support = sort(collect(CUDAapi.devices_for_llvm(version)))

    ptx_support = CUDAapi.isas_for_llvm(version)
    push!(ptx_support, v"6.0") # JuliaLang/julia#23817
    ptx_support = sort(collect(ptx_support))

    @debug("LLVM support", targets=target_support, isas=ptx_support)
    return target_support, ptx_support
end

function cuda_support(driver_version, toolkit_version)
    @debug("Using CUDA driver $driver_version and toolkit $toolkit_version")

    # the toolkit version as reported contains major.minor.patch,
    # but the version number returned by libcuda is only major.minor.
    toolkit_version = VersionNumber(toolkit_version.major, toolkit_version.minor)
    if toolkit_version > driver_version
        build_error("""
            CUDA $(toolkit_version.major).$(toolkit_version.minor) is not supported by
            your driver (which supports up to $(driver_version.major).$(driver_version.minor))""")
    end

    driver_target_support = CUDAapi.devices_for_cuda(driver_version)
    toolkit_target_support = CUDAapi.devices_for_cuda(toolkit_version)
    target_support = sort(collect(driver_target_support ∩ toolkit_target_support))

    driver_ptx_support = CUDAapi.isas_for_cuda(driver_version)
    toolkit_ptx_support = CUDAapi.isas_for_cuda(toolkit_version)
    ptx_support = sort(collect(driver_ptx_support ∩ toolkit_ptx_support))

    @debug("CUDA driver support", version=driver_version,
           targets=driver_target_support, isas=driver_ptx_support)
    @debug("CUDA toolkit support", version=toolkit_version,
           targets=toolkit_target_support, isas=toolkit_ptx_support)

    return target_support, ptx_support
end

## main

# development often happens in lockstep with other packages,
# so check-out the master branch of those packages.
using Pkg
if haskey(ENV, "GITLAB_CI")
  for package in ("CUDAdrv", "LLVM")
    Pkg.add(PackageSpec(name=package, rev="master"))
  end
end

const config_path = joinpath(@__DIR__, "ext.jl")
const previous_config_path = config_path * ".bak"

function write_ext(config, path)
    open(path, "w") do io
        println(io, "# autogenerated file, do not edit")
        for (key,val) in config
            println(io, "const $key = $(repr(val))")
        end
    end
end

function read_ext(path)
    config = Dict{Symbol,Any}()
    r = r"^const (\w+) = (.+)$"
    open(path, "r") do io
        for line in eachline(io)
            m = match(r, line)
            if m != nothing
                config[Symbol(m.captures[1])] = eval(Meta.parse(m.captures[2]))
            end
        end
    end
    return config
end

function main()
    ispath(config_path) && mv(config_path, previous_config_path; force=true)
    config = Dict{Symbol,Any}(:configured => false)
    write_ext(config, config_path)

    ## gather info

    ### LLVM.jl

    LLVM.configured || build_error("Dependent package LLVM.jl has not been built successfully")

    LLVM.libllvm_system && build_error("CUDAnative.jl requires LLVM.jl to be built against Julia's LLVM library, not a system-provided one")

    llvm_version = LLVM.version()
    llvm_targets, llvm_isas = llvm_support(llvm_version)

    ### julia

    julia_llvm_version = Base.libllvm_version
    if julia_llvm_version != llvm_version
        build_error("LLVM $llvm_version incompatible with Julia's LLVM $julia_llvm_version")
    end

    ### CUDA

    CUDAdrv.configured || build_error("Dependent package CUDAdrv.jl has not been built successfully")

    toolkit_dirs = find_toolkit()
    cuda_toolkit_version = find_toolkit_version(toolkit_dirs)

    config[:cuda_driver_version] = CUDAdrv.version()
    cuda_targets, cuda_isas = cuda_support(config[:cuda_driver_version], cuda_toolkit_version)

    config[:target_support] = sort(collect(llvm_targets ∩ cuda_targets))
    isempty(config[:target_support]) && build_error("Your toolchain does not support any device target")

    config[:ptx_support] = sort(collect(llvm_isas ∩ cuda_isas))
    isempty(config[:target_support]) && build_error("Your toolchain does not support any PTX ISA")

    @debug("CUDAnative support", targets=config[:target_support], isas=config[:ptx_support])

    # discover other CUDA toolkit artifacts
    ## required
    config[:libdevice] = find_libdevice(config[:target_support], toolkit_dirs)
    config[:libdevice] == nothing && build_error("Available CUDA toolchain does not provide libdevice")
    ## optional
    config[:nvdisasm] = find_cuda_binary("nvdisasm", toolkit_dirs)
    config[:ptxas] = find_cuda_binary("ptxas", toolkit_dirs)

    config[:configured] = true

    ## (re)generate ext.jl

    if isfile(previous_config_path)
        previous_config = read_ext(previous_config_path)

        if config == previous_config
            mv(previous_config_path, config_path; force=true)
            return
        end
    end

    write_ext(config, config_path)

    return
end

main()


text:
module commands

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using Unidecode, Test

@test unidecode("α") == "alpha"
@test unidecode("🍫") == ":chocolate_bar:"


text:
module SimpleReplicable

using Test
using Hyperspecialize

global A = Set{Type}([])

f(::Real) = false
g(::Signed, ::Unsigned) = false
h(::Signed, ::Signed, ::Unsigned) = false

@testset "Simple Replicable" begin

  # First, checking if the hyperspecialization works like concretization
  empty!(A)

  @replicable begin
    push!(A, @hyperspecialize Real)
  end

  @test A == Set{Type}(Hyperspecialize.concretesubtypes(Real))

  empty!(A)

  @concretize NotAType []

  @replicable begin
    push!(A, @hyperspecialize NotAType)
  end

  @test A == Set{Type}()

  @test_throws ErrorException (@concretize NotAType Int8)

  @widen NotAType Int8

  @test A == Set{Type}([Int8])

  empty!(A)

  # Can we widen and then get replicables to follow?

  @replicable function f(::@hyperspecialize(NotAType))
    return true
  end

  @test !f(1.0)
  @test f(Int8(1))

  @widen NotAType Float64

  @test f(1.0)

  # What happens if a replicable has no hyperspecializes?

  @replicable function f(::String)
    return true
  end

  @test f("hello")

  # How does specialization interact with multiple hyperspecials?

  @concretize Signed (Int8, Int16)
  @concretize Unsigned []

  @replicable function g(::@hyperspecialize(Signed), ::@hyperspecialize(Unsigned))
    return true
  end

  @test !g(Int8(1), UInt8(1))
  @test !g(Int16(1), UInt8(1))
  @test !g(Int8(1), UInt16(1))
  @test !g(Int16(1), UInt16(1))

  @widen Unsigned [UInt8, UInt16]

  @test g(Int8(1), UInt8(1))
  @test g(Int16(1), UInt8(1))
  @test g(Int8(1), UInt16(1))
  @test g(Int16(1), UInt16(1))

  @test !h(Int8(1), Int8(1), UInt8(1))
  @test !h(Int16(1), Int8(1), UInt8(1))
  @test !h(Int8(1), Int16(1), UInt8(1))
  @test !h(Int16(1), Int16(1), UInt8(1))
  @test !h(Int8(1), Int8(1), UInt16(1))
  @test !h(Int16(1), Int8(1), UInt16(1))
  @test !h(Int8(1), Int16(1), UInt16(1))
  @test !h(Int16(1), Int16(1), UInt16(1))

  @replicable function h(::@hyperspecialize(Signed), ::@hyperspecialize(Signed), ::@hyperspecialize(Unsigned))
    return true
  end

  @test h(Int8(1), Int8(1), UInt8(1))
  @test h(Int16(1), Int8(1), UInt8(1))
  @test h(Int8(1), Int16(1), UInt8(1))
  @test h(Int16(1), Int16(1), UInt8(1))
  @test h(Int8(1), Int8(1), UInt16(1))
  @test h(Int16(1), Int8(1), UInt16(1))
  @test h(Int8(1), Int16(1), UInt16(1))
  @test h(Int16(1), Int16(1), UInt16(1))

  @test !h(Int32(1), Int16(1), UInt16(1))
  @test !h(Int16(1), Int16(1), UInt64(1))

end

end #module


text:
"""
Attach `k` pre-allocated children to node `p` and repeat on its descendants up to `t` times.
`k` is increased by Δk at every forking time.
"""
function kfurcate_symmetric!(p::AbstractNode, k, i, t, tree, Δk=0)
    if t > 0
        for _ in 1:k
            i += 1
            c = tree.nodes[i]
            setlabel!(c, string(i))
            link!(p => c)
            i = kfurcate_symmetric!(c, k + Δk, i, t - 1, tree, Δk)
        end
    end

    return i
end

"""
    symmetric_tree([Type{<:AbstractTree}], n, k = 2)
    symmetric_tree([Type{<:AbstractTree}]; depth, k = 2)

Create a symmetric tree with forks of `k` degree (`k` = 2 for a dichotomic tree).

The size of the tree can be set either by the number `n` of desired outer nodes, or by the
desired `depth` d of the tree (number of nodes from the root to an outer node). The relation
between the three parameters is

    k^d = n
    k, d, n ∈ ℕ

# Examples
```jldoctest
julia>symmetric_tree(Any, 9, k = 3)
RTree{Any}: 13 nodes
```
"""
function symmetric_tree(::Type{ND}, ::Type{BD}; depth, k=2) where {ND,BD}
    N = sum(k .^ (depth:-1:1)) + 1
    tree = createtree(ND, BD, N)
    tree.root = tree.nodes[1]
    setlabel!(tree.root, "1")
    kfurcate_symmetric!(getroot(tree), k, 1, depth, tree, 0)

    return tree
end

function symmetric_tree(::Type{ND}, ::Type{BD}, n; k = 2) where {ND,BD}
    d = log(k, n)    ## Number of times to fork terminal nodes
    isinteger(d) ||
        @error "`n` is not a natural number power of `k`"

    return symmetric_tree(ND, BD; depth = Int(d), k = k)
end

symmetric_tree(n; k=2) = symmetric_tree(Dict, Dict, n; k = k)

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
using GeoPhyInv
using Test
using SparseArrays

mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];

ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(10));

ageom=AGeom(mgrid, :xwell, SSrcs(10), Recs(10));

@test (ageom ∈ mgrid)

update!(ageom[1], Srcs(), [0,1], [10,20],);
update!(ageom[1], Recs(), [0,0], [10,20],);
update!(ageom, SSrcs(), [0,1], [10,20], );

ageom2=vcat(ageom, ageom);


text:
module GridapTypesTests

using Test
using Gridap.Helpers

import Gridap.Helpers: operate

struct Foo{A} <: GridapType
  bar::A
end

object = Foo(nothing)

@test sprint(show,"text/plain",object) == "$(nameof(typeof(object)))()"

function operate(op,a::Foo)
  Foo(op(a.bar))
end

function operate(op,a::Foo,b::Foo)
  Foo(op(a.bar,b.bar))
end

function operate(op,a::Number,b::Foo)
  operate(op,Foo(a),b)
end

function operate(op,a::Foo,b::Number)
  operate(op,a,Foo(b))
end

a = Foo(1)

b = -a
@test b.bar == -1

c = a + b

@test c.bar == 0

c = a - 2
@test c.bar == -1

c = 2 - a
@test c.bar == 1

using Gridap.TensorValues

c = inner(2,a)
@test c.bar == 2

c = a / b
@test c.bar == -1

end # module


text:
tuple_dot(a, b) = sum(a .+ b)

"""
function returning true at nozzle position
"""
function isjet(p, nozzle_cen, nozzle_len, nozzle_rad)
    (abs(p[1] - nozzle_cen[1]) <= nozzle_len / 2) &
    ((p[2] - nozzle_cen[2])^2 +
    (p[3] - nozzle_cen[3])^2 <= nozzle_rad .^ 2) != 0
end

function restrict_kernel(psi, isjet, kvec, pos, omgterm)
    if isjet
        amp = @fastmath abs.(psi)
        phase = tuple_dot(kvec, pos) - omgterm
        @fastmath amp .* exp(Complex64(0f0, 1f0) .* phase)
    else
        psi
    end
end

"""
Restrict the velocity in restrict_positions
"""
function restrict_velocity!(isf, psi, kvec, restrict_positions, omgterm = 0f0)
    psi .= restrict_kernel.(psi, restrict_positions, (kvec,), isf.positions, omgterm)
end

"""
Add particles around a jet on the CPU
"""
function add_particles!(particles::Vector, range, nozzle_cen, nozzle_rad)
    randstate = GPUArrays.cached_state(particles)
    for i in range
        rt = rand(Float32) * 2f0 * pi

        spawn_offset = ((
            rand(Float32),
            rand(Float32),
            rand(Float32),
        ) .* 0.1f0) .- 0.05f0

        spawned = ((
            nozzle_cen[1] - 0.1f0,
            nozzle_cen[2] + 0.9f0 * nozzle_rad * cos(rt),
            nozzle_cen[3] + 0.9f0 * nozzle_rad * sin(rt)
        )) .+ spawn_offset

        particles[i] = spawned
    end
    return
end

"""
Add particles around a jet on the GPU!
"""
function add_particles!(particles, range, nozzle_cen, nozzle_rad)
    randstate = GPUArrays.cached_state(particles)
    N = length(range)
    gpu_call(
        add_particles!, particles,
        (
            randstate, particles, UInt32(first(range)), UInt32(N),
            nozzle_cen, nozzle_rad, Float32(pi)
        ),
        N
    )
    return
end

function add_particles!(state, randstate, particles, offset, n, nozzle_cen, nozzle_rad, pi)
    ilin = linear_index(state)
    ilin > n && return
    rt = gpu_rand(Float32, state, randstate) * 2f0 * pi

    spawn_offset = ((
        gpu_rand(Float32, state, randstate),
        gpu_rand(Float32, state, randstate),
        gpu_rand(Float32, state, randstate),
    ) .* 0.1f0) .- 0.05f0

    spawned = ((
        nozzle_cen[1] - 0.1f0,
        nozzle_cen[2] + 0.9f0 * nozzle_rad * cos(rt),
        nozzle_cen[3] + 0.9f0 * nozzle_rad * sin(rt)
    )) .+ spawn_offset

    particles[ilin + offset] = spawned
    return
end

export add_particles!, restrict_velocity!, isjet


text:
using DesignModels
using Test

import DesignModels.Design, DesignModels.Facet

@testset "DesignModels.jl" begin
    design = Design()
    facet = Facet()
    # println(convert_design("../examples/ascii.stl"; compact=false))
    for f in [
        "ascii",
        "primitive closed box",
        "primitive floor",
        "primitive open box"
        ]
        convert_design("../examples/$f.stl"; file = "../examples/$f.json")
    end

end


text:
using SparseArrays

struct optimization_parameters
    equation
    optim
    eq_params::Array{<:Real, 1}
    k_sweeps::Integer
    k_batch::Integer
    epochs::Integer
    BCs::Array{String, 1}
    L::Integer
    name::String
    symmetry::String
    error::String
end

struct preconditioner
    A::SparseMatrixCSC{<:Real, <:Integer}
    Ind::Array{Any,1}
    BCs::Array{String, 1}
    L::Integer
    name::String
    symmetry::String
end

struct T_preconditioner
    A::SparseMatrixCSC{<:Real, <:Integer}
    C::Array{Array{SparseMatrixCSC{Float64,Int64},1},1}
    Ind::Array{Any,1}
    BCs::Array{String, 1}
    L::Integer
    name::String
    symmetry::String
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
# NOTE: You pass in weights using a row vector, e.g. [500 300 100], then I convert it to a 1-dim column vector
function alex_cifar_tutorial_nn(count_classes::Int)
    model_name = "alex_cifar_tutorial"

    # Returns a function that returns a function to create model
    nn_factory = function()
      conv1_layer = ConvolutionLayer(name="conv1", n_filter=32, kernel=(5,5), pad=(2,2), stride=(1,1),
                                       filter_regu = NoRegu(),
                                       #filter_init=GaussianInitializer(std=0.0001),
                                       bottoms=[:data], tops=[:conv1])
      pool1_layer = PoolingLayer(name="pool1", kernel=(3,3), stride=(2,2), neuron=Neurons.ReLU(),
                                       bottoms=[:conv1], tops=[:pool1])
      norm1_layer = LRNLayer(name="norm1", kernel=3, scale=5e-5, power=0.75, mode=LRNMode.WithinChannel(),
                                       bottoms=[:pool1], tops=[:norm1])

      conv2_layer = ConvolutionLayer(name="conv2", n_filter=32, kernel=(5,5), pad=(2,2), stride=(1,1),
                                       filter_regu = NoRegu(),
                                       #filter_init=GaussianInitializer(std=0.01),
                                       bottoms=[:norm1], tops=[:conv2], neuron=Neurons.ReLU())
      pool2_layer = PoolingLayer(name="pool2", kernel=(3,3), stride=(2,2), pooling=Pooling.Mean(),
                                       bottoms=[:conv2], tops=[:pool2])
      norm2_layer = LRNLayer(name="norm2", kernel=3, scale=5e-5, power=0.75, mode=LRNMode.WithinChannel(),
                                       bottoms=[:pool2], tops=[:norm2])

      conv3_layer = ConvolutionLayer(name="conv3", n_filter=64, kernel=(5,5), pad=(2,2), stride=(1,1),
                                       filter_regu = NoRegu(),
                                       #filter_init=GaussianInitializer(std=0.01),
                                       bottoms=[:norm2], tops=[:conv3], neuron=Neurons.ReLU())
      pool3_layer = PoolingLayer(name="pool3", kernel=(3,3), stride=(2,2), pooling=Pooling.Mean(),
                                       bottoms=[:conv3], tops=[:pool3])

      ip1_layer   = InnerProductLayer(name="ip1", output_dim=count_classes,
                                       #weight_init=GaussianInitializer(std=0.01),
                                       weight_regu=NoRegu(),
                                       bottoms=[:pool3], tops=[:ip1])

       pred_prob_layer = SoftmaxLayer( name = "pred_probs", bottoms = [:ip1],tops= [:pred_probs])
       storage_layer = MemoryOutputLayer( name = "storage", bottoms = [:pred_probs])
      loss_layer  = SoftmaxLossLayer(name="loss", bottoms=[:ip1, :label])
      acc_layer   = AccuracyLayer(name="test-accuracy", bottoms=[:ip1, :label])

      common_layers = [conv1_layer, pool1_layer, norm1_layer, conv2_layer, pool2_layer, norm2_layer,
                                       conv3_layer, pool3_layer, ip1_layer]

      return (common_layers, pred_prob_layer, storage_layer, loss_layer, acc_layer)
    end

    return (nn_factory, model_name)
end


text:
cd(@__DIR__); include("setups/rgc100.jl")
gr(dpi = 200)

## construct full HGLET and Lapped-HGLET dictionaries
@time HGLET_dic = HGLET_dictionary(GP, G_Sig; method = :L)
@time LPHGLET_dic = LPHGLET_dictionary(GP, G_Sig; method = :L, ϵ = 0.3)

## generate figures
# pre-selected (j, k, l)s
T = [(2, 1, 1), (2, 1, 12), (3, 1, 3), (3, 1, 5),
    (4, 5, 6), (6, 23, 3), (7, 5, 2), (8, 15, 1)]

for (j, k, l) in T
    𝚽H = HGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
    𝚽lH = LPHGLET_dic[GP.rs[k, j]:(GP.rs[k + 1, j] - 1), j, :]'
    p1 = gplot(W, X; width = 1)
        scatter_gplot!(X; marker = 𝚽H[:, l], plotOrder = :s2l)
        plot!(xlims = [-180, 220], xlabel = "x(μm)", ylabel = "y(μm)", frame = :box, cbar = false)
        plot!(left_margin = 3mm, bottom_margin = 3mm, clims = (-0.05, 0.05))
    p2 = gplot(W, X; width = 1)
        scatter_gplot!(X; marker = 𝚽lH[:, l], plotOrder = :s2l)
        plot!(xlims = [-180, 220], xlabel = "x(μm)", ylabel = "y(μm)", frame = :box, cbar = false)
        plot!(left_margin = 3mm, bottom_margin = 3mm, clims = (-0.05, 0.05))
    plt = plot(p1, p2, size = (800, 400), layout = Plots.grid(1, 2))
    savefig(plt, "../figs/RGC100_HGLET_LPHGLET_j$(j-1)k$(k-1)l$(l-1).png")
end


text:
function test_binary_crossentropy_loss_layer(backend::Backend, tensor_dim, T, eps)
  println("-- Testing BinaryCrossEntropyLossLayer on $(typeof(backend)){$T}...")

  dims = abs(rand(Int,tensor_dim)) % 6 + 2

  println("    > $dims")

  #dims_label = copy(dims); dims_label[op_dim] = 1; dims_label = tuple(dims_label...)
  dims = tuple(dims...)

  prob = rand(T, dims)

  label = rand(Int, dims) .> 0.5
  label = convert(Array{T}, label)

  prob_blob = make_blob(backend, prob)
  label_blob = make_blob(backend, label)
  inputs = Blob[prob_blob, label_blob]
  weight = 0.25
  layer = BinaryCrossEntropyLossLayer(bottoms=[:pred, :labels], weight=weight)
  state = setup(backend, layer, inputs, Blob[])

  forward(backend, state, inputs)

  expected_loss = convert(T, 0)

  for i = 1:prod(dims)
      expected_loss += -log(vec(prob)[i])*vec(label)[i]
      expected_loss += -log(1 - vec(prob)[i])*vec(1 - label)[i]
  end

  expected_loss /= dims[end]
  expected_loss *= weight

  @test -eps < 1 - state.loss/expected_loss < eps

  diff_blob2  = make_blob(backend, prob)
  diff_blob1  = make_blob(backend, prob)
  diffs = Blob[diff_blob1, diff_blob2]
  backward(backend, state, inputs, diffs)
  grad_pred = -weight * (label./prob - (1-label)./(1-prob) ) / dims[end]
  diff = similar(grad_pred)

  copy!(diff, diffs[1])

  @test all(-eps .< 1 - grad_pred./diff .< eps)

  grad_label = -weight * log(prob./(1-prob)) / dims[end]
  diff = similar(grad_pred)
  copy!(diff, diffs[2])

  @test all(-eps .< grad_label - diff .< eps)

  shutdown(backend, state)
end

function test_binary_crossentropy_loss_layer(backend::Backend, T, eps)
  for i in [2,4,5]
      test_binary_crossentropy_loss_layer(backend, i, T, eps)
  end
end

function test_binary_crossentropy_loss_layer(backend::Backend)
  test_binary_crossentropy_loss_layer(backend, Float32, 1e-5)
  test_binary_crossentropy_loss_layer(backend, Float64, 1e-6)
end

if test_cpu
  test_binary_crossentropy_loss_layer(backend_cpu)
end
if test_gpu
  test_binary_crossentropy_loss_layer(backend_gpu)
end