text:
import MitosisStochasticDiffEq as MSDE
using Mitosis
using StochasticDiffEq
using Test, Random
using Statistics
using LinearAlgebra

"""
    conjugate_posterior(Y, Ξ)

Sample the posterior distribution of the conjugate drift parameters from path `Y`,
prior precision matrix `Ξ` with non-conjugate parameters fixed in the model.
Adjusted from http://www.math.chalmers.se/~smoritz/journal/2018/01/19/parameter-inference-for-a-simple-sir-model/
for Bridge.jl
"""
function conjugate_posterior(ts, Ys, Ξ, p, g)
    paramgrad(t, u) = [u, 1]
    paramintercept(t, u) = 0
    t, y = ts[1], Ys[1]
    ϕ = paramgrad(t, y)
    mu = zero(ϕ)
    G = zero(mu*mu')

    for i in 1:length(Ys)-1
        ϕ = paramgrad(t, y)'
        Gϕ = pinv(g(y, p, t)*g(y, p, t)')*ϕ # a is sigma*sigma'. Todo: smoothing like this is very slow
        zi = ϕ'*Gϕ
        t2, y2 = ts[i + 1], Ys[i + 1]
        dy = y2 - y
        ds = t2 - t
        #@show size(mu), size(Gϕ'), (dy - paramintercept(t, y)*ds)
        mu = mu + Gϕ'*(dy - paramintercept(t, y)*ds)
        t, y = t2, y2
        G = G + zi*ds
    end
    Mitosis.Gaussian{(:F,:Γ)}(mu, G + Ξ)#, mulist
end

@testset "regression tests" begin

  K = 1000

  # define SDE function
  f(du,u,p,t) = du .= p[1]*u .+ p[2]
  g(du,u,p,t) = du .= p[3]
  foop(u,p,t) = p[1]*u .+ p[2]
  goop(u,p,t) = p[3]

  # paramjac
  function f_jac(J,u,p,t)
    J[1,1] = u[1]
    J[1,2] = true
    nothing
  end

  function f_jacoop(u,p,t)
    [u[1]  true]
  end

  # intercept
  function ϕ0(du,u,p,t)
    du .= false
  end

  function ϕ0oop(u,p,t)
    [zero(eltype(u))]
  end

  ϕprototype = zeros((1,2))
  yprototype = zeros(1)

  # time span
  tstart = 0.0
  tend = 100.0
  dt = 0.001
  trange = tstart:dt:tend

  # initial condition
  u0 = 1.1

  # set true model parameters
  par = [-0.3, 0.2, 0.5]

  # define SDE kernel
  sdekernel = MSDE.SDEKernel(foop,goop,trange,par)
  sdekernel2 = MSDE.SDEKernel(f,g,trange,par)

  # sample using MSDE and EM default
  Random.seed!(100)
  uend, (ts, u, noise) = MSDE.sample(sdekernel, u0, save_noise=true)
  Random.seed!(100)
  # for later AD check
  uend2, (ts2, u2, noise2) = MSDE.sample(sdekernel, [u0], save_noise=true)
  @show uend
  @show length(u)

  R = MSDE.Regression!(sdekernel,yprototype,
      paramjac_prototype=ϕprototype,paramjac=f_jac,intercept=ϕ0)
  R2 = MSDE.Regression(sdekernel,paramjac=f_jacoop,intercept=ϕ0oop)
  R3 = MSDE.Regression!(sdekernel2,yprototype,
      paramjac_prototype=ϕprototype,paramjac=f_jac,intercept=ϕ0)
  R4 = MSDE.Regression(sdekernel2,
      dyprototype=yprototype,paramjac=f_jacoop,intercept=ϕ0oop)

  G = MSDE.conjugate(R, u, 0.1*I(2), ts)
  G2 = MSDE.conjugate(R2, u, 0.1*I(2), ts)
  G3 = MSDE.conjugate(R3, u, 0.1*I(2), ts)
  G4 = MSDE.conjugate(R4, u, 0.1*I(2), ts)
  Gtest = conjugate_posterior(ts, u, 0.1*I(2), par, goop)

  @testset "Regression! oop tests" begin
    @test G ≈ Gtest rtol=1e-10
    @test G.F ≈ Gtest.F rtol=1e-10
    @test G.Γ ≈ Gtest.Γ rtol=1e-10
  end
  @testset "Regression oop tests" begin
    @test G2 ≈ Gtest rtol=1e-10
    @test G2.F ≈ Gtest.F rtol=1e-10
    @test G2.Γ ≈ Gtest.Γ rtol=1e-10
  end
  @testset "Regression! iip tests" begin
    @test G3 ≈ Gtest rtol=1e-10
    @test G3.F ≈ Gtest.F rtol=1e-10
    @test G3.Γ ≈ Gtest.Γ rtol=1e-10
  end
  @testset "Regression iip tests" begin
    @test G4 ≈ Gtest rtol=1e-10
    @test G4.F ≈ Gtest.F rtol=1e-10
    @test G4.Γ ≈ Gtest.Γ rtol=1e-10
  end
  @info Gtest.F
  @info Gtest.Γ

  # test samples
  mu = G.F
  Gamma = G.Γ
  WL = (cholesky(Hermitian(Gamma)).U)'

  Random.seed!(1)
  Π = []
  Π2 = []
  for i=1:K
    th° = WL'\(randn(size(mu))+WL\mu)
    push!(Π,th°)
    th° = WL'\(randn(size(mu))+WL\mu)
    push!(Π2,th°)
  end

  mu = G2.F
  Gamma = G2.Γ
  WL = (cholesky(Hermitian(Gamma)).U)'

  Random.seed!(1)
  for i=1:K
    th° = WL'\(randn(size(mu))+WL\mu)
    push!(Π2,th°)
  end
  @testset "regression samples tests" begin
    @test par[1:2] ≈ mean(Π) rtol=0.2
    @test par[1:2] ≈ mean(Π2) rtol=0.2
    @test mean(Π) ≈ mean(Π2) atol=0.1
  end
# using Plots
# pl = scatter(first.(Π), last.(Π), markersize=1, c=:blue, label="posterior samples")
# scatter!(first.(Π2), last.(Π2), markersize=1, c=:green, label="posterior samples")
# scatter!([par[1]], [par[2]], color="red", label="truth")
# savefig(pl, "regression.png")

  # test with ForwardDiff
  @testset "AD for Jacobian tests" begin
    using ForwardDiff
    pf = MSDE.ParamJacobianWrapper(sdekernel.f,first(sdekernel.trange),[u0])
    pf2 = MSDE.ParamJacobianWrapper2(sdekernel2.f,first(sdekernel2.trange),[u0])
    f_jac(ϕprototype,[u0],par,first(sdekernel.trange))
    @test f_jacoop(u0,par,first(sdekernel.trange)) == ϕprototype
    @test ForwardDiff.jacobian(pf, par[1:2]) == ϕprototype
    @test ForwardDiff.jacobian(pf2, par[1:2]) == ϕprototype

    # check θ function
    RAD = MSDE.Regression!(sdekernel,yprototype,
      paramjac_prototype=ϕprototype,intercept=ϕ0,θ=par[1:2])
    GAD = MSDE.conjugate(RAD, u, 0.1*I(2), ts)
    @test G ≈ GAD rtol=1e-10

    #check with intercept = nothing
    RAD = MSDE.Regression!(sdekernel,yprototype,
      paramjac_prototype=ϕprototype,θ=par[1:2])
    RAD2 = MSDE.Regression(sdekernel,θ=par[1:2],yprototype=yprototype)
    GAD = MSDE.conjugate(RAD, u, 0.1*I(2), ts)
    GAD2 = MSDE.conjugate(RAD2, u2, 0.1*I(2), ts2)
    @test G ≈ GAD rtol=1e-10
    @test G ≈ GAD2 rtol=1e-10

    #check AD for inplace function
    RAD3 = MSDE.Regression!(sdekernel2,yprototype,
        paramjac_prototype=ϕprototype,θ=par[1:2])
    RAD4 = MSDE.Regression(sdekernel2,
        dyprototype=yprototype,yprototype=yprototype,θ=par[1:2])
    GAD3 = MSDE.conjugate(RAD3, u, 0.1*I(2), ts)
    GAD4 = MSDE.conjugate(RAD4, u2, 0.1*I(2), ts2)

    @test G ≈ GAD3 rtol=1e-10
    @test G ≈ GAD4 rtol=1e-10
  end
end

@testset "regression SDE Problem tests" begin
  # define SDE function
  foop(u,p,t) = p[1]*u .+ p[2]
  goop(u,p,t) = p[3]

  # paramjac
  function f_jac(J,u,p,t)
    J[1,1] = u[1]
    J[1,2] = true
    nothing
  end

  function f_jacoop(u,p,t)
    [u[1]  true]
  end

  # intercept
  function ϕ0(du,u,p,t)
    du .= false
  end

  function ϕ0oop(u,p,t)
    [zero(eltype(u))]
  end

  ϕprototype = zeros((1,2))
  yprototype = zeros(1)

  # time span
  tstart = 0.0
  tend = 100.0
  dt = 0.01
  trange = tstart:dt:tend

  # initial condition
  u0 = 1.1

  # set true model parameters
  par = [-0.3, 0.2, 0.5]

  # define SDE Problem and sample using EM
  Random.seed!(100)
  prob = SDEProblem(foop, goop, u0, (tstart,tend), par)
  sol = solve(prob, EM(), dt = dt)

  @show sol[end]

  R = MSDE.Regression!(prob,yprototype,
     paramjac_prototype=ϕprototype,paramjac=f_jac,intercept=ϕ0)
  R2 = MSDE.Regression(prob,paramjac=f_jacoop,intercept=ϕ0oop)

  G = MSDE.conjugate(R, sol, 0.1*I(2))
  G2 = MSDE.conjugate(R2, sol, 0.1*I(2))
  G3 = conjugate_posterior(sol.t, sol.u, 0.1*I(2), par, goop)

  @testset "iip tests" begin
    @test G ≈ G3 rtol=1e-10
    @test G.F ≈ G3.F rtol=1e-10
    @test G.Γ ≈ G3.Γ rtol=1e-10
  end
  @testset "oop tests" begin
    @test G2 ≈ G3 rtol=1e-10
    @test G2.F ≈ G3.F rtol=1e-10
    @test G2.Γ ≈ G3.Γ rtol=1e-10
  end
  @info G3.F
  @info G3.Γ
end

@testset "scalar Brusselator" begin
  # define SDE function for Brusselator, test inplace version and AD
  using DiffEqNoiseProcess, Zygote, DiffEqSensitivity
  function brusselator_f!(du,u,p,t)
    @inbounds begin
      du[1] = (p[1]-1)*u[1]+p[1]*u[1]^2+(u[1]+1)^2*u[2]
      du[2] = -p[1]*u[1]-p[1]*u[1]^2-(u[1]+1)^2*u[2]
    end
    nothing
  end

  function brusselator_f(u,p,t)
    dx1 = (p[1]-1)*u[1]+p[1]*u[1]^2+(u[1]+1)^2*u[2]
    dx2 = -p[1]*u[1]-p[1]*u[1]^2-(u[1]+1)^2*u[2]

    return [dx1,dx2]
  end

  function scalar_noise!(du,u,p,t)
    @inbounds begin
      du[1] = p[2]*u[1]*(1+u[1])
      du[2] = -p[2]*u[1]*(1+u[1])
     end
     nothing
  end

  function scalar_noise(u,p,t)
    dx1 = p[2]*u[1]*(1+u[1])
    dx2 = -p[2]*u[1]*(1+u[1])

    return [dx1,dx2]
  end

  # matrix form
  function scalar_noise_matrix(u,p,t)
    dx1 = p[2]*u[1]*(1+u[1])
    dx2 = -p[2]*u[1]*(1+u[1])

    reshape([ dx1, dx2],2,1)
  end

  # paramjac
  function f_jac!(J,u,p,t)
    J[1,1] = u[1]+u[1]^2
    J[2,1] = -(u[1]+u[1]^2)
    nothing
  end

  function f_jac(u,p,t)
    [ u[1]+u[1]^2
      -(u[1]+u[1]^2)]
  end

  # intercept
  function ϕ0!(du,u,p,t)
    du[1] = -u[1]+(u[1]+1)^2*u[2]
    du[2] = -(u[1]+1)^2*u[2]
  end

  function ϕ0(u,p,t)
    dx1 = -u[1]+(u[1]+1)^2*u[2]
    dx2 = -(u[1]+1)^2*u[2]
    [dx1, dx2]
  end

  yprototype = zeros(2)
  ϕprototype = zeros((2,1))

  # fix seeds
  seed = 100
  Random.seed!(seed)
  W = WienerProcess(0.0,0.0,nothing)

  u0 = [-0.1,0.0]
  tspan = (0.0,100.0)
  p = [1.9,0.1] # p[1] is only in the drift

  prob = SDEProblem(brusselator_f!,scalar_noise!,u0,tspan,p,noise=W)
  sol = solve(prob, EM(), dt = 0.01)

  @test size(scalar_noise_matrix(sol.u[1], p, 0.0)) == (2,1)

  testvec = randn(2)
  o1 = MSDE.outer_(scalar_noise_matrix(testvec, p, 0.0))
  o2 = MSDE.outer_(reshape(scalar_noise(testvec,p,0.0),2,1))
  @test o1==o2

  # oop regression with matrix form for diffusion to test
  sdekernel1 = MSDE.SDEKernel(brusselator_f,scalar_noise_matrix,tspan,p)
  R1 = MSDE.Regression(sdekernel1,paramjac=f_jac,intercept=ϕ0)
  # oop regression with vector form as in StochasticDiffEq
  sdekernel2 = MSDE.SDEKernel(brusselator_f,scalar_noise,tspan,p)
  R2 = MSDE.Regression(sdekernel2,paramjac=f_jac,intercept=ϕ0, m=1)
  # inplace regression with vector form as in StochasticDiffEq
  sdekernel3 = MSDE.SDEKernel(brusselator_f!,scalar_noise!,tspan,p)
  R3 = MSDE.Regression!(sdekernel3,yprototype,
     paramjac_prototype=ϕprototype,paramjac=f_jac!,intercept=ϕ0!,m=1)

  G1 = MSDE.conjugate(R1, sol, 0.1*I(1))
  G2 = MSDE.conjugate(R2, sol, 0.1*I(1))
  G3 = MSDE.conjugate(R3, sol, 0.1*I(1))

  @test G1.F ≈ G2.F rtol=1e-10
  @test G1.Γ ≈ G2.Γ rtol=1e-10
  @test G1.F ≈ G3.F rtol=1e-10
  @test G1.Γ ≈ G3.Γ rtol=1e-10
  @test G2.F ≈ G3.F rtol=1e-10
  @test G2.Γ ≈ G3.Γ rtol=1e-10

  # Estimate
  p̂ = mean(G1)
  se = sqrt.(diag(cov(G1)))
  display(map((p̂, se, p) -> "$(round(p̂, digits=3)) ± $(round(se, digits=3)) (true: $p)", p̂, se, p))
  @test p̂[1] ≈ p[1] rtol=se[1]

  # use NoiseWrapper to reproduce same trajectory
  # Wrep = NoiseWrapper(sol.W)
  ts = sol.t
  gaussian = Mitosis.Gaussian{(:F,:Γ)}(zeros(eltype(Matrix(0.1*I(1))), size(Matrix(0.1*I(1)), 1)), Matrix(0.1*I(1)))

  function loss(p,u0,R)
    _prob = remake(prob,p = p, u0 = u0, noise=W)
    _sol = (solve(_prob, EM(), dt = 0.01)).u
    G = MSDE.conjugate(R, _sol, gaussian, ts)
    p̂ = mean(G)[1]
    se = sqrt.(diag(cov(G)))[1]
    se + (p̂-p[1])^2
  end
  @test loss(p,u0,R1) != zero(loss(p,u0,R1))
  @test loss(p,u0,R2) != zero(loss(p,u0,R2))
  @test loss(p,u0,R3) != zero(loss(p,u0,R3))

  @test loss(p,u0,R1) ≈ loss(p,u0,R2) atol=1e-3
  @test loss(p,u0,R2) ≈ loss(p,u0,R3) atol=1e-3

  du01,dp1 = Zygote.gradient((u0,p)->loss(p,u0,R1),u0,p)
  @test du01 != zero(du01)
  @test dp1 != zero(dp1)
end


text:
# # Binary (or 0-1) knapsack problem
# Given a knapsack of some capacity $C$ and $n$ objects with object $i$ having weight $w_i$ and profit $p_i$, the goal is to choose some subset of the objects that can fit in the knapsack (i.e. the sum of their weights is no more than $C$) while maximizing profit.
#
# This can be formulated as a mixed-integer program as:
#
# $$
# \begin{array}{ll}
#   \text{maximize} & x' p \\
#     \text{subject to} & x \in \{0, 1\} \\
#   & w' x \leq C \\
# \end{array}
# $$
#
# where $x$ is a vector is size $n$ where $x_i$ is one if we chose to keep the object in the knapsack, 0 otherwise.

## Data taken from http://people.sc.fsu.edu/~jburkardt/datasets/knapsack_01/knapsack_01.html
w = [23; 31; 29; 44; 53; 38; 63; 85; 89; 82]
C = 165
p =  [92; 57; 49; 68; 60; 43; 67; 84; 87; 72];
n = length(w)

#-

using Convex, GLPK
x = Variable(n, :Bin)
problem = maximize(dot(p, x), dot(w, x) <= C)
solve!(problem, GLPK.Optimizer)
evaluate(x)

ext: jl


text:
module CuSynchronization

using CUDAnative.LLVM
using CUDAnative.LLVM.Interop
using CUDAnative

@generated function unsafe_volatile_load(ptr::Ptr{T}) where T
    eltyp = convert(LLVMType, T)
    T_ptr = convert(LLVMType, Ptr{T})
    T_actual_ptr = LLVM.PointerType(eltyp)

    # create a function
    param_types = [T_ptr]
    llvm_f, _ = create_function(eltyp, param_types)

    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ld = load!(builder, ptr)
        LLVM.API.LLVMSetVolatile(LLVM.ref(ld), LLVM.True)
        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{Ptr{T}}, :(ptr,))
end

@generated function unsafe_volatile_store!(ptr::Ptr{T}, val::T) where T
    eltyp = convert(LLVMType, T)
    T_ptr = convert(LLVMType, Ptr{T})
    T_actual_ptr = LLVM.PointerType(eltyp)

    # create a function
    param_types = [T_ptr, eltyp]
    llvm_f, _ = create_function(LLVM.VoidType(JuliaContext()), param_types)

    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        val = parameters(llvm_f)[2]
        st = store!(builder, val, ptr)
        LLVM.API.LLVMSetVolatile(LLVM.ref(st), LLVM.True)
        ret!(builder)
    end

    call_function(llvm_f, Cvoid, Tuple{Ptr{T}, T}, :(ptr, val))
end

import CUDAnative: DevicePtr

struct Semaphore{T, AS}
    sem::DevicePtr{T, AS}
    value::T
end

@inline function Base.getindex(sem::Semaphore{T}) where T
    ptr = Base.unsafe_convert(Ptr{T}, sem.sem)
    unsafe_volatile_load(ptr)
end

@inline function Base.setindex!(sem::Semaphore{T}, val) where T
    ptr = Base.unsafe_convert(Ptr{T}, sem.sem)
    unsafe_volatile_store!(ptr, convert(T, val))
end

@inline function wait(sem::Semaphore{T}) where T
    while true
        if sem[] == sem.value
            break
        end
        threadfence_block()
    end
    sem[] = 2 % T # finalize
    threadfence_system()
    return nothing
end

end # CuSynchronization


text:
import PowerModelsACDC;
const _PMACDC = PowerModelsACDC;
import PowerModels;
const _PM = PowerModels;
import InfrastructureModels;
const _IM = InfrastructureModels;
import JuMP
import Gurobi
using MAT
using XLSX
using JLD2
using Statistics
include("basencont_nw.jl")

Total_sample = 500  # sample per year
total_yr = 3# the years in horizon, data coming from excels
period = "multi" # single or multi

Prot_system = "Permanentloss"
curtailed_gen = [1] #geneartor numbers # change also constraint max(), generating power,file name
syncarea = 2
max_curt = 0
@load "scenario_500.jld2"

    file = "./test/data/4bus_OPF_PLdim.m"
    data_sp = _PM.parse_file(file)
    _PMACDC.process_additional_data!(data_sp)

    include("bkgrnd_cal.jl")
    gurobi = JuMP.optimizer_with_attributes(Gurobi.Optimizer, "Presolve" => -1)

    no_nw = 2 * Total_sample*total_yr
    data_mp = multi_network(file, no_nw)
    year = ["NAT_2025_Generation","NAT_2035_Generation", "NAT_2045_Generation"]
    year_num = ["2025", "2035", "2045"]
    # year = ["NAT_2050_Generation"]
    yr = 1

    # @load "cluster.jld2"

    data_ps = mp_datainputs_nocl_PLdim(deepcopy(data_mp), Total_sample, year, year_num, scenario_500, file)
    data_cont, Cont_list, base_list = mp_contignecy_nocl(deepcopy(data_ps), Total_sample*total_yr, 1)
    year_base = year_base_networks(Total_sample, total_yr, base_list)

    if Prot_system == "FS_HDCCB" || Prot_system == "FS_MDCCB"
    s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => true, "NSprotection" => false,
    "Permanentloss" => false, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea, "year_base" => year_base, "total_yr" => total_yr)
    elseif Prot_system == "NS_CB" || Prot_system == "NS_FB"
    s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => false, "NSprotection" => true,
    "Permanentloss" => false, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea,  "year_base" => year_base, "total_yr" => total_yr)
    else
    s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false, "FSprotection" => false, "NSprotection" => false,
    "Permanentloss" => true, "Cont_list" => Cont_list,"base_list" => base_list,"Total_sample" =>Total_sample, "curtailed_gen" => curtailed_gen, "max_curt" => max_curt, "syncarea" => syncarea, "year_base" => year_base, "total_yr" => total_yr)
    end

    if period == "single"   push!(s, "multiperiod"=> false)
    elseif period == "multi"  push!(s, "multiperiod"=> true )
    end

    @assert ( s["FSprotection"] == true && (Prot_system == "FS_HDCCB" || Prot_system == "FS_MDCCB") ) || (s["NSprotection"] == true && (Prot_system == "NS_CB" || Prot_system == "NS_FB"))|| (Prot_system == "Permanentloss" && s["NSprotection"] == false && s["FSprotection"] == false)
    conv_rate = Int(data_cont["nw"]["1"]["convdc"]["1"]["Pacmax"]*100)
    if occursin("4bus", file)
        filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\4bus\\injection\\wthFCRlim\\",conv_rate,"MW\\",Prot_system,"_sPLdim.xlsx")
    elseif occursin("6bus", file)
        filepath = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\plots\\OPF\\4bus\\injection\\wthFCRlim\\",conv_rate,"MW\\",Prot_system,"_sPLdim.xlsx")
    end

    for (n,nw) in data_cont["nw"]
                for (r,reserves) in nw["reserves"]
                     if Prot_system == "FS_HDCCB"; reserves["Tcl"] = 0.148
                     elseif Prot_system == "FS_MDCCB"; reserves["Tcl"] = 0.293
                     elseif Prot_system == "NS_FB"; reserves["Tcl"] = 0.100
                      elseif Prot_system == "NS_CB"; reserves["Tcl"] = 0.150
                     elseif Prot_system == "Permanentloss"; reserves["Tcl"] = 100
                     end
                 end
     end

    resultDC1 = _PMACDC.run_acdcscopf_nocl_PLdim(data_cont, _PM.DCPPowerModel, gurobi, multinetwork = true;  setting = s)
    display_keyindictionary_OPF_PLdim(resultDC1, "isbuilt", "Pgg")
    # display(curtailed_gen)
    curtail, maxFFR, maxFCR, meanFFR, meanFCR = curtailment(data_cont, base_list, resultDC1, curtailed_gen)
     XLSX.openxlsx(filepath, mode="w") do xf
        sheet = xf[1]
        sheet["$(string("A",1))"] = data_cont["nw"]["1"]["reserves"]["2"]["H"]
        sheet["$(string("A",2))"] = resultDC1["solution"]["nw"]["1"]["FFR_Reserves"]
        sheet["$(string("A",3))"] = resultDC1["solution"]["nw"]["1"]["FCR_Reserves"]
        sheet["$(string("A",4))"] = resultDC1["solution"]["nw"]["1"]["Gen_cost"]
        sheet["$(string("A",5))"] = sum(resultDC1["solution"]["nw"]["1"]["Curt"])
        # sheet["$(string("A",5))"] = sum(curtail)
        sheet["$(string("A",6))"] = resultDC1["objective"]
        sheet["$(string("A",7))"] = maxFFR
        sheet["$(string("A",8))"] = maxFCR
        sheet["$(string("A",9))"] = resultDC1["objective_lb"]
        sheet["$(string("A",10))"] = meanFFR
        sheet["$(string("A",11))"] = meanFCR
    end

FCR_reserve_hr_2025 = Any[]; FCR_reserve_hr_2030 = Any[]; FCR_reserve_hr_2035 = Any[]; FCR_reserve_hr_2040 = Any[]; FCR_reserve_hr_2045 = Any[]; FCR_reserve_hr_2050 = Any[]
FFR_reserve_hr_2025 = Any[]; FFR_reserve_hr_2030 = Any[]; FFR_reserve_hr_2035 = Any[]; FFR_reserve_hr_2040 = Any[]; FFR_reserve_hr_2045 = Any[]; FFR_reserve_hr_2050 = Any[]
Inertia_2025 = Any[]; Inertia_2030 = Any[]; Inertia_2035 = Any[]; Inertia_2040 = Any[]; Inertia_2045 = Any[]; Inertia_2050 = Any[]

for b in year_base[1]
     push!(Inertia_2025, data_cont["nw"]["$b"]["reserves"]["2"]["H"])
     push!(FFR_reserve_hr_2025, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pff"])
     push!(FCR_reserve_hr_2025, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pgg"])
end

for b in year_base[2]
     push!(Inertia_2030, data_cont["nw"]["$b"]["reserves"]["2"]["H"])
     push!(FFR_reserve_hr_2030, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pff"])
     push!(FCR_reserve_hr_2030, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pgg"])
end

for b in year_base[3]
     push!(Inertia_2035, data_cont["nw"]["$b"]["reserves"]["2"]["H"])
     push!(FFR_reserve_hr_2035, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pff"])
     push!(FCR_reserve_hr_2035, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pgg"])
end
# for b in year_base[4]
#      push!(Inertia_2040, data_cont["nw"]["$b"]["reserves"]["2"]["H"])
#      push!(FFR_reserve_hr_2040, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pff"])
#      push!(FCR_reserve_hr_2040, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pgg"])
# end
# for b in year_base[5]
#      push!(Inertia_2045, data_cont["nw"]["$b"]["reserves"]["2"]["H"])
#      push!(FFR_reserve_hr_2045, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pff"])
#      push!(FCR_reserve_hr_2045, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pgg"])
# end
# for b in year_base[6]
#      push!(Inertia_2050, data_cont["nw"]["$b"]["reserves"]["2"]["H"])
#      push!(FFR_reserve_hr_2050, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pff"])
#      push!(FCR_reserve_hr_2050, resultDC1["solution"]["nw"]["$(b+1)"]["reserves"]["2"]["Pgg"])
# end


text:
using AllPlugins
using Test

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


text:
""" Code to compute Mie scattering parameters based on MieScatter julia
package

https://github.com/dronir/MieScatter.jl

"""

module MieParams
using FastGaussQuadrature

export compute

# Unfortunately the julia package manager does not really work for dependencies
# on unregistered packages.  And it seems that nothing really works for
# sub-packages so I could also not use them.  This was the only half-satisfactory
# way.
include("MieScatter/src/MieScatter.jl")
using .MieScatter

include("constants.jl")
const co = constants

"""
    g, ω₀, Qext = compute(r, λ, m; order=1000)

Compute the parameters `g`, `ω₀`, `Qext` for given radius of spherical particle,
`r`, wavelength `λ` and refractive index m (generally a complex number).
"""
function compute(r::Number, λ, m; order=1000)
    μ, w = gausslegendre(order)
    theta = acos.(μ)

    S, _, _, _ = compute_mie(size_parameter(r, λ), m, theta)

    # For some reason compute_mie returns Qext and Qback only when we ask
    # for forward and / or backward scattering
    _, Qsca, Qext, Qback = compute_mie(size_parameter(r, λ), m, [0, π])

    ω₀ = Qsca / Qext
    g = sum(μ .* w .* S[:, 1]) / sum(w .* S[:, 1])

    g, ω₀, Qext
end

function compute(r::AbstractVector, λ, m; kw...)
    w = compute.(r, λ, m; kw...)

    g, ω0, Qext = [[itm[i] for itm in w] for i in 1:3]
    g, ω0, Qext
end

end


text:
module Agents

export step!, AbstractModel, AbstractAgent, AbstractSpace, batchrunner, data_collector, grid, gridsize, move_agent!, add_agent!, coord_to_vertex, vertex_to_coord, get_node_contents, node_neighbors, nagents, return_activation_order, random_activation, as_added, visualize_data, add_agent_single!, move_agent_single!, kill_agent!, find_empty_nodes, find_empty_nodes_coords, id_to_agent, write_to_file, visualize_2D_agent_distribution, Random, combine_columns!, Node_iter, empty_nodes, is_empty, pick_empty, partial_activation, dummystep, visualize_1DCA, visualize_2DCA, SimpleGraph, batchrunner_parallel

using Distributed
using LightGraphs
using DataFrames
using VegaLite
using Random
import Base.Iterators.product
import StatsBase
using DataVoyager
using CSV
using ColorTypes
using GraphPlot
using Compose
using Cairo, Fontconfig
import Base.iterate
import Base.length
include("agents_component.jl")
include("model_component.jl")
include("space.jl")
include("scheduler.jl")
include("data_collector.jl")
include("batch_runner.jl")
include("visualization.jl")
include("CA1D.jl")
include("CA2D.jl")

end # module


text:
function _make_kkt_solver!(ws::COSMO.Workspace)

	ws.kkt_solver = ws.settings.kkt_solver(SparseMatrixCSC(ws.p.P),
						SparseMatrixCSC(ws.p.A),
						ws.settings.sigma,
						ws.ρvec)

end

function setup!(ws::COSMO.Workspace)

  allocate_set_memory!(ws)
	# scale problem data
	if ws.settings.scaling != 0
		scale_ruiz!(ws)
	end

	set_rho_vec!(ws)

	# create a KKT solver object populated with our data
	if(ws.flags.FACTOR_LHS)
		if ws.settings.verbose_timing
			ws.times.factor_time += @elapsed _make_kkt_solver!(ws)
		else
			_make_kkt_solver!(ws)
		end

	end
end

function allocate_set_memory!(ws)
  for (k, C) in enumerate(ws.p.C.sets)
    allocate_memory!(C)
  end
end


text:
# ---------------------------------------------------------------------------- #
#
#   intrico.jl
#
#   Geometry generation for lattice structures for AM
#
#   intrico
#   Spring 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

__precompile__()

"""
    intrico

Julia package to generate CAD of lattice structure for additive manufacturing.

Max Opgenoord

Spring/Fall 2018
"""

module intrico

    using Requires

    using StaticArrays

    using divido

    # Mesh datastructures
    export MeshF, MeshF2D, MeshF3D
    include("mesh/meshF.jl")

    # Lattice datastructure
    export Lattice
    include("lattice/lattice.jl")

    # Generate CAD
    export genCAD, genSTL, genSTLmodel, genCAD
    include("cad/ioSTL.jl")
    include("cad/genSTL.jl")
    include("cad/genSTLmodel.jl")
    @require egads include("cad/genCAD.jl")

end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# NDArray functions dynamically imported from libmxnet

function _invoke_mxfunction(func_handle::MX_handle, use_vars, scalars, mut_vars; kwargs...)
  names = String[string(entry[1]) for entry in kwargs]
  args = String[string(entry[2]) for entry in kwargs]
  @mxcall(:MXFuncInvokeEx,
          (MX_handle, Ptr{MX_handle}, Ptr{MX_float}, Ptr{MX_handle}, Cint, char_pp, char_pp),
          func_handle, use_vars, scalars, mut_vars, length(names), names, args)
end

@enum(LIBMX_FUNC_TYPE_MASK,
  NDARRAY_ARG_BEFORE_SCALAR = 1,
  ACCEPT_EMPTY_MUTATE_TARGET = (1 << 2)
)

# Import corresponding math functions from base so the automatically defined libmxnet
# functions can overload them
import Base: sqrt

"""
The libxmnet APIs are automatically imported from `libmxnet.so`. The functions listed
here operate on `NDArray` objects. The arguments to the functions are typically ordered
as

```julia
  func_name(arg_in1, arg_in2, ..., scalar1, scalar2, ..., arg_out1, arg_out2, ...)
```

unless `NDARRAY_ARG_BEFORE_SCALAR` is not set. In this case, the scalars are put before the input arguments:

```julia
  func_name(scalar1, scalar2, ..., arg_in1, arg_in2, ..., arg_out1, arg_out2, ...)
```

If `ACCEPT_EMPTY_MUTATE_TARGET` is set. An overloaded function without the output arguments will also be defined:

```julia
  func_name(arg_in1, arg_in2, ..., scalar1, scalar2, ...)
```

Upon calling, the output arguments will be automatically initialized with empty NDArrays.

Those functions always return the output arguments. If there is only one output (the typical situation), that
object (`NDArray`) is returned. Otherwise, a tuple containing all the outputs will be returned.
"""
function _get_ndarray_function_def(name::String)
  func_name = Symbol(name)

  func_def = quote
    function $func_name(::Type{<:NDArray}, args::NDArray...; out=nothing, kwargs...)
      if out != nothing
        output_vars = out
        if isa(output_vars, NDArray)
          output_vars = NDArray[output_vars]
        end
        num_outputs = length(output_vars)
      else
        output_vars = NDArray[]
        num_outputs = 0
      end

      args = collect(args)  # tuple to list
      if length(args) == 0
        args = MX_handle[]
      end

      output_handles_pp = if length(output_vars) > 0
        [map(x -> x.handle, output_vars)]
      else
        [Ptr{MX_handle}(C_NULL)]
      end
      num_outputs_p = [convert(Cint, num_outputs)]

      kw_keys_str = String[string(x[1]) for x in kwargs]
      kw_vals_str = String[dump_mx_param(x[2]) for x in kwargs]

      op_handle = _get_cached_libmx_op_handle($(name))
      @mxcall(:MXImperativeInvoke,
              (MX_handle, Cint, Ptr{MX_handle},
               Ptr{Cint}, Ptr{Ptr{MX_handle}},
               Cint, char_pp, char_pp),
              op_handle, length(args), args,
              num_outputs_p, output_handles_pp,
              length(kwargs), kw_keys_str, kw_vals_str)

      if out == nothing
        n = num_outputs_p[]
        hdls = unsafe_wrap(Array{MX_handle}, output_handles_pp[], n)
        xs = NDArray[NDArray(MX_NDArrayHandle(x)) for x in hdls]
        if n == 1
          return xs[]
        else
          return xs
        end
      else
        return out
      end
    end
  end

  func_def2 = quote
    function $func_name(args::NDArray...; out=nothing, kwargs...)
      $func_name(NDArray, args...; out=out, kwargs...)
    end
  end

  return func_def, func_def2
end

const _op_import_bl = [  # import black list; do not import these funcs
    "_full",   # we already have `mx.fill`
    "_ones",   # we already have `mx.ones`
    "_zeros",  # we already have `mx.zeros`
    "clip",
    "expand_dims",

    # arithmetic
    "_plus",
    "_minus",
    "_mod",
    "_mod_scalar",
    "_rmod_scalar",

    "dot",
    "max",
    "max_axis",
    "mean",
    "min",
    "min_axis",
    "prod",
    "reshape",
    "sum",
    "transpose",

    # trigonometric
    "sin",
    "cos",
    "tan",
    "arcsin",
    "arccos",
    "arctan",

    # hyperbolic
    "sinh",
    "cosh",
    "tanh",
    "arcsinh",
    "arccosh",
    "arctanh",

    # activation
    "sigmoid",
    "relu",
    "softmax",
    "log_softmax",

    # broadcast
    "broadcast_add",
    "broadcast_plus",
    "broadcast_minus",
    "broadcast_sub",
    "broadcast_mul",
    "broadcast_div",
    "broadcast_mod",
    "broadcast_power",
    "broadcast_equal",
    "broadcast_not_equal",
    "broadcast_greater",
    "broadcast_greater_equal",
    "broadcast_lesser",
    "broadcast_lesser_equal",
    "broadcast_maximum",
    "broadcast_minimum",
    "broadcast_to",
    "broadcast_axis",
    "broadcast_axes",
    "broadcast_hypot",

    # reduction
    "argmax",
    "argmin",
]

macro _import_ndarray_functions()
  names = filter(n -> ∉(lowercase(n), _op_import_bl), _get_libmx_op_names())

  func_exprs = map(names) do name
    op_handle = _get_libmx_op_handle(name)

    desc, key_narg = _get_libmx_op_description(name, op_handle)
    func_def, func_def2 = _get_ndarray_function_def(name)

    func_name = Symbol(name)

    import_expr = _import_expr(func_name)

    quote
      $import_expr
      $func_def
      @doc $desc
      $func_def2
    end
  end

  esc(quote
    $(func_exprs...)
  end)
end

@_import_ndarray_functions


text:
export isboolean

function isboolean(str::AbstractString, loose=false::Bool)::Bool
    strictBooleans = ["true", "false", "1", "0"]
    looseBooleans = [strictBooleans..., "yes", "no"]

    return loose ? lowercase(str) in looseBooleans : str in strictBooleans
end


text:
using CSV, DataFrames

println("Hello SWUFE!")

max_stars_repo_name: sukanka/LaTeX-SWUFE-Bachelor-Thesis

max_issues_repo_name: sukanka/LaTeX-SWUFE-Bachelor-Thesis

max_forks_repo_name: sukanka/LaTeX-SWUFE-Bachelor-Thesis


text:
using GeoInterface, ADCPDataProcessing, MapTiles

struct CreekPoint <: AbstractPoint
    creek::Creek
end

GeoInterface.coordinates(cf::CreekPoint) = _get_geom(cf.creek)

function _get_geom(creek::Creek)
    M = ADCPDataProcessing.metadataload(creek)
    lat = M["gps"]["lat"]
    lon = M["gps"]["lon"]
    [lon,lat]
end

function _get_coords(creek::Creek)
    M = ADCPDataProcessing.metadataload(creek)
    E = M["gps"]["east"]
    N = M["gps"]["north"]
    [E,N]
end

function _get_bbox(m::MultiPoint,α=0.05)
    X = coordinates(m)
    Z = [map(x->x[1],X) map(x->x[2],X)]
    Zm = minimum(Z,1)
    Zp = maximum(Z,1)

    Zc = (Zm+Zp)/2
    Zd = (Zp-Zm)/2
    [(Zc+(1+α)*Zd)'; (Zc-(1+α)*Zd)']
end

function fetchtilesmp(m::MultiPoint,z::Int,maxtiles::Int,provider::MapTiles.AbstractProvider,α=0.05)
    MapTiles.fetchmap(_get_bbox(m)...,z,maxtiles=maxtiles,provider=provider)
end

@userplot CreekMap

@recipe function f(cm::CreekMap)
    creeks,z,maxtiles,provider = cm.args

    m = MultiPoint(CreekPoint.(creeks))

    b = fetchtilesmp(m,z,maxtiles,provider)

    framestyle --> :none
    ticks --> nothing
    aspect_ratio --> 1
    yflip --> true

    @series begin
        b
    end

    @series begin
        MapTiles.project(b,m)
    end
end


text:
using Base.Test
import Ogre
import WaterData

module test_structure_resources
import Ogre
pressure = Ogre.pressure
temperature = Ogre.temperature

type PressureEOS <: Ogre.EOS end
Ogre.@addEOSCall PressureEOS
(::PressureEOS)(P) = 4100. + 0.00161*(P^0.541)
Ogre.istempdependent(::PressureEOS) = false

type PTEOS <: Ogre.EOS end
Ogre.@addEOSCall PTEOS
(::PTEOS)(P, T) = 4100. + 0.00161*(P^0.541)
# the effect of T on ρ here is arbitrary - we just want increasing
# T to mean decreasing ρ
Ogre.istempdependent(::PTEOS) = true
end

@testset "Structure equations" begin
    res = test_structure_resources

    # realistic values taken from PREM at r=1000 km from the centre
    mass = 5.4e22      # kg enclosed
    radius = 1e6       # m
    pressure = 340e9   # Pa

    @testset "with no temperature dependence" begin
        realistic_values = Ogre.ValueSet(mass, radius, pressure)

        @testset "using an analytic EOS" begin
            eos = res.PressureEOS()

            masscontinuity = Ogre.MassContinuity(eos)
            pressurebalance = Ogre.PressureBalance()

            @testset "Zero gradients at the r=0 limit" begin
                zero_values = zero(Ogre.MassRadiusPressure)
                @test pressurebalance(zero_values) == 0
                @test masscontinuity(zero_values) == 0
            end

            @testset "Correct signs for the structure equations" begin
                # change in pressure is negative outwards
                @test pressurebalance(realistic_values) < 0
                # change in mass is positive outwards
                @test masscontinuity(realistic_values) > 0
            end

            @testset "Attempting to go over the singularity at r=0 returns zero" begin
                negative_radius = Ogre.ValueSet(5.4e22, -1e6, 340e9)
                negative_pressure = Ogre.ValueSet(5.4e22, 1e6, -340e9)
                negative_mass = Ogre.ValueSet(-5.4e22, 1e6, 340e9)
                map([negative_radius, negative_pressure, negative_mass]) do vs
                    @test pressurebalance(vs) == 0
                    @test masscontinuity(vs) == 0
                end
            end
        end
    end

    @testset "with temperature dependence" begin
        # realistic values taken from PREM at r=1000 km from the centre
        mass = 5.4e22      # kg enclosed
        radius = 1e6       # m
        pressure = 340e9   # Pa
        temperature = 5000 # K, this value is only rough
        realistic_values = Ogre.ValueSet(mass, radius, pressure, temperature)

        @testset "using an analytic EOS" begin
            eos = res.PTEOS()

            # arbitrary heat capacity
            Cₚ = WaterData.ConstantHeatCapacity(100)

            masscontinuity = Ogre.MassContinuity(eos)
            pressurebalance = Ogre.PressureBalance()
            temperaturegradient = Ogre.TemperatureGradient(eos, Cₚ)

            @testset "Zero gradients at the r=0 limit" begin
                zero_values = zero(Ogre.PhysicalValues)
                @test pressurebalance(zero_values) == 0
                @test masscontinuity(zero_values) == 0
                @test temperaturegradient(zero_values) == 0
            end

            @testset "Correct signs for the structure equations" begin
                # change in pressure is negative outwards
                @test pressurebalance(realistic_values) < 0
                # change in mass is positive outwards
                @test masscontinuity(realistic_values) > 0
                # change in temperature is negative outwards
                @test temperaturegradient(realistic_values) < 0
            end

            @testset "Attempting to go over the singularity at r=0 returns zero" begin
                negative_radius = Ogre.ValueSet(5.4e22, -1e6, 340e9, 5000.)
                negative_pressure = Ogre.ValueSet(5.4e22, 1e6, -340e9, 5000.)
                negative_mass = Ogre.ValueSet(-5.4e22, 1e6, 340e9, 5000.)
                negative_temperature = Ogre.ValueSet(5.4e22, -1e6, 340e9, -5000.)
                map([negative_radius, negative_pressure,
                     negative_mass, negative_temperature]) do vs
                         @test pressurebalance(vs) == 0
                         @test masscontinuity(vs) == 0
                         @test temperaturegradient(vs) == 0
                     end
            end
        end
    end
end

@testset "Planetary structure types" begin
    res = test_structure_resources
    @testset "Boundary values" begin
        bvs = Ogre.BoundaryValues(1,2,3)
        bvsT = Ogre.BoundaryValues(1,2,3,4)
        @test Ogre.radius(bvs) == 2
        @test Ogre.temperature(bvsT) == 4
    end

    @testset "Centre and surface values" begin
        mc, ms = 1, 2
        rc, rs = 3, 4
        Pc, Ps = 5, 6
        Tc, Ts = 7, 8
        npoints = 10

        m = linspace(ms, mc, npoints)
        r = linspace(rs, rc, npoints)
        P = linspace(Ps, Pc, npoints)
        T = linspace(Ts, Tc, npoints)

        s1 = Ogre.PlanetStructure(m, r, P)
        s2 = Ogre.PlanetStructure(m, r, P, T)

        @test Ogre.mass(Ogre.centre(s1)) == mc
        @test Ogre.mass(Ogre.centre(s2)) == mc
        @test Ogre.mass(Ogre.surface(s1)) == ms
        @test Ogre.mass(Ogre.surface(s2)) == ms
        @test Ogre.radius(Ogre.centre(s1)) == rc
        @test Ogre.radius(Ogre.centre(s2)) == rc
        @test Ogre.radius(Ogre.surface(s1)) == rs
        @test Ogre.radius(Ogre.surface(s2)) == rs
        @test Ogre.pressure(Ogre.centre(s1)) == Pc
        @test Ogre.pressure(Ogre.centre(s2)) == Pc
        @test Ogre.pressure(Ogre.surface(s1)) == Ps
        @test Ogre.pressure(Ogre.surface(s2)) == Ps
        @test Ogre.temperature(Ogre.centre(s2)) == Tc
        @test Ogre.temperature(Ogre.surface(s2)) == Ts
    end

    @testset "Planet system and solution setup" begin
        M = 5.972e24
        R = 6.3781e6
        Psurf = 1e5
        solution_grid = linspace(M, 0, 5)
        radius_bracket = [0, 10] * M

        pressurebalance = Ogre.PressureBalance()

        @testset "No temperature dependence" begin
            bvs = Ogre.BoundaryValues(M, R, Psurf)
            eos = res.PressureEOS()
            masscontinuity = Ogre.MassContinuity(eos)
            structure = Ogre.EquationSet([masscontinuity,
                                          pressurebalance])
            system = Ogre.PlanetSystem(M, eos, bvs, solution_grid,
                                       radius_bracket)
            structure = Ogre.blank_structure(system)

            @test length(Ogre.mass(structure)) == 5
            @test length(Ogre.radius(structure)) == 5
            @test length(Ogre.pressure(structure)) == 5
            @test_throws MethodError Ogre.temperature(structure)
        end

        @testset "Temperature dependence" begin
            Tsurf = 300
            bvs = Ogre.BoundaryValues(M, R, Psurf, Tsurf)
            eos = res.PTEOS()
            Cₚ = WaterData.ConstantHeatCapacity(1000)
            masscontinuity = Ogre.MassContinuity(eos)
            temperaturegradient = Ogre.TemperatureGradient(eos, Cₚ)
            structure = Ogre.EquationSet([masscontinuity,
                                          pressurebalance,
                                          temperaturegradient])
            system = Ogre.PlanetSystem(M, eos, Cₚ, bvs, solution_grid,
                                       radius_bracket)
            structure = Ogre.blank_structure(system)

            @test length(Ogre.mass(structure)) == 5
            @test length(Ogre.radius(structure)) == 5
            @test length(Ogre.pressure(structure)) == 5
            @test length(Ogre.temperature(structure)) == 5
        end
    end
end


text:
include("markov.jl")

include("savecosts.jl")

using PyPlot

function qvalueiteration(cost_akt)
    nb_time_steps = 169
    nb_stages = 18
    nb_actions = 2*3^3
    qval = zeros(nb_time_steps, nb_stages, nb_actions)
    # need to fill qval with cost values
    opt_action = zeros(nb_time_steps, nb_stages)

    alpha = 0.9
    beta = 1
    while alpha > 0.01
        print(alpha, "\n")
        for t = 2:nb_time_steps
		    for k = 1:nb_stages
		    	state = ind_to_vect(k,3,3)[1]
		    	nb_actions = state[1] == 1 ? 2*3^3 : 3^3
                min_q = 0
		    	for a = 1:nb_actions
		    		actions = ind_to_vect(a,4,3)
		    		if actions[1] == 1
                        min_q = minimum(qval[t-1,:,:])
                    end
		    		if actions[1] == 2
                        min_q = minimum(qval[t-1,10:18,28:54])
                    end
                    qval[t,k,a] = qval[t,k,a] + alpha*(cost_akt[a,k,t] + min_q - qval[t,k,a])
                    opt_action[t,k] = indmin(qval[t,k,:])
               end
           end
        end
        # hyperbolic decrease
        alpha = alpha*beta/(alpha + beta)
    end
    return opt_action
end

# tic()

# cost_akt = loadcost()
# actions = qvalueiteration(cost_akt)

# function run_and_save_bestactions()
#     cost_akt = loadcost()
#     opt_actions = qvalueiteration(cost_akt)
#     save("./data/opt_actions.jld", "opt_actions", opt_actions)
# end

# toc()

# run_and_save_bestactions()

raw = load("./data/opt_actions.jld")
opt_actions = raw["opt_actions"]

nb_time_steps = size(opt_actions)[1]
nb_stages = size(opt_actions)[2]

state_min = 1
power_cut = zeros(nb_time_steps)

for t = 2:nb_time_steps
    iter_state = rand(state_min:18)
    # once we set the load OFF, we have to keep it OFF till the end
    # if (iter_state[1] > 9)
    #     state_min = 10
    # end
    iter_action = ind_to_vect(opt_actions[t,iter_state], 4, 3)
    action_h3 = iter_action[4]
    # print("House 3 : ", action_h3, "\n")
    if action_h3 == 1
        power_cut[t] = 0
    elseif action_h3 == 2
        power_cut[t] = 0.5
    else
        power_cut[t] = 1
    end
end

data = matread("./data/data_4nodes.mat")
p_pred = data["pSolarMax"][:, 3]

p_pred = p_pred*1000

p_grid = power_cut.*p_pred

fig = plot()

times = [i for i=1:nb_time_steps]
# plot!(times, p_grid,
#    label=["Power cut"],
#    xlabel=("Time step"),
#    ylabel=("Power in kW"))
# plot!(times, p_grid, marker="o", label=["Power forecast"])

# histogram(p_grid, nbins = 20)

bar(times, p_grid)
xticks([],[])
xlabel("Time step")
ylabel("Curtailment [kW]")
show()

# show()
savefig(fig, "../fig.pdf")

max_stars_repo_name: qlete/ANManagement

max_issues_repo_name: qlete/ANManagement

max_forks_repo_name: qlete/ANManagement


text:
using CompScienceMeshes, BEAST
using LinearAlgebra

Γ = meshcuboid(0.5, 1.0, 1.0, 0.075)
CompScienceMeshes.translate!(Γ, point(-0.25, -0.5, -0.5))
# Γ = meshsphere(1.0, 0.075)
X = raviartthomas(Γ)
@show numfunctions(X)

κ,  η  = 6.0, 1.0
κ′, η′ = 2.4κ, η/2.4

# κ = 1.0
# η = 1.0
# κ′ = 1.5*κ
# η′ = η/1.5

T  = Maxwell3D.singlelayer(wavenumber=κ)
T′ = Maxwell3D.singlelayer(wavenumber=κ′)
K  = Maxwell3D.doublelayer(wavenumber=κ)
K′ = Maxwell3D.doublelayer(wavenumber=κ′)

# d = normalize(x̂+ŷ+ẑ)
d = ẑ
# p = normalize(d × ẑ)
p = x̂
E = Maxwell3D.planewave(direction=d, polarization=p, wavenumber=κ)
H = -1/(im*κ*η)*curl(E)

e, h = (n × E) × n, (n × H) × n

@hilbertspace j m
@hilbertspace k l

α, α′ = 1/η, 1/η′
pmchwt = @discretise(
    (η*T+η′*T′)[k,j] -      (K+K′)[k,m] +
         (K+K′)[l,j] + (α*T+α′*T′)[l,m] == -e[k] - h[l],
    j∈X, m∈X, k∈X, l∈X)

u = solve(pmchwt)

nX = numfunctions(X)
uj = u[1:nX]
um = u[nX+1:end]

Θ, Φ = range(0.0,stop=2π,length=100), 0.0
ffpoints = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]

# Don't forget the far field comprises two contributions
ffm = potential(MWFarField3D(κ*im), ffpoints, um, X)
ffj = potential(MWFarField3D(κ*im), ffpoints, uj, X)
ff = η*im*κ*ffj + im*κ*cross.(ffpoints, ffj)

using Plots
plot(xlabel="theta")
plot!(Θ,norm.(ff),label="far field")

import Plotly
using LinearAlgebra
fcrj, _ = facecurrents(uj,X)
fcrm, _ = facecurrents(um,X)
Plotly.plot(patch(Γ, norm.(fcrj)))
Plotly.plot(patch(Γ, norm.(fcrm)))

function nearfield(um,uj,Xm,Xj,κ,η,points,
    Einc=(x->point(0,0,0)),
    Hinc=(x->point(0,0,0)))

    K = BEAST.MWDoubleLayerField3D(wavenumber=κ)
    T = BEAST.MWSingleLayerField3D(wavenumber=κ)

    Em = potential(K, points, um, Xm)
    Ej = potential(T, points, uj, Xj)
    E = -Em + η * Ej + Einc.(points)

    Hm = potential(T, points, um, Xm)
    Hj = potential(K, points, uj, Xj)
    H = 1/η*Hm + Hj + Hinc.(points)

    return E, H
end

Z = range(-1,1,length=100)
Y = range(-1,1,length=100)
nfpoints = [point(0,y,z) for z in Z, y in Y]

# nfm_in = potential(BEAST.MWDoubleLayerField3D(wavenumber=κ′), nfpoints, um, X)
# nfj_in = potential(BEAST.MWSingleLayerField3D(wavenumber=κ′), nfpoints, uj, X)
# E_in = nfm_in - η′ * nfj_in

# nfm_ex = potential(BEAST.MWDoubleLayerField3D(wavenumber=κ), nfpoints, um, X)
# nfj_ex = potential(BEAST.MWSingleLayerField3D(wavenumber=κ), nfpoints, uj, X)
# E_ex = -nfm_ex + η * nfj_ex + E.(nfpoints)

E_ex, H_ex = nearfield(um,uj,X,X,κ,η,nfpoints,E,H)
E_in, H_in = nearfield(-um,-uj,X,X,κ′,η′,nfpoints)

E_tot = E_in + E_ex
H_tot = H_in + H_ex

contour(real.(getindex.(E_tot,1)))
heatmap(Z, Y, real.(getindex.(E_tot,1)))
plot(real.(getindex.(E_tot[:,51],1)))

# Compute also the near magnetic field
# Hm_ex = 1/η*potential(BEAST.MWSingleLayerField3D(wavenumber=κ), nfpoints, um, X)
# Hj_ex = potential(BEAST.MWDoubleLayerField3D(wavenumber=κ), nfpoints, uj, X)
# H_ex = Hm_ex + Hj_ex + H.(nfpoints)

# Hm_in = -1/η′*potential(BEAST.MWSingleLayerField3D(wavenumber=κ′), nfpoints, um, X)
# Hj_in = -potential(BEAST.MWDoubleLayerField3D(wavenumber=κ′), nfpoints, uj, X)
# H_in = Hm_in + Hj_in
# Hnf = H_ex + H_in

contour(real.(getindex.(H_tot,2)))
heatmap(Z, Y, real.(getindex.(H_tot,2)))
plot!(real.(getindex.(H_tot[:,51],2)))


text:
module MTL

import Libdl
using CEnum

const cmt_lib = joinpath(@__DIR__, "..", "..", "deps", "libcmt_lib.dylib")

## source code includes

# Basic types
include("libcmt.jl")
include("error.jl")
# convert from 1 based indexing to 0 based indexing
Base.convert(::Type{NsRange}, range::UnitRange{T}) where T <: Integer =
	NsRange(first(range), length(range))
# used for byte ranges.
Base.convert(::Type{NsRange}, range::StepRange{T}) where T <: Integer =
	NsRange(first(range)-step(range), length(range)*step(range))

# low-level autogeneraed wrappers
export MtSize, MtlDim, MtlDim3
include("libcmt.jl")

include("helpers.jl")

# julia wrappers
include("storage_type.jl")
include("resource.jl")
include("device.jl")
include("compile-options.jl")
include("library.jl")
include("function.jl")
include("events.jl")
include("heap.jl")
include("buffer.jl")
include("cmd-queue.jl")
include("cmd-buffer.jl")
include("compute-pipeline-state.jl")
include("buffer_argument.jl")

include("command_enc.jl")
include("command_enc_blit.jl")
include("command_enc_compute.jl")

include("reflection/reflection_compute.jl")

include("size.jl")

end # module


text:
using JWTs
using Test
using JSON
using MbedTLS

const test_payload_data = [
    JSON.parse("""{
        "jti": "0b821616-0a5f-47f3-af00-8caf03619303",
        "exp": 1543351759,
        "nbf": 0,
        "iat": 1543315759,
        "iss": "https://example.com/auth/",
        "aud": "portal",
        "sub": "b1df5448-a16b-4a13-b03b-2213d56ea1b5",
        "typ": "Bearer",
        "azp": "portal",
        "auth_time": 1543315759,
        "session_state": "f196425d-226b-4e6d-bc81-feecb276f424",
        "acr": "1",
        "allowed-origins": [ "" ],
        "realm_access": { "roles": [ "uma_authorization" ] },
        "resource_access": {
            "broker": { "roles": [ "read-token" ] },
            "account": { "roles": [ "manage-account", "manage-account-links", "view-profile" ] }
        },
        "preferred_username": "chhhhhhhhhhhhhhhhhhhhhhhhhaaaaaaaaaaaaabbb"
    }"""),
    JSON.parse("""{
        "iss": "https://auth2.juliacomputing.io/dex",
        "sub": "ChUxjfgsajfurjsjdut0483672kdhgstgy283jssZQ",
        "aud": "example-audience",
        "exp": 1536080651,
        "iat": 1535994251,
        "nonce": "1777777777777aaaaaaaaabbbbbbbbbb",
        "at_hash": "222222-G-JJJJJJJJJJJJJ",
        "email": "user@example.com",
        "email_verified": true,
        "name": "Example User"
    }""")
]

function print_header(msg)
    println("")
    println("-"^60)
    println(msg)
    println("-"^60)
end

function test_and_get_keyset(url)
    print_header("keyset: $url")

    keyset = JWKSet(url)
    @test length(keyset.keys) == 0

    refresh!(keyset)
    @test length(keyset.keys) > 0
    for (k,v) in keyset.keys
        println("    ", k, " => ", v.key)
    end

    keyset
end

function test_in_mem_keyset(template)
    print_header("keyset: $template")
    keyset = JWKSet(JSON.parse(read(template, String))["keys"])
    @test length(keyset.keys) == 4
    for (k,v) in keyset.keys
        println("    ", k, " => ", v.key)
    end
end

function test_signing_keys(keyset, signingkeyset)
    for k in keys(keyset.keys)
        for d in test_payload_data
            jwt = JWT(; payload=d)
            @test !issigned(jwt)
            sign!(jwt, signingkeyset, k)
            @test issigned(jwt)
            @test isvalid(jwt)
            @test isverified(jwt)

            println("    JWT: ", jwt)
            jwt2 = JWT(; jwt=string(jwt))
            @test issigned(jwt2)
            @test !isverified(jwt2)
            @test isvalid(jwt2) === nothing
            @test validate!(jwt, keyset, k)
            @test issigned(jwt)
            @test isvalid(jwt)
            @test isverified(jwt)

            jwt2 = JWT(; jwt=string(jwt))
            @test issigned(jwt2)
            @test !isverified(jwt2)
            @test isvalid(jwt2) === nothing
            invalidkey = first(filter(x->x!=k, keys(keyset.keys)))
            @test !validate!(jwt2, keyset, invalidkey)
            @test issigned(jwt2)
            @test !isvalid(jwt2)
            @test isverified(jwt2)
        end
    end
end

function test_signing_asymmetric_keys(keyset_url)
    print_header("signing asymmetric keys")
    keyset = JWKSet(keyset_url)
    refresh!(keyset)
    signingkeyset = deepcopy(keyset)
    for k in keys(signingkeyset.keys)
        keyfile = joinpath(dirname(keyset_url), "$k.private.pem")
        if startswith(keyfile, "file://")
            keyfile = keyfile[8:end]
        end
        signingkeyset.keys[k] = JWKRSA(signingkeyset.keys[k].kind, MbedTLS.parse_keyfile(keyfile))
    end
    test_signing_keys(keyset, signingkeyset)
end

function test_signing_symmetric_keys(keyset_url)
    print_header("signing symmetric keys")
    keyset = test_and_get_keyset(keyset_url)
    test_signing_keys(keyset, keyset)
end

test_and_get_keyset("https://www.googleapis.com/oauth2/v3/certs")
test_signing_symmetric_keys("file://" * joinpath(@__DIR__, "keys", "oct", "jwkkey.json"))
test_in_mem_keyset(joinpath(@__DIR__, "keys", "oct", "jwkkey.json"))
test_signing_asymmetric_keys("file://" * joinpath(@__DIR__, "keys", "rsa", "jwkkey.json"))


text:
function statesyms(x::Array{String,1})
  y = zeros(Sym,length(x))
  for i in eachindex(x)
    y[i] = symbols(x[i],real=true)
  end
  return y
end

"""
Create variables that are functions of time.  Used to define the positional DOF of a body.
"""
macro StateVars(x...)
  show(x)
  q = Expr(:block)
  ss = []
  for s in x
    push!(ss,s)

    # Set the symbol equal to a symbolic function of time
    push!(q.args,
      Expr(:(=), s,
        Expr(:call,
          Expr(:call, :(SymFunction), string(s))
        ,:t)
      )
    )

  end
  eval(Main, q)
end

"""
Create variables that depend on time
"""
function AddTimeDependence(varname::String)
  return SymFunction(varname)(t)
end
function AddTimeDependence(x::Sym)
  s = string(x)
  m = match(r"\(",s)
  if m != nothing
    s = s[1:m.offset-1]
  end
  return SymFunction(s)(t)
end
function AddTimeDependence(varvector::Vector)
  v = zeros(Sym,size(varvector))
  for i in eachindex(varvector)
    v[i] = AddTimeDependence(varvector[i])
  end
  return v
end
function AddTimeDependence(vardict::Dict)
  d = deepcopy(vardict)
  for key in keys(d)
    d[key] = AddTimeDependence(vardict[key])
  end
  return d
end

"""
Remove dependence on time from SymPy variables and assume they are real.
"""
function RemoveTimeDependence(x::Sym)
  s = string(x)
  m = match(r"\(",s)
  if m != nothing
    s = s[1:m.offset-1]
  end
  return symbols(s, real = true)
end
function RemoveTimeDependence(x::Array{Sym,1})
  v = zeros(Sym,size(x))
  for i in x
    v[i] = RemoveTimeDependence(x[i])
  end
  return v
end
function RemoveTimeDependence(x::Dict)
  d = deepcopy(x)
  for key in keys(d)
    d[key] = RemoveTimeDependence(x[key])
  end
  return d
end


text:
#
#
# From base, but copied here to make sure we don't fail bacause base changed
function my_gcd(a::T, b::T) where T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128}
    a == 0 && return abs(b)
    b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && throw(OverflowError("gcd($a, $b) overflows"))
    r % T
end

# make sure source code is unavailable for this method:
include_string(@__MODULE__, """
function my_gcd_noinfo(a::T, b::T) where T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128}
    a == 0 && return abs(b)
    b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && throw(OverflowError("gcd(\$a, \$b) overflows"))
    r % T
end
""", "nope.jl")

@testset "UI" begin
    if Sys.isunix() && VERSION >= v"1.1.0"
        Debugger._print_full_path[] = false
        using TerminalRegressionTests

        function run_terminal_test(frame, commands, validation)
            TerminalRegressionTests.automated_test(joinpath(@__DIR__, validation), commands) do emuterm
                repl = REPL.LineEditREPL(emuterm, true)
                repl.interface = REPL.setup_interface(repl)
                repl.specialdisplay = REPL.REPLDisplay(repl)
                RunDebugger(frame, repl, emuterm)
            end
        end

        CTRL_C = "\x3"
        EOT = "\x4"
        UP_ARROW = "\e[A"

        run_terminal_test(@make_frame(my_gcd(10, 20)),
                          ["n\n","`", "a\n", UP_ARROW, UP_ARROW, CTRL_C, EOT],
                          "ui/history_gcd.multiout")

        if v"1.1">= VERSION < v"1.2"
            run_terminal_test(@make_frame(my_gcd_noinfo(10, 20)),
                            ["n\n","`", "a\n", UP_ARROW, UP_ARROW, CTRL_C, EOT],
                             "ui/history_noinfo.multiout")
        else
            @warn "Skipping tests for IR display due to mismatched Julia versions."
        end
    else
        @warn "Skipping UI tests on non unix systems"
    end
end

# Completions
function test_complete(c, s)
    c, r, s = Debugger.completions(c, s, lastindex(s))
    return unique!(map(REPL.REPLCompletions.completion_text, c)), r, s
end

module F
    local_var = 1
    f(x) = x
end

@testset "REPL completions" begin
    frame = JuliaInterpreter.enter_call_expr(:($(F.f)(1)))
    state = dummy_state(frame)
    prov = Debugger.DebugCompletionProvider(state)

    c, r = test_complete(prov, "local")
    @test "local_var" in c
end


text:
using BenchmarkTools
using DataFrames
using CSV
using Statistics
using StatsBase: sample, quantile
using Plots
using Plots: colormap
using Base.Threads: @threads

using Revise
using EvoTrees
using EvoTrees: sigmoid, logit
using EvoTrees: softmax
using EvoTrees: update_grads!, get_gain, TrainNode, SplitInfo, get_edges, binarize, find_bags, grow_tree, find_split_static!
using EvoTrees: pred_leaf, softmax
using Flux: onehot

# prepare a dataset
iris = CSV.read("./data/iris.csv")
names(iris)

features = iris[[:PetalLength, :PetalWidth, :SepalLength, :SepalWidth]]
X = convert(Matrix, features)
Y = iris[:Species]
values = sort(unique(Y))
dict = Dict{String, Int}(values[i] => i for i in 1:length(values))
Y = map((x) -> dict[x], Y)

# train-eval split
𝑖 = collect(1:size(X,1))
𝑗 = collect(1:size(X,2))
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]
𝑖 = collect(1:size(X_train,1))
# scatter(X_train[:,1], X_train[:,2], color=Y_train, legend=nothing)
# scatter(X_eval[:,1], X_eval[:,2], color=Y_eval, legend=nothing)

##################################
# Step by step development
##################################
# set parameters
params1 = EvoTreeRegressor(
    loss=:softmax, metric=:mlogloss,
    nrounds=1, nbins=16,
    λ = 0.0, γ=0.0, η=0.3,
    max_depth = 3, min_weight = 1.0,
    rowsample=1.0, colsample=1.0,
    K = 3, seed=44)

# initial info
@time δ, δ² = zeros(SVector{params1.K, Float64}, size(X_train, 1)), zeros(SVector{params1.K, Float64}, size(X_train, 1))
𝑤 = zeros(SVector{1, Float64}, size(X_train, 1)) .+ 1
pred = zeros(SVector{params1.K,Float64}, size(X_train,1))
@time update_grads!(params1.loss, params1.α, pred, Y_train, δ, δ², 𝑤)
∑δ, ∑δ², ∑𝑤 = sum(δ[𝑖]), sum(δ²[𝑖]), sum(𝑤[𝑖])
@time gain = get_gain(params1.loss, ∑δ, ∑δ², ∑𝑤, params1.λ)

# initialize train_nodes
train_nodes = Vector{TrainNode{params1.K, Float64, BitSet, Array{Int64, 1}, Int}}(undef, 2^params1.max_depth-1)
for node in 1:2^params1.max_depth-1
    train_nodes[node] = TrainNode(0, SVector{params1.K, Float64}(fill(-Inf, params1.K)), SVector{params1.K, Float64}(fill(-Inf, params1.K)), SVector{1, Float64}(fill(-Inf, 1)), -Inf, BitSet([0]), [0])
    # train_nodes[feat] = TrainNode(0, fill(-Inf, params1.K), fill(-Inf, params1.K), -Inf, -Inf, BitSet([0]), [0])
end

# initializde node splits info and tracks - colsample size (𝑗)
splits = Vector{SplitInfo{params1.K, Float64, Int}}(undef, size(𝑗, 1))
for feat in 1:size(𝑗, 1)
    splits[feat] = SplitInfo{params1.K, Float64, Int}(gain, SVector{params1.K, Float64}(zeros(params1.K)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{1, Float64}(zeros(1)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{1, Float64}(zeros(1)), -Inf, -Inf, 0, feat, 0.0)
end

# binarize data and create bags
@time edges = get_edges(X_train, params1.nbins)
@time X_bin = binarize(X_train, edges)
@time bags = Vector{Vector{BitSet}}(undef, size(𝑗, 1))
function prep(X_bin, bags)
    @threads for feat in 1:size(𝑗, 1)
         bags[feat] = find_bags(X_bin[:,feat])
    end
    return bags
end
@time bags = prep(X_bin, bags)

# initialize histograms
feat=1
hist_δ = Vector{Vector{SVector{params1.K, Float64}}}(undef, size(𝑗, 1))
hist_δ² = Vector{Vector{SVector{params1.K, Float64}}}(undef, size(𝑗, 1))
hist_𝑤 = Vector{Vector{SVector{1, Float64}}}(undef, size(𝑗, 1))
for feat in 1:size(𝑗, 1)
    hist_δ[feat] = zeros(SVector{params1.K, Float64}, length(bags[feat]))
    hist_δ²[feat] = zeros(SVector{params1.K, Float64}, length(bags[feat]))
    hist_𝑤[feat] = zeros(SVector{1, Float64}, length(bags[feat]))
end

# grow single tree
@time train_nodes[1] = TrainNode(1, ∑δ, ∑δ², ∑𝑤, gain, BitSet(𝑖), 𝑗)
@time pred_leaf_ = pred_leaf(params1.loss, train_nodes[1], params1, δ²)
# @btime pred_leaf_ = pred_leaf($params1.loss, $train_nodes[1], $params1, $δ²)
@time tree = grow_tree(bags, δ, δ², 𝑤, hist_δ, hist_δ², hist_𝑤, params1, train_nodes, splits, edges, X_bin)

# feat = 1
# typeof(bags[feat][1])
# train_nodes[1] = TrainNode(1, ∑δ, ∑δ², ∑𝑤, gain, BitSet(𝑖), 𝑗)
# find_split_turbo!(bags[feat], view(X_bin,:,feat), δ, δ², 𝑤, ∑δ, ∑δ², ∑𝑤, params1, splits[feat], tracks[feat], edges[feat], train_nodes[1].𝑖)

pred = predict(tree, X_train, params1.K)
for i in eachindex(pred)
    pred[i] = exp.(pred[i]) / sum(exp.(pred[i]))
end
pred_int = zeros(Int, length(pred))
for i in eachindex(pred)
    pred_int[i] = findmax(pred[i])[2]
end
sum(pred_int .== Y_train)

params1 = EvoTreeRegressor(
    loss=:softmax, metric=:mlogloss,
    nrounds=20, nbins=16,
    λ = 0.0, γ=1e-5, η=0.3,
    max_depth = 3, min_weight = 1.0,
    rowsample=1.0, colsample=1.0,
    K = 3, seed=44)

@time model = grow_gbtree(X_train, Y_train, params1, print_every_n = Inf)
# @time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 1)

sum(Y_train.==3)/length(Y_train)
@time pred_train = predict(model, X_train)
@time pred_eval = predict(model, X_eval)

maximum(pred_train)
minimum(pred_train)

pred_train_int = zeros(Int, length(Y_train))
for i in 1:size(pred_train, 1)
    pred_train_int[i] = findmax(pred_train[i,:])[2]
end

pred_eval_int = zeros(Int, length(Y_eval))
for i in 1:size(pred_eval, 1)
    pred_eval_int[i] = findmax(pred_eval[i,:])[2]
end
sum(pred_train_int .== Y_train), sum(pred_eval_int .== Y_eval)


text:
import HTTP

mutable struct MmapStream # Mmap based
    mmap_ary::Vector{UInt8}
    seekloc::Int
    size::Int
    function MmapStream(filepath::AbstractString)
        size = filesize(filepath)
        new(mmap(filepath), 0, size)
     end
end

read_seek_nb(fobj::MmapStream, seek, nb) = fobj.mmap_ary[seek+1:seek+nb]

function Base.read(fobj::MmapStream, nb::Integer)
    stop = min(fobj.seekloc + nb, fobj.size)
    b = fobj.mmap_ary[fobj.seekloc+1 : stop]
    fobj.seekloc += nb
    return b
end

function Base.close(fobj::MmapStream) # no-op
    nothing
end

# SciToken discovery https://zenodo.org/record/3937438
function _find_scitoken()
    op1 = get(ENV, "BEARER_TOKEN", "")
    op2 = get(ENV, "BEARER_TOKEN_FILE", "")
    op3 = get(ENV, "XDG_RUNTIME_DIR", "")
    uid = @static if Sys.iswindows()
            "julia"
        else
            strip(read(`id -u`, String))
        end
    op3_file = joinpath(op3, "bt_u$uid")
    op4_file = "/tmp/bt_u$uid"
    token = if !isempty(op1)
        op1
    elseif !isempty(op2)
        read(op2, String)
    elseif !isempty(op3) && isfile(op3_file)
        read(op3_file, String)
    elseif isfile(op4_file)
        read(op4_file, String)
    else
        ""
    end
    return strip(token)
end

mutable struct HTTPStream
    uri::HTTP.URI
    seekloc::Int
    size::Int
    multipart::Bool
    scitoken::String
    function HTTPStream(uri::AbstractString; scitoken = _find_scitoken())
        #TODO: determin multipart support
        test = HTTP.request("GET", uri,
        ("Range" => "bytes=0-3", "User-Agent" => "UnROOTjl", "Authorization" => "Bearer $scitoken")
        )
        @assert test.status==206 "bad network or wrong server"
        @assert String(test.body)=="root" "not a root file"
        multipart = false
        local v
        for pair in test.headers
            if lowercase(pair[1]) == "content-range"
                v = pair[2]
                break
            end
        end
        size = parse(Int, match(r"/(\d+)", v).captures[1])
        new(HTTP.URI(uri), 0, size, multipart, scitoken)
    end
end

const SourceStream = Union{MmapStream, HTTPStream}

function Base.read(fobj::SourceStream, ::Type{T}) where T
    return reinterpret(T, read(fobj, sizeof(T)))[1] # TODO: use `only`
end

function Base.position(fobj::SourceStream)
    fobj.seekloc
end

function Base.seek(fobj::SourceStream, loc)
    fobj.seekloc = loc
    return fobj
end

function Base.skip(fobj::SourceStream, stride)
    fobj.seekloc += stride
    return fobj
end

function Base.seekstart(fobj::SourceStream)
    fobj.seekloc = 0
    return fobj
end

function Base.close(fobj::HTTPStream) # no-op
    nothing
end

function Base.read(fobj::HTTPStream, nb::Integer)
    @debug nb
    b = read_seek_nb(fobj, fobj.seekloc, nb)
    fobj.seekloc += nb
    return b
end

function read_seek_nb(fobj::HTTPStream, seek, nb)
    stop = seek+nb-1
    hd = ("Range" => "bytes=$(seek)-$stop", "Authorization" => "Bearer $(fobj.scitoken)")
    b = HTTP.request(HTTP.stack(), "GET", fobj.uri, hd, UInt8[]).body
    return b
end

function Base.read(fobj::SourceStream)
    read(fobj, fobj.size - fobj.seekloc + 1)
end


text:
# So that this test can be run independently
using Cairo
using Base64

include("shape_functions.jl")

@testset "Test that writing images to a Julia IO object works" begin
    c = CairoRGBSurface(256,256)
    cr = CairoContext(c)
    ddots4(cr,256,246,1.0,3000)
    buf = IOBuffer()
    pipe = Base64EncodePipe(buf)
    write_to_png(c,pipe)
    close(pipe)

    # Catch short writes

    str = String(take!(buf))
    str_data = codeunits(str)

    @test length(str_data) > 200
end


text:
Base.IndexStyle(x::SummarizedExperiment) = IndexCartesian()

"""
    getindex(x::SummarizedExperiment, i, j)

Subset `x` by the rows or columns based on `i` and `j`, respectively.
Types for the arguments to `i` and `j` are similar to those for arrays:

- An integer `Vector` containing indices.
- An `Int` containing a single index.
- A boolean `Vector` of length equal to the relevant dimension,
  indicating whether each entry of that dimension should be retained.
- A `:` operator to retain the entirety of a dimension's extent.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10);

julia> x[1,:]
1x10 SummarizedExperiment
  assays(3): foo bar whee
  rownames: Gene1
  rowdata(2): name Type
  colnames: Patient1 Patient2 ... Patient9 Patient10
  coldata(3): name Treatment Response
  metadata(1): version

julia> x[:,1:5]
20x5 SummarizedExperiment
  assays(3): foo bar whee
  rownames: Gene1 Gene2 ... Gene19 Gene20
  rowdata(2): name Type
  colnames: Patient1 Patient2 ... Patient4 Patient5
  coldata(3): name Treatment Response
  metadata(1): version

julia> keep = [ i > 5 for i in 1:size(x)[1] ];

julia> x[keep,1:2]
15x2 SummarizedExperiment
  assays(3): foo bar whee
  rownames: Gene6 Gene7 ... Gene19 Gene20
  rowdata(2): name Type
  colnames: Patient1 Patient2
  coldata(3): name Treatment Response
  metadata(1): version
```
"""
function Base.getindex(x::SummarizedExperiment, i, j)
    indices = to_indices(x, (i, j))

    if isa(indices[1], Int)
        index_rows = [indices[1]]
    else
        index_rows = indices[1]
    end

    if isa(indices[2], Int)
        index_cols = [indices[2]]
    else
        index_cols = indices[2]
    end

    new_rowdata = x.rowdata[index_rows,:]
    new_coldata = x.coldata[index_cols,:]

    new_assays = OrderedDict{String,AbstractArray}();
    for (key, val) in x.assays
        subdex = Array{Any,1}(undef, length(size(val)))
        fill!(subdex, :)
        subdex[1] = index_rows
        subdex[2] = index_cols
        new_assays[key] = getindex(val, subdex...)
    end

    return SummarizedExperiment(new_assays, new_rowdata, new_coldata, copy(x.metadata))
end


text:
using RBNF
using MLStyle
using PrettyPrint

# Using code from https://github.com/thautwarm/RBNF.jl
struct QASMLang end

second((a, b)) = b
second(vec::V) where V <: AbstractArray = vec[2]

struct Struct_bin
    l
    op :: RBNF.Token
    r
end

RBNF.@parser QASMLang begin
    # define ignorances
    ignore{includes, space, comments}

    @grammar
    # define grammars
    mainprogram := ["OPENQASM", ver=nnreal, ';', prog=program]
    program     = statement{*}
    statement   = (decl | gate | opaque | ifstmt | barrier | qop)
    # stmts
    ifstmt      := ["if", '(', l=id, "==", r=nninteger, ')', gate_name=id, ['(', [args=explist].?, ')'].?, outs=mixedlist, ';']
    opaque      := ["opaque", id=id, ['(', [arglist1=idlist].?, ')'].? , arglist2=idlist, ';']
    barrier     := ["barrier", value=mixedlist, ';']
    decl        := [regtype=("qreg" | "creg"), id=id, '[', int=nninteger, ']', ';']

    # gate
    gate        := [decl=gatedecl, [goplist=goplist].?, '}']
    gatedecl    := ["gate", id=id, ['(', [args=idlist].?, ')'].?, (outs=idlist), '{']

    goplist     = (barrier_ids | uop){*}
    barrier_ids := ["barrier", ids=idlist, ';'] # not impl
    # qop
    qop         = (measure | reset | uop)
    reset       := ["reset", arg=argument, ';'] # not impl
    measure     := ["measure", arg1=argument, "->", arg2=argument, ';'] # not impl

    uop         = (u | cx | h | x | y | z | s | sdg | t | tdg | rx | ry | rz | crx | cry | crz | iduop)
    u          := ['U', '(', in1=nnexp, ',', in2=nnexp, ',', in3=nnexp, ')', out=argument, ';']
    cx         := ["CX", out1=argument, ',', out2=argument, ';']
    ch         := ["ch", out1=argument, ',', out2=argument, ';']
    h          := ['h',  out=argument, ';']
    x          := ['x',  out=argument, ';']
    y          := ['y',  out=argument, ';']
    z          := ['z',  out=argument, ';']
    s          := ['s',  out=argument, ';']
    sdg        := ["sdg",  out=argument, ';']
    t          := ['t',  out=argument, ';']
    tdg        := ["tdg",  out=argument, ';']
    rx         := ["rx",  '(', in=nnexp, ')', out=argument, ';']
    ry         := ["ry",  '(', in=nnexp, ')', out=argument, ';']
    rz         := ["rz",  '(', in=nnexp, ')', out=argument, ';']
    crx        := ["crx",  '(', in=nnexp, ')', out1=argument, ',', out2=argument, ';']
    cry        := ["cry",  '(', in=nnexp, ')', out1=argument, ',', out2=argument, ';']
    crz        := ["crz",  '(', in=nnexp, ')', out1=argument, ',', out2=argument, ';']

    iduop      := [gate_name=id, ['(', [args=explist].?, ')'].?, outs=mixedlist, ';']

    idlist     := [hd=id, [',', tl=idlist].?]

    mixedlist  := [hd=argument, [',', tl=mixedlist].?]

    argument   := [id=id, ['[', (arg=nninteger), ']'].?]

    explist    := [hd=nnexp, [',', tl=explist].?]
    pi         := "pi"
    atom       =  (nnreal | nninteger | fnexp | pi | id) | (['(', nnexp, ')'] % second) | neg
    fnexp      := [fn=fn, '(', arg=nnexp, ')']
    neg        := ['-', value=nnexp]

    nnexp        = [l=mul,  [op=('+' |'-'), r=nnexp].?] => _.op === nothing ? _.l : Struct_bin(_.l, _.op, _.r)
    mul        = [l=atom, [op=('*' | '/'), r=mul].?] => _.op === nothing ? _.l : Struct_bin(_.l, _.op, _.r)
    fn         = ("sin" | "cos" | "tan" | "exp" | "ln" | "sqrt")

    # define tokens
    @token
    includes  := r"\Ginclude .*;"
    id        := r"\G[a-z]{1}[A-Za-z0-9_]*"
    nnreal    := r"\G([0-9]+\.[0-9]*|[0-9]*\.[0.9]+)([eE][-+]?[0-9]+)?"
    nninteger := r"\G([1-9]+[0-9]*|0)"
    space     := r"\G\s+"
    comments  := r"\G//.*"
end

Token = RBNF.Token

function lex(src :: String)
    RBNF.runlexer(QASMLang, src)
end

function parse_qasm(tokens :: Vector{Token{A} where A})
    ast, ctx = RBNF.runparser(mainprogram, tokens)
    ast
end