text:
info("Testing ", MOD.CompositionClass)
@test MOD.iscomposable(MLKTest.TestClass(), MLKTest.TestKernel()) == false
for class_obj in composition_classes
    info("Testing ", class_obj)

    # Test constructors
    for T in FloatingPointTypes
        default_floats, default_others = all_default_args[class_obj]
        default_args = (T[default_floats...]..., default_others...)
        fields = fieldnames(class_obj)
        k = (class_obj)(default_args...)
        @test eltype(k) == T

        for i in eachindex(fields)
            @test getfield(k, fields[i]).value == default_args[i]
        end

        for z in (zero(T),one(T))
            f = all_phifunctions[class_obj]
            @test_approx_eq MOD.phi(k, z) f(default_args..., z)
        end
    end

    # Test conversions
    for T in FloatingPointTypes
        phi = (class_obj)()
        for U in FloatingPointTypes
            @test U == eltype(convert(CompositionClass{U}, phi))
        end
    end

    # Test equality
    phi1 = (class_obj)()
    for class_obj2 in composition_classes
        phi2 = (class_obj2)()
        if class_obj == class_obj2
            @test phi1 == phi2
        else
            @test phi1 != phi2
        end
    end

    # Test properties
    properties = all_kernelproperties[class_obj]
    phi = (class_obj)()

    @test MOD.ismercer(phi)        == properties[1]
    @test MOD.isnegdef(phi)        == properties[2]
    @test MOD.attainsnegative(phi) == properties[3]
    @test MOD.attainszero(phi)     == properties[4]
    @test MOD.attainspositive(phi) == properties[5]

    # Test iscomposable() rules
    valid_kernel_objs = composition_pairs[class_obj]
    for kernel_obj in additive_kernels
        isvalid = get(valid_kernel_objs, kernel_obj, false)
        k = (kernel_obj)()
        @test MOD.iscomposable(phi, k) == isvalid
    end

    # Test that output does not create error
    show(DevNull, phi)
end

info("Testing ", MOD.KernelComposition)
for class_obj in composition_classes
    phi = (class_obj)()
    valid_kernel_objs = composition_pairs[class_obj]
    for kernel_obj in additive_kernels
        kappa = (kernel_obj)()
        if get(valid_kernel_objs, kernel_obj, false)
            psi = phi ∘ kappa
            @test MOD.ismercer(psi)        == MOD.ismercer(phi)
            @test MOD.isnegdef(psi)        == MOD.isnegdef(phi)
            @test MOD.attainsnegative(psi) == MOD.attainsnegative(phi)
            @test MOD.attainszero(psi)     == MOD.attainszero(phi)
            @test MOD.attainspositive(psi) == MOD.attainspositive(phi)
            show(DevNull, psi)
        else
            @test_throws ErrorException phi ∘ kappa
        end
    end
end

for kernel_obj in composition_kernels
    @test typeof((kernel_obj)()) <: KernelComposition
end

for kernel_obj in additive_kernels
    kappa = (kernel_obj)()

    if MOD.iscomposable(PolynomialClass(), kappa)
        @test typeof((kappa^2).phi) <: PolynomialClass
    else
        @test_throws ErrorException kappa^2
    end

    if MOD.iscomposable(PowerClass(), kappa)
        @test typeof((kappa^0.5).phi) <: PowerClass
    else
        @test_throws ErrorException kappa^0.5
    end

    if MOD.iscomposable(ExponentiatedClass(), kappa)
        @test typeof((exp(kappa)).phi) <: ExponentiatedClass
    else
        @test_throws ErrorException exp(kappa)
    end

    if MOD.iscomposable(SigmoidClass(), kappa)
        @test typeof((tanh(kappa)).phi) <: SigmoidClass
    else
        @test_throws ErrorException tanh(kappa)
    end

end


text:
# only prints if str is nonempty
function printiferror(log,str)
	str = strip(str)
	isempty(str) || println(log, "ERROR: ", str)
end

# only prints if str is nonempty
function printifwarning(log,str)
	str = strip(str)
	isempty(str) || println(log, "WARNING: ", str)
end

# only prints if str is nonempty
function printifinfo(log,str)
	str = strip(str)
	isempty(str) || println(log, str)
end

# annoying helper functions that are needed since I cannot redirect stdout/err to IOBuffers
# and working with Pipes directly is difficult as they might block, waiting for input,
# and it's hard to know when that will happen.
# Should be replaced by a better solution using IOBuffers or Pipes directly when Julia has better support.

# TODO: fix using redirect_stdout/redirect_stderr ???

mutable struct DiskBuffer
	filename::String
	stream#::IOStream
end
DiskBuffer() = DiskBuffer(tempname(),nothing)

openbuf(db::DiskBuffer) = db.stream = open(db.filename, "w")
function closebuf(db::DiskBuffer)::String
	db.stream == nothing && return ""
	close(db.stream)
	str = read(db.filename, String)
	rm(db.filename) # delete file
	str
end

ext: jl


text:
using Scientific_Programming_in_Julia
using Base: tail

import Scientific_Programming_in_Julia.ReverseDiff

using Base: tail
struct Dense{M<:TrackedMatrix, V<:TrackedVector, F}
    W::M
    b::V
    f::F
end

Dense(W::Matrix, b::Vector, f) = Dense(track(W),track(b),f)
Dense(in,out,f=identity) = Dense(rand(out,in) .-0.5, rand(out) .-0.5, f)
(m::Dense)(x) = m.f(m.W*x + m.b)

struct Chain{T<:Tuple}
    layers::T
end
Chain(ls...) = Chain(ls)

(m::Chain)(x) = applychain(m.layers, x)
applychain(ls::Tuple{}, x) = x
applychain(ls::Tuple, x) = applychain(tail(ls), first(ls)(x))

params(m::Dense) = (W=m.W, b=m.b)
params(m::Chain) = params(m.layers)
params(ls::Tuple) = (params(first(ls))..., params(tail(ls))...)
params(ls::Tuple{}) = ls

function ReverseDiff.reset!(m::Chain)
    map(ReverseDiff.reset!, m.layers)
    nothing
end

function ReverseDiff.reset!(m::Dense)
    ReverseDiff.reset!(m.W)
    ReverseDiff.reset!(m.b)
    nothing
end

function train_step!(loss, network, xs, ys; λ=0.01)
    l = loss(network, xs, ys)
    l.grad = 1.0
    for w in params(network)
        w.data .-= λ .* accum!(w)
    end
    l
end

# task
f(x,y) = y^2 + sin(x)
f(xs) = [f(xs[1],xs[2])]

# data
xs = map(1:30) do i
    x = rand(-4:0.1:4)
    y = rand(-2:0.1:2)
    [x,y]
end
ys = track.(f.(xs))
xs = track.(xs)

hdim = 15
network = Chain(
    Dense(2,hdim,σ),
    Dense(hdim,hdim,σ),
    Dense(hdim,1)
)
function loss(network,xs,ys)
    errs = hcat([abs2(network(x)-y) for (x,y) in zip(xs,ys)]...)
    sum(errs)
end

using Plots
color_scheme = cgrad(:RdYlBu_5, rev=true)

forward(network,x,y) = network(track([x,y])).data[1]

# training
anim = @animate for i in 1:2000
    ReverseDiff.reset!.(xs)
    ReverseDiff.reset!.(ys)
    ReverseDiff.reset!(network)
    l = train_step!(loss, network, xs, ys, λ=0.003)
    if mod1(i,50) == 1
        @info i l
        p1 = contour(-4:0.3:4, -2:0.3:2, f, fill=true, c=color_scheme, xlabel="x", ylabel="y", title="Truth")
        p2 = contour(-4:0.3:4, -2:0.3:2, (x,y)->forward(network,x,y), fill=true, c=color_scheme, xlabel="x",title="Iteration: $i")
        p = plot(p1,p2,size=(1200,400)) |> display
    end
end every 50

gif(anim, "anim.gif", fps=15)
error()

for i in 1:20000
    ReverseDiff.reset!.(xs)
    ReverseDiff.reset!.(ys)
    ReverseDiff.reset!(network)
    l = train_step!(loss, network, xs, ys, λ=0.002)
    if mod1(i,1000) == 1
        @info i l
        p1 = contour(-4:0.3:4, -2:0.3:2, f, fill=true, c=color_scheme, xlabel="x", ylabel="y", title="Truth")
        p2 = contour(-4:0.3:4, -2:0.3:2, (x,y)->forward(network,x,y), fill=true, c=color_scheme, xlabel="x",title="Iteration: $i")
        plot(p1,p2) |> display
    end
end


text:
using Agents: vertices
using MetaGraphs: add_vertex!, has_edge
using Agents, AgentsPlots, Plots, LightGraphs, MetaGraphs, SimpleWeightedGraphs, GraphPlot, GraphRecipes, NetworkLayout
using Tables, DataFrames
using CSV, JSON, Serialization, DelimitedFiles
using BritishNationalGrid
using ZipFile, Shapefile
using Random
using Match
using LinearAlgebra
using StatsBase
using Distributions
using StatsPlots
using SparseArrays
using Laplacians
using DataStructures
using RollingFunctions
using BenchmarkTools
using Statistics
using Logging,LoggingExtras, LoggingFacilities
using Dates,TimeZones

csv_custom_nodes = "/home/mep53/workspace/data/bt/metro-access/"*"nodes_t2_t3_150618_CutA.csv"
csv_custom_links = "/home/mep53/workspace/data/bt/metro-access/"*"New21CNCoreTopology.csv"
shp_geo_uk = "/home/mep53/workspace/uk/geo/"*"NUTS_UK_2018.shp"

include("netManAbm.jl")
loggers = Dict()

simplified_logger(logger) = TransformerLogger(logger) do log
    log.
    merge(log, (; message = "$(log.message)"
    , file="", line="", _module=""
    ))
end

function single_run_with_logging(config)
    run_label = get_run_label(config)
    io = stdout
    logger = MessageOnlyLogger(io,Logging.Info)
    with_logger(logger) do
        start_time = now()
        log_info("$start_time: start $run_label")
        single_run(config)
        end_time = now()
        log_info("$end_time: end $run_label. Elapsed: $((end_time - start_time))")
    end
end

function single_run_with_file_logging(config)

    run_label = get_run_label(config)

    io = open( data_dir * run_label * "_log.log", "w+")
    # logger = SimpleLogger(io,Logging.Debug)

    # logger = FormatLogger() do io, args
    #     log_info(io, args._module, " | ", "[", args.level, "] ", args.message)
    # end;
    # logger = MinLevelLogger(FileLogger(run_label* "_testing.log"), Logging.Info) |> simplified_logger

    # logger = OneLineTransformerLogger(MinLevelLogger(FileLogger( data_dir * run_label* ".log"), Logging.Info)#|> OneLineTransformerLogger
    # logger = SimpleLogger(stdout, Logging.Debug) |> OneLineTransformerLogger
    logger = MessageOnlyLogger(io,Logging.Info)
    with_logger(logger) do
        start_time = now()
        log_info("$start_time: start $run_label")
        single_run(config)
        end_time = now()
        log_info("$end_time: end $run_label. Elapsed: $((end_time - start_time))")
    end
    flush(io)
    close(io)
end

data_dir = "data/"
plots_dir = "plots/"
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100
#configs = load_run_configs()

services = [(1,7),(4,1),(5,14),(12,8)]

scenario_2_ruls = [65,78,84,86,93,49,90]
scenario_2_ids = [2,6,9,13,10,15,16]

steps = 145
configs = []
#push!(configs,first(load_run_configs(16,0,services,steps)))
#push!(configs,first(load_run_configs(16,1,services,steps)))
push!(configs,first(load_run_configs(16,2,services,steps)))

#Logging.disable_logging(Logging.Info)
#enable logs
Logging.disable_logging(Logging.BelowMinLevel)
Threads.@threads for config in configs
    single_run_with_logging(config)
    #single_run_with_file_logging(config)
end


text:
# This file includes the structure of Branch-and-Bound tree.

include_dependency("../dd/core.jl")
include_dependency("../dd/ordering.jl")
include_dependency("../dd/mathprog/dp_specification.jl")
include_dependency("../dd/mathprog/function_analysis.jl")

"""
Type to store node information in b&b
"""
struct NodeInfo
    parent_id::Int                  # id of the parent node in the tree (array)
    children_id::Tuple{Int, Int}      # id of the children nodes in the tree
    depth::Int                         # depth of the node in the tree
end

"""
Type to store model information in b&b
"""
struct ModelInfo
    m::JuMP.Model                       # the model after adding cuts in the node
    #    best_primal_sol::Array{Float64, 1}          # the best feasible solution at this node
    #    obj_primal::Float64                         # objective lower (primal) bound obtained at the node
    #    best_dual_sol::Array{Float64, 1}           # the optimal solution of the relaxation at this node
    #    obj_dual::Float64                         # objective upper (dual) bound obtained at the node
end

"""
Type to store DD information in b&b
"""
struct DDInfo
    dd::Array{DecisionDiagram,1}        # array of dds of all constraints in the model
    initial::Array{Float64, 1}              # array of intial states for constraints
    transition::Array{Array{UVFunc, 1}, 1}     # 1st dim: constraints, 2nd dim: variables, contenct: transtiion function
    con_function::Array{Function, 1}            # stores the function representing variable part of constraints
    domain::Array{Tuple{Float64, Float64, Symbol}, 1}   # array of variable domains
    order::Array{Ordering, 1}                   # ordering of variables in dd, for constraints
    lb_ub::Array{Array{Tuple{Float64, Float64}, 1},1}   # 1st dim: constraint, 2nd dim: variable, content: (lowerbound, upperbound) of the univariate function over the variable range defined in specs.
end

"""
Type to store Branching information in b&b
"""
struct BranchInfo
    #    branch_policy::?                   # which of the two policies is used in this branch form the parent node: variable range (arc label) branch, constraint value (node state) branch, and on which variables or constraints?
    #    refine_policy::Array{Bool, 1}        # do we reconstruct the dds in this branch (having more width) or we just leave it as it is after branching
    #    modified_node_layers::Array{Array{Array{Bool, 1}, 1}, 1}        # 1st dim: constraints, 2nd dim and 3rd dim: dd node layers, content: whether the node is removed (0) or not from the dd of the parent node
    #    modified_arc_layers::Array{Array{Array{Bool, 1}, 1}, 1}     # similar to the modified_node_layers, but for arcs
end

"""
Type to store branch-and-bound nodes
"""
struct BBNode
    node_info::NodeInfo
    model_info::ModelInfo
    dd_info::DDInfo
    branch_info::BranchInfo
end

"""
Type to store branch-and-bound tree
"""
struct BBTree
    original_model::JuMP.Model      # the main model of the problem
    LP_model::JuMP.Model            # the LP model used in the ECP
    tree::Array{BBNode, 1}          # an array that stores branch-and-bound nodes (works as a tree)
end

"""
Create a BB tree from a model

# Input
- `m1`: Original model
- `m2`: LP relaxation of the model used in the ECP
"""
function create_BB_tree(m1::JuMP.Model, m2::JuMP.Model)::BBTree
    bb_tree = BBTree(m1, m2, Array{BBNode, 1}())
    return bb_tree
end

"""
create root node in branch-and-bound

# Input
- `bb_tree`: the newly constructed object of type BBTree, which contains a model and its tree array field must be empty.
- `prop`: Keyword that contains an initial assessment for properties of the univariate functions (optional)
"""
function create_BB_root!(bb_tree::BBTree; prop::Array{Array{Array{UVFuncProp, 1}, 1}, 1} = [])
    @assert(length(bb_tree.tree) == 0)          # ensure that the tree array field must be empty

    n::Int = MathProgBase.numvar(bb_tree.original_model)      # number of variables in model
    c_n::Int = JuMP.numnlconstr(bb_tree.original_model) + MathProgBase.numquadconstr(bb_tree.original_model)      # number of nonlinear constraints in model

    node_info = NodeInfo(0, (0,0), 0)
    branch_info = BranchInfo()
    model_info = ModelInfo(copy(bb_tree.LP_model))
    dd_info = DDInfo(Array{DecisionDiagram,1}(  c_n), Array{Float64, 1}(  c_n), Array{Array{UVFunc, 1}, 1}(  c_n), Array{Function, 1}(  c_n), Array{Tuple{Float64, Float64, Symbol}, 1}(  n), Array{Ordering, 1}(  c_n), Array{Array{Tuple{Float64, Float64}, 1},1}(  c_n))

    # Determine constraint decomposition, ordering and variable domains
    eval = evaluate_separable_constraints(bb_tree.original_model)
    create_constraint_domain_function!(dd_info.domain, bb_tree.original_model)
    for i=1:c_n
        dd_info.order[i] = NoOrdering()
    end

    # Create specs for decision diagram
    for i::Int=1:c_n
        dd_info.initial[i] = create_constraint_initial_state(eval[i])
        dd_info.con_function[i] = create_constraint_single_function(eval[i])

        dd_info.transition[i] = Array{UVFunc, 1}(n)
        if isassigned(prop, i)
            create_constraint_transition_function!(dd_info.transition[i], eval[i], prp = prop[i])            # if functional property is given as input
        else
            create_constraint_transition_function!(dd_info.transition[i], eval[i])                 # no functional properties is given as input
        end

        dd_info.lb_ub[i] = Array{Tuple{Float64, Float64}, 1}(  n)
        for j::Int=1:n
            if isassigned(dd_info.transition[i], j)
                dd_info.lb_ub[i][j] = (uvf_lb(dd_info.transition[i][j].f, dd_info.transition[i][j].prop, dd_info.domain[j]), uvf_ub(dd_info.transition[i][j].f, dd_info.transition[i][j].prop, dd_info.domain[j]))
            end
        end

        println("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
        # Construct decision diagram
        dd_info.dd[i] = construct_DD(dd_info.transition[i], dd_info.initial[i], dd_info.domain, dd_info.lb_ub[i], ordering=dd_info.order[i], reduced_arc=true)
        println("&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&")
    end

    push!(bb_tree.tree, BBNode(node_info, model_info, dd_info, branch_info))

end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
# eliminate children
function collapse(ht::HTensorNode{T}) where {T}
   if isleaf(ht); return ht; end
   core = coreshape(ht)
   nchildren = length(ht.children)
   newsubshape = Int[]
   newchildren = HTensorNode{T}[]
   for i in 1:nchildren
      child = ht.children[i]
      core = tensordot(core, child.data, i)
      append!(newsubshape, child.subshape)
      append!(newchildren, child.children)
   end
  #println(STDERR, "!!! collapse: $(ht.subshape) -> $(newsubshape)")
   core = reshape(core, newsubshape..., ht.nbasis)
   HTensorNode(core, newchildren)
end

# eliminate children recursively
function collapse_all(ht::HTensorNode)
   if isleaf(ht); return ht; end
   newchildren = map(collapse_all, ht.children)
   collapse(HTensorNode(ht.data, newchildren))
end


text:
function RegressionResults(arg0::Vector{jdouble}, arg1::Vector{Vector{jdouble}}, arg2::jboolean, arg3::jlong, arg4::jint, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jboolean, arg9::jboolean)
    return RegressionResults((Vector{jdouble}, Vector{Vector{jdouble}}, jboolean, jlong, jint, jdouble, jdouble, jdouble, jboolean, jboolean), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_adjusted_r_squared(obj::RegressionResults)
    return jcall(obj, "getAdjustedRSquared", jdouble, ())
end

function get_covariance_of_parameters(obj::RegressionResults, arg0::jint, arg1::jint)
    return jcall(obj, "getCovarianceOfParameters", jdouble, (jint, jint), arg0, arg1)
end

function get_error_sum_squares(obj::RegressionResults)
    return jcall(obj, "getErrorSumSquares", jdouble, ())
end

function get_mean_square_error(obj::RegressionResults)
    return jcall(obj, "getMeanSquareError", jdouble, ())
end

function get_n(obj::RegressionResults)
    return jcall(obj, "getN", jlong, ())
end

function get_number_of_parameters(obj::RegressionResults)
    return jcall(obj, "getNumberOfParameters", jint, ())
end

function get_parameter_estimate(obj::RegressionResults, arg0::jint)
    return jcall(obj, "getParameterEstimate", jdouble, (jint,), arg0)
end

function get_parameter_estimates(obj::RegressionResults)
    return jcall(obj, "getParameterEstimates", Vector{jdouble}, ())
end

function get_r_squared(obj::RegressionResults)
    return jcall(obj, "getRSquared", jdouble, ())
end

function get_regression_sum_squares(obj::RegressionResults)
    return jcall(obj, "getRegressionSumSquares", jdouble, ())
end

function get_std_error_of_estimate(obj::RegressionResults, arg0::jint)
    return jcall(obj, "getStdErrorOfEstimate", jdouble, (jint,), arg0)
end

function get_std_error_of_estimates(obj::RegressionResults)
    return jcall(obj, "getStdErrorOfEstimates", Vector{jdouble}, ())
end

function get_total_sum_squares(obj::RegressionResults)
    return jcall(obj, "getTotalSumSquares", jdouble, ())
end

function has_intercept(obj::RegressionResults)
    return jcall(obj, "hasIntercept", jboolean, ())
end


text:
function fill_wavelets!(iss::Int64, wavelets::Array{Array{Float64,2},2}, srcwav::Array{SrcWav}, sflags::Vector{Int64})

	npw = size(wavelets,1)
	nt = size(wavelets,2)
	δt = step(srcwav[1][1].grid)
	for ipw=1:npw
		ns, snfield = size(wavelets[ipw,1]) # ns may vary with ipw
		for ifield=1:snfield, is=1:ns
			snt = length(srcwav[ipw][1].grid)
			if(sflags[ipw] == 0)
				# nothing # just put zeros, no sources added
				for it=1:nt
					wavelets[ipw,it][is,ifield] = 0.0
				end
			elseif(sflags[ipw] == 1)
				"ϕ[t] = s[t]"
				for it=1:snt
					source_term = srcwav[ipw][iss].d[ifield][it,is]
					wavelets[ipw,it][is,ifield] = source_term
				end
			elseif(sflags[ipw] == -1)
				"source sink for 1"
				"ϕ[t] = -s[t]"
				for it=2:snt
					source_term = srcwav[ipw][iss].d[ifield][it,is]
					wavelets[ipw,nt-it+2][is,ifield] = -1.0*source_term
				end
			elseif(sflags[ipw] == 2)
				"ϕ[t] = ∑₁ᵗ⁻¹ s[t]"
				source_term_stack = 0.0;
				if(ifield == 1)
					for it=1:snt-1
						source_term_stack += (srcwav[ipw][iss].d[ifield][it,is] .* δt)
						wavelets[ipw,it+1][is,ifield] = source_term_stack
					end
				else
					for it=2:snt-1
						source_term_stack += (((srcwav[ipw][iss].d[ifield][it,is] .* δt) +
						   (srcwav[ipw][iss].d[ifield][it-1,is] .* δt)) * 0.5)
						wavelets[ipw,it+1][is,ifield] = source_term_stack
					end

				end
				if(nt > snt)
					wavelets[ipw,snt+1:end][is,ifield] = wavelets[ipw,snt][is,ifield]
				end
			elseif(sflags[ipw] == -2)
				"source sink for 2"
				# multiplication with -1 for subtraction #time reversal
				# as the source wavelet has to be subtracted before the propagation step, I shift here by one sample"
				source_term_stack = 0.0;
				for it=1:snt-1
					source_term_stack += (srcwav[ipw][iss].d[ifield][it,is] .* δt)
					wavelets[ipw,nt-it+1][is,ifield] = -1.0 * source_term_stack
				end
				if(nt > snt)
					nt_diff = nt-snt
					wavelets[ipw,1:nt_diff+1][is,ifield] = wavelets[ipw,nt_diff+2][is,ifield]
				end
			end
		end
	end
end

struct Source_B1 end
struct Source_B0 end

# This routine ABSOLUTELY should not allocate any memory, called inside time loop.
@inbounds @fastmath function add_source!(it::Int64, issp::Int64, iss::Int64, pac::PFdtdc, pass::Vector{PFdtdss}, pap::PFdtdp, ::Source_B1)
	# aliases
	p=pap.p;
	wavelets=pass[issp].wavelets
	ageom=pac.ageom
	isx1=pass[issp].isx1
	isx2=pass[issp].isx2
	isz1=pass[issp].isz1
	isz2=pass[issp].isz2
	ssprayw=pass[issp].ssprayw
	modttI=pac.modttI
	modrr=pac.modrrvz
	"""
	adding source to pressure field at [it]
	"""
	for ipw in pac.activepw
		sfields=pac.sfields[ipw]
		isfields=pac.isfields[ipw]
	if(pac.sflags[ipw] ≠ 0) # only if sflags is non-zero
		pw=p[ipw]
		for (iff, ifield) in enumerate(isfields)
		@simd for is = 1:ageom[ipw][iss].ns
			"""
			use wavelets at [it], i.e., sum of source terms
			until [it-1]
			division of source term with δx and δz (see Jan's fdelmodc manual)
			"""
			source_term = wavelets[ipw,it][is, iff] * pac.δt * pac.δxI * pac.δzI

			"""
			multiplication with modttI
			"""
			pw[isz1[ipw][is], isx1[ipw][is],ifield] +=
				source(source_term,ssprayw[ipw][1,is], pac, isz1[ipw][is], isx1[ipw][is],eval(sfields[iff])())
			pw[isz1[ipw][is], isx2[ipw][is],ifield] +=
				source(source_term,ssprayw[ipw][2,is], pac, isz1[ipw][is], isx2[ipw][is],eval(sfields[iff])())
			pw[isz2[ipw][is], isx1[ipw][is],ifield] +=
				source(source_term,ssprayw[ipw][3,is], pac, isz2[ipw][is], isx1[ipw][is],eval(sfields[iff])())
			pw[isz2[ipw][is], isx2[ipw][is],ifield] +=
				source(source_term,ssprayw[ipw][4,is], pac, isz2[ipw][is], isx2[ipw][is],eval(sfields[iff])())
		end
		end
	end
	end
end

# This routine ABSOLUTELY should not allocate any memory, called inside time loop.
@inbounds @fastmath function add_source!(it::Int64, issp::Int64, iss::Int64, pac::PFdtdc, pass::Vector{PFdtdss}, pap::PFdtdp, ::Source_B0)
	# aliases
	p=pap.p;
	wavelets=pass[issp].wavelets
	ageom=pac.ageom
	isx1=pass[issp].isx1
	isz1=pass[issp].isz1
	ssprayw=pass[issp].ssprayw
	modttI=pac.modttI
	"""
	adding source to pressure field at [it]
	"""
	for ipw in pac.activepw
		sfields=pac.sfields[ipw]
		isfields=pac.isfields[ipw]
	if(pac.sflags[ipw] ≠ 0) # only if sflags is non-zero
		pw=p[ipw]
		for (iff, ifield) in enumerate(isfields)
		@simd for is = 1:ageom[ipw].ns[iss]
			"""
			use wavelets at [it], i.e., sum of source terms
			until [it-1]
			division of source term with δx and δz (see Jan's fdelmodc manual)
			"""
			source_term = wavelets[ipw,it][is, iff] * pac.δt * pac.δxI * pac.δzI

			"""
			multiplication with modttI
			"""
			pw[isz1[ipw][is], isx1[ipw][is],ifield] += source(source_term,1.0,pac,isz1[ipw][is],isx1[ipw][is],eval(sfields[iff])())
		end
	end
	end
	end
end

# on pressure grid
source(source_term, spray, pac, iz, ix, ::P) = source_term * spray * pac.modttI[iz,ix]

# on Vx grid
source(source_term, spray, pac, iz, ix, ::Vx) = source_term * spray * pac.modrrvx[iz,ix]

# on Vz grid
source(source_term, spray, pac, iz, ix, ::Vz) = source_term * spray * pac.modrrvz[iz,ix]


text:
# Copyright 2019, Patricia Levi, Elias Kueppper, and contributors

## dependencies
using StatsPlots, DataFrames, CSV

"""
    groupedbar(var::OptVariable{T,N};
                    title::String="",
                    export_name="",
                    kwargs...) where {T,N}
Convert the provided OptVariable into a groupedbar-plot and export it
"""
import StatsPlots.groupedbar
function groupedbar(var::OptVariable{T,N};
                    title::String="",
                    export_name="",
                    kwargs...) where {T,N}
    if N==2
        outplot = groupedbar(var.axes[1], var.data; labels=var.axes[2], title = title, kwargs...)
    end
    if export_name!=""
        savefig(outplot, export_name*".png")
    end
    return outplot
end

"""
        dataframe(var::OptVariable{T,N};
                    export_name="") where {T,N}
Convert the provided OptVariable into a `DataFrame` and export it
"""
function dataframe(var::OptVariable{T,N};
                    export_name="") where {T,N}
    if N==2
        table=DataFrame(var.data)
        names!(table,Symbol.(var.axes[2]))
        table[Symbol(var.axes_names[1])]=var.axes[1]
    end
    showall(table)
    if export_name!=""
        CSV.write(export_name*".csv",table)
    end
end

"""
        plotcapacity(results::Array{OptResult,1},
                    outputname::String;
                    kwargs...)
Plot the capacity of the `results`-Array and form a table
"""
function plotcapacity(results::Array{OptResult,1},
                    outputname::String;
                    kwargs...)
    desc=Array{String,1}()
    techs=Array{String,1}()
    for result in results
        push!(desc, result.config["descriptor"])
        unique!(union!(techs, result.sets["tech"]["all"]))
    end

    cap=OptVariable{Number}(undef, techs, desc; axes_names=["tech", "results"])
    cap.=0
    for result in results
        for tech in result.sets["tech"]["node"]
            cap[tech, result.config["descriptor"]]=sum(result.variables["CAP"][tech,:,:])
        end
        if haskey(result.sets["tech"], "line")
            lines=load_cep_data_provided(result.config["region"]).lines
            for tech in result.sets["tech"]["line"]
                c=0
                for line in axes(lines, "line")
                    c+=sum(result.variables["TRANS"][tech,:,line])*lines[tech,line].length
                end
                cap[tech, result.config["descriptor"]]=c./1e3
            end
        end
    end
    dataframe(cap; export_name=outputname*"capdata")
    return groupedbar(cap; title="Installed Capacities", export_name=outputname*"capdata", yaxis="Capacities [MW] and Transmission [GW*km]", kwargs...)
end

"""
        plotcost(results::Array{OptResult,1},
                    outputname::String;
                    kwargs...)
Plot the costs of the `results`-Array and form a table
"""
function plotcost(results::Array{OptResult,1},
                outputname::String;
                kwargs...)
    desc=Array{String,1}()
    for result in results
        push!(desc, result.config["descriptor"])
    end

    cost=OptVariable{Number}(undef, desc, ["var", "cap_fix"]; axes_names=["results", "cost_type"])
    for result in results
        for cost_type in axes(cost, "cost_type")
            cost[result.config["descriptor"], cost_type]=sum(result.variables["COST"][cost_type,first(result.sets["impact"]["mon"]),:])
        end
    end
    dataframe(cost; export_name=outputname*"costdata")
    return groupedbar(cost; title="Costs", export_name=outputname*"costdata", bar_position=:stack, yaxis="Cost [$(first(first(results).sets["impact"]["mon"]))]", kwargs...)
end

"""
        plotgen(results::Array{OptResult,1},
                    outputname::String;
                    kwargs...)
Plot the gen of the `results`-Array and form a table
"""
function plotgen(results::Array{OptResult,1},
                outputname::String;
                kwargs...)
    desc=Array{String,1}()
    techs=Array{String,1}()
    for result in results
        push!(desc, result.config["descriptor"])
        unique!(union!(techs, result.sets["tech"]["power"]))
    end

    gen=OptVariable{Number}(undef, techs, desc; axes_names=["results", "tech"])
    gen.=0
    for result in results
        for tech in result.sets["tech"]["power"]
            g=0
            for k in result.sets["time_K"]["all"]
                for t in result.sets["time_T_period"]["all"]
                    g+=sum(result.variables["GEN"][tech,"electricity",t,k,:])*result.config["time_series"]["delta_t"][t,k]*result.config["time_series"]["delta_t"][k]
                end
            end
            gen[tech, result.config["descriptor"]]=g
        end
    end
    dataframe(gen; export_name=outputname*"gendata")
    return groupedbar(gen; title="Total Generation, Demand, Loss", export_name=outputname*"gendata", yaxis="Energy [MWh]", kwargs...)
end


text:
using Dates

struct JTestResult
    result::Bool
    file::String
    line::Int
    expr::String
    time::Dates.Nanosecond
end

elapsed_seconds(j::JTestResult)::Float64 = Float64(j.time.value) / Float64(1e9) # nanoseconds in a second

macro jtest(args...)
    f = string(__source__.file)
    l = Int(__source__.line)
    quote
        local elapsedtime = Dates.Nanosecond(time_ns())
        local val = $(esc(args...))
        elapsedtime = Dates.Nanosecond(time_ns()) - elapsedtime
        JTestResult(val, $f, $l, string($args[1]), elapsedtime)
    end
end


text:
# 2bit Reader
# ===========

mutable struct Reader{T<:IO} <: BioGenerics.IO.AbstractReader
    # input stream
    input::T
    # sequence names
    names::Vector{String}
    # file offsets
    offsets::Vector{UInt32}
    # byte-swapped or not
    swap::Bool
end

"""
    TwoBit.Reader(input::IO)

Create a data reader of the 2bit file format.

# Arguments
* `input`: data source
"""
function Reader(input::IO)
    seqcount, swap = readheader!(input)
    names, offsets = readindex!(input, seqcount, swap)
    @assert seqcount == length(names) == length(offsets)
    return Reader(input, names, offsets, swap)
end

function BioGenerics.IO.stream(reader::Reader)
    return reader.input
end

function Base.eltype(::Type{Reader{T}}) where {T}
    return Record
end

function Base.length(reader::Reader)
    return length(reader.names)
end

function Base.iterate(reader::Reader, i::Int=1)
    if i > lastindex(reader.names)
        return nothing
    else
        return reader[i], i + 1
    end
end

"""
    seqnames(reader::Reader)::Vector{String}

Get the sequence names.

Sequences are stored in this order.
"""
function seqnames(reader::Reader)
    return copy(reader.names)
end

# Random access
# -------------

function Base.checkbounds(reader::Reader, i::Integer)
    if 1 ≤ i ≤ lastindex(reader)
        return true
    end
    throw(BoundsError(reader, i))
end

function Base.checkbounds(reader::Reader, r::AbstractRange)
    if isempty(r) || (1 ≤ first(r) && last(r) ≤ lastindex(reader))
        return true
    end
    throw(BoundsError(reader, r))
end

function Base.lastindex(reader::Reader)
    return length(reader)
end

function Base.getindex(reader::Reader, i::Integer)
    checkbounds(reader, i)
    seek(reader.input, reader.offsets[i])
    record = Record()
    read!(reader, record)
    return record
end

function Base.getindex(reader::Reader, r::AbstractRange)
    checkbounds(reader, r)
    return [reader[i] for i in r]
end

function Base.getindex(reader::Reader, name::AbstractString)
    i = findfirst(isequal(name), reader.names)
    if i === nothing
        throw(KeyError(name))
    end
    return reader[i]
end

function Base.getindex(reader::Reader, names::AbstractVector{S}) where S<:AbstractString
    return [reader[name] for name in names]
end

# Readers
# -------

const SIGNATURE = UInt32(0x1A412743)

function readheader!(input)
    signature = read(input, UInt32)
    if signature == SIGNATURE
        swap = false
    elseif signature == bswap(SIGNATURE)
        swap = true
    else
        error("invalid 2bit file signature")
    end
    version  = read32(input, swap)
    seqcount = read32(input, swap)
    reserved = read32(input, swap)
    @assert version == reserved == 0
    return seqcount, swap
end

function readindex!(input, seqcount, swap)
    names = String[]
    offsets = UInt32[]
    for i in 1:seqcount
        namesize = read(input, UInt8)
        name = read(input, namesize)
        offset = read32(input, swap)
        push!(names, String(name))
        push!(offsets, offset)
    end
    return names, offsets
end

function Base.read!(reader::Reader, record::Record)
    initialize!(record)
    record.dnasize = read32(reader)
    record.blockcount = read32(reader)
    read32!(reader, record.blockstarts, record.blockcount)
    read32!(reader, record.blocksizes, record.blockcount)
    record.maskedblockcount = read32(reader)
    read32!(reader, record.maskedblockstarts, record.maskedblockcount)
    read32!(reader, record.maskedblocksizes, record.maskedblockcount)
    record.reserved = read32(reader)
    datalen = cld(record.dnasize, 4)
    if length(record.packeddna) < datalen
        resize!(record.packeddna, datalen)
    end
    unsafe_read(reader.input, pointer(record.packeddna), datalen)
    record.filled = true
    return record
end

# read 32 bits and swap bytes if necessary
function read32(reader::Reader)
    return read32(reader.input, reader.swap)
end

function read32!(reader::Reader, out, n)
    resize!(out, n)
    read!(reader.input, out)
    if reader.swap
        map!(bswap, out)
    end
    return out
end

function read32(input::IO, swap)
    x = read(input, UInt32)
    if swap
        x = bswap(x)
    end
    return x
end


text:
using Test, YaoBlocks, YaoBase, LuxurySparse, YaoArrayRegister

@test_throws LocationConflictError swap(2, 1, 1)
@test swap(2, 1, 2) isa Swap
@test mat(swap(2, 1, 2)) ≈ PermMatrix([1, 3, 2, 4], ones(1 << 2))
@test mat(swap(4, 1, 3)) ≈ linop2dense(s -> instruct!(Val(2), s, Val(:SWAP), (1, 3)), 4)

@test swap(1, 2)(3) == swap(3, 1, 2)

@testset "check apply" begin
    r = rand_state(2)
    @test apply!(copy(r), YaoBlocks.ConstGate.SWAP) ≈ apply!(copy(r), swap(2, 1, 2))
end

@test occupied_locs(swap(3, 1, 2)) == (1, 2)
@test isunitary(swap(3, 1, 2))
@test ishermitian(swap(3, 1, 2))
@test isreflexive(swap(3, 1, 2))

@testset "pswap gate" begin
    pb = pswap(6, 2, 4, 0.0)
    @test pb isa PSwap{Float64} && pb.n == 6
    @test pb == pswap(2, 4, 0.0)(6)
    reg = rand_state(6)
    @test copy(reg) |> pb ≈ invoke(apply!, Tuple{ArrayReg,PutBlock}, copy(reg), pb)
    @test copy(reg) |> pb ≈ reg

    dispatch!(pb, π)
    @test copy(reg) |> pb ≈ -im * (copy(reg) |> swap(6, 2, 4))
    @test copy(reg) |> pb |> isnormalized
    pb = dispatch(pb, π)
    @test copy(reg) |> pb ≈ -im * (copy(reg) |> swap(6, 2, 4))
    @test copy(reg) |> pb |> isnormalized

    pb = pswap(6, 2, 4, 0.0)
    dispatch!(pb, :random)
    @test copy(reg) |> pb ≈ invoke(apply!, Tuple{ArrayReg,PutBlock}, copy(reg), pb)
    pb = dispatch(pb, :random)
    @test copy(reg) |> pb ≈ invoke(apply!, Tuple{ArrayReg,PutBlock}, copy(reg), pb)
end


text:
#=
author: mingchuno
smart permutation inspired by python-high-perf
=#
array = [0 ,2 ,3, 4, 5, 6, 7 ,8, 9]
sarray = Set(array)

for last4_arr in combinations(array, 4), efgh_arr in permutations(last4_arr)
  e = efgh_arr[1]
  if e == 0 continue end # fail fast

  f = efgh_arr[2]

  g = efgh_arr[3]
  if g == 0 continue end # fail fast

  h = efgh_arr[4]

  ef = e*10 + f
  gh = g*10 + h

  if ef + gh == 111
    remain_arr = collect(setdiff(sarray, Set(last4_arr)))
    for first4_arr in combinations(remain_arr, 4), abcd_arr in permutations(first4_arr)
      a = abcd_arr[1]
      if a == 0 continue end # fail fast

      b = abcd_arr[2]

      c = abcd_arr[3]
      if c == 0 continue end # fail fast

      d = abcd_arr[4]
      ab = a*10 + b
      cd = c*10 + d
      if ab - cd == ef
        println("($(ab) - $(cd) = $(ef)) + $(gh) = 111")
      end
    end
  end
end

max_stars_repo_name: jackhftang/ABCDEFGHPPP

max_issues_repo_name: jackhftang/ABCDEFGHPPP

max_forks_repo_name: jackhftang/ABCDEFGHPPP


text:
using Multilane
using StatsBase

solvers = Dict{String, Any}()

dpws = DPWSolver(depth=20,
                 n_iterations=100,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())

solvers["dpw_100"] = deepcopy(dpws)
dpws.n_iterations = 300
solvers["dpw_300"] = deepcopy(dpws)
dpws.n_iterations = 500
solvers["dpw_500"] = deepcopy(dpws)

curve = TestSet(lambda=[1., 10., 17.78, 31.62, 56.23, 100.], N=500)

tests = []
for (s,_) in solvers
    push!(tests, TestSet(curve, solver_key=s))
end

objects = gen_initials(tests)

@show objects["param_table"]
objects["solvers"] = solvers

results = evaluate(tests, objects)


text:
@testset "Neighborhoods" begin
  @testset "BallNeighborhood" begin
    # Euclidean metric
    n = BallNeighborhood(0.5)
    @test isneighbor(n, [0.], [0.])
    @test !isneighbor(n, [0.], [1.])

    n = BallNeighborhood(1.0)
    @test isneighbor(n, [0.,0.], [0.,0.])
    @test isneighbor(n, [0.,0.], [1.,0.])
    @test isneighbor(n, [0.,0.], [0.,1.])

    # Chebyshev metric
    n = BallNeighborhood(0.5, Chebyshev())
    @test isneighbor(n, [0.], [0.])
    @test !isneighbor(n, [0.], [1.])

    for r in [1.,2.,3.,4.,5.]
      n = BallNeighborhood(r, Chebyshev())
      for i in 0.0:1.0:r, j in 0.0:1.0:r
        @test isneighbor(n, [0.,0.], [i,j])
      end
    end

    n = BallNeighborhood(1.0)
    @test sprint(show, n) == "BallNeighborhood(1.0, Euclidean(0.0))"
    @test sprint(show, MIME"text/plain"(), n) == "BallNeighborhood\n  radius: 1.0\n  metric: Euclidean(0.0)"
  end
end


text:
struct ARParams <: AbstractParams
    c::Bool
    p::Vector{<:Integer}
end

ARParams(c::Bool, p::Integer) = ARParams(c, collect(1:p))
ARParams(p::Integer) = ARParams(true, p)

function fit(params::ARParams, z::AbstractVector{T}) where T
    p = params.p
    P = maximum(p)
    if p == collect(1:P)
        ϕ, σ2 = levinson_durbin(z, P)
        μ = params.c ? mean(z) * (1 - sum(ϕ)) : zero(T)
    else
        μ, ϕ, θ, β, σ2 = least_squares(z, params.c, params.p)
    end
    return ARModel{P}(μ, ϕ, σ2)
end

const ARModel{p,T} = ARMAXModel{p,0,0,T}
ARModel{p}(μ::T,ϕ::AbstractVector{T},σ2::T) where {p,T} = ARMAXModel{p,0,0,T}(μ, ϕ, T[], T[], σ2)

isstationary(model::ARModel) = all(norm.(roots(Polynomial([1;.-model.ϕ]))) .> 1)

isinvertible(model::ARModel) = true

function forecast(model::M, z::AbstractVector{T}) where {p,T,M <: ARModel{p,T}}
    N = length(z)
    zhat = model.μ
    for i in 1:min(N, p)
        zhat += model.ϕ[i] * z[N - i + 1]
    end
    return zhat
end

forecast(model::M, z::AbstractVector{T}, a::AbstractVector{T}) where {p,T,M <: ARModel{p,T}} = forecast(model, z)

least_squares(z::AbstractVector{T}, c::Bool, p::Vector{<:Integer}) where T = least_squares(z, T[], zeros(T, 0, 0), c, p, Int[])
least_squares(z::AbstractVector{T}, c::Bool, p::Integer) where T = least_squares(z, c, collect(1:p))
least_squares(z::AbstractVector, p::Integer) = least_squares(z, true, p)

function yule_walker(z::AbstractVector, p::Integer)
    ρ = map(i -> autocorrelation(z, i), 1:p)
    Γ = autocorrelation_matrix(z, p)
    ϕ = cholesky(Γ) \ ρ
    m = mean(z)
    γ0 = autocovariance(z, 0)
    σ2 = (1 - dot(ϕ, ρ)) * γ0
    return ϕ, σ2
end

function levinson_durbin!(ϕ::AbstractVector{T}, σ2::Ref{T}, ρ::AbstractVector{T}, p::Integer) where {T}
    num = ρ[p]
    den = one(T)
    for j in 1:p - 1
        num -= ϕ[j] * ρ[p - j]
        den -= ϕ[j] * ρ[j]
    end
    ϕ[p] = num / den
    for j in 1:div(p - 1, 2)
        tmp = ϕ[j]
        ϕ[j] -= ϕ[p] * ϕ[p - j]
        ϕ[p - j] -= ϕ[p] * tmp
    end
        if isodd(p - 1) ϕ[div(p - 1, 2) + 1] *= (one(T) - ϕ[p]) end
    σ2[] *= (one(T) - ϕ[p]^2)
end

function levinson_durbin(z::AbstractVector, p::Integer)
    ρ = map(i -> autocorrelation(z, i), 1:p)
    T = typeof(zero(eltype(z)) / 1)
    ϕ = zeros(T, p)
    σ2 = Ref{T}(autocovariance(z, 0))
    for i in 1:p
        levinson_durbin!(ϕ, σ2, ρ, i)
    end
    return ϕ, σ2[]
end


text:
module MP3

# package code goes here
using SampledSignals
using FileIO
using FixedPointNumbers
using LinearAlgebra
using mpg123_jll
using LAME_jll

# methods to override
import SampledSignals: nchannels, nframes, samplerate, unsafe_read!, unsafe_write

# re-export
export Hz, kHz, s

struct MP3INFO
    nframes::Int64
    nchannels::Int32
    samplerate::Int64
    datatype::DataType
end

"""create an MP3INFO object from given audio buffer"""
function MP3INFO(buf::SampleBuf{T}) where {T}
    MP3INFO(nframes(buf), nchannels(buf), samplerate(buf), T)
end

include("readers.jl")
include("writers.jl")

function __init__()
    initialize_readers()
    initialize_writers()

    # MP3 files with ID3v1 (or no tags) and ID3v2 tags have different headers
    magic = (UInt8[0xff, 0xfb], UInt8[0x49, 0x44, 0x33])
    add_format(format"MP3", magic, [".mp3"], [:MP3])
end

end # module


text:
using QuantumPropagators
using QuantumControlBase
using Krotov
using QuantumControl
using Documenter

DocMeta.setdocmeta!(QuantumControl, :DocTestSetup, :(using QuantumControl); recursive=true)

makedocs(;
    modules=[QuantumPropagators, QuantumControlBase, QuantumControl, Krotov],
    checkdocs = :exports,
    authors="Michael Goerz <mail@michaelgoerz.net>, Alastair Marshall <alastair@nvision-imaging.com>, and contributors",
    repo="https://github.com/quantumcontrol-jl/QuantumControl.jl/blob/{commit}{path}#{line}",
    sitename="QuantumControl.jl",
    format=Documenter.HTML(;
        prettyurls=true,
        canonical="https://QuantumControl-jl.github.io/QuantumControl.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/QuantumControl-jl/QuantumControl.jl",
)


text:
type AdiaComputer <: QuComput
    ###########################################
    #  Eternal
    ###########################################
    HB::AbstractMatrix      # Base Hamiltonian
    HP::AbstractMatrix      # Problem Hamiltonian
    maxtime::Real           # max evolution time
    n::Int64           # number of bits
    dt::Real                # time step

    ###########################################
    # variable for current state
    ###########################################
    location::Real          # time location
    state::AbstractVector
    eigens::AbstractVecOrMat

    ###########################################
    # measure
    ###########################################

    prob::Real

    ###########################################
    # other
    ###########################################

    nev::Int

    # n is the number of bits
    # maxtime is the max evolution time
    function AdiaComputer{M,N}(ins::Instance{M,N},n::Int,maxtime::Real;dt=1e-2,nev=6)
        HB = bHamilton(n)
        HP = pHamilton(ins,n)
        # set time location to be the beginning
        location = 0
        # prepare the initial state
        state = convert(Array{Complex,1},[1/sqrt(2^n) for i=1:2^n])

        # initialize states
        eigens = eigs(HB,nev=nev,which=:SM)[1].'

        # adjust nev if the number of bits is smaller than 3
        if n<3
            nev = 2^n
            warn("adjusting nev to $(nev)\n")
        end

        # probility
        prob = norm([x==findmin(diag(HP))[2]?1:0 for x=1:2^n])^2

        new(HB,HP,maxtime,n,dt,location,state,eigens,prob,nev)
    end
end

function Hamiltonian(Hs::AdiaComputer)
    return (1-Hs.location)*Hs.HB+Hs.location*Hs.HP
end

export AdiaComputer

max_stars_repo_name: Roger-luo/QuCmp

max_issues_repo_name: Roger-luo/QuCmp

max_forks_repo_name: Roger-luo/QuCmp


text:
################################################################################
#
# Testcases for Polynomials.jl are added here.
#
################################################################################

using Polynomials
using FluxHELayers.Polynomials: Poly2, Poly3, Poly4, Poly5, Poly6, Poly7, Poly8, Poly9, Poly10

@testset "Poly2" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8])
    poly_test = Poly2(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.1
end

@testset "Poly3" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8, 7])
    poly_test = Poly3(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.1
end

@testset "Poly4" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8, 7, 6])
    poly_test = Poly4(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.1
end

@testset "Poly5" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8, 7, 6, 5])
    poly_test = Poly5(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.1
end

@testset "Poly6" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8, 7, 6, 5, 4])
    poly_test = Poly6(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.3
end

@testset "Poly7" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8, 7, 6, 5, 4, 3])
    poly_test = Poly7(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.1
end

@testset "Poly8" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8, 7, 6, 5, 4, 3, 2])
    poly_test = Poly8(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.1
end

@testset "Poly9" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8, 7, 6, 5, 4, 3, 2, 1])
    poly_test = Poly9(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.1
end

@testset "Poly10" begin
    initial_scale = 2.0^40

    poly = Polynomial([10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 1])
    poly_test = Poly10(tuple(poly.coeffs...))

    cases = [1.0, 1.25, 1.5, 1.75, 2.0]
    expected = map(x -> poly(x), cases)

    vector = encrypt_array(cases, initial_scale, SEALParams)
    output = poly_test(vector)
    output = decrypt_array(output, SEALParams)

    @test output[1:5] ≈ expected atol=0.1
end


text:
#Script that matches oligotype sequences to database of full length Synechococcus sequences

using BioSequences
using BioAlignments
using JLD2
using DataFrames, CSV

using Gadfly, Colors, Compose, ColorSchemes, Cairo, Fontconfig
include("alignment_functions.jl")
## load in that reference database!

#unclear exactly how this database was made in the first place....oops!
@load "/home/kristen/Documents/V6V8_analysis/analysis_products/marine_syn_updated_db.jld2" synref

#must change string sequences back to BioSequence types!
synref.sequence = DNASequence.(synref[!,:sequence])
synref.Length = length.(synref[!,:sequence])

## we need to do some clean up before can use it!

#exclude sequences that are less than full length: goes from 232 seqs to 229
ii=findall(synref[!,:Length] .> 700)
synref=synref[ii,:]

#to find out how many have unambiguous clades;
#findall((synref[!,:clade] .== "?") .| (synref[!,:clade] .== "-"))

## first shorten to V6-V8 region!
synref.V6V8=Array{Any,1}(missing,size(synref,1)) #keep the alignments for future reference!
for q=1:size(synref,1)
    synref[q,:V6V8]=find_V6V8_region(synref[q,:sequence])
end

synref.V6V8_length=map(x -> length(x),synref[!,:V6V8])

ii=findall(synref[!,:V6V8_length] .>= 300) #exclude shorter sequences - goes down to 225 sequences
synref=synref[ii,:]
#the double check:
#kk=findall(isempty.(synref[!,:V6V8]))

## alright, soooo, we need to show two types of aligment relationships:

# how doe strain reps in the database differ from each other?
# how are these connected to the oligotype sequences?

#matrix appears to be the easiest way to do this!

#we could just add the oligotype sequences in and designate as new 'clade'....

# ii=findall(synrefB[!,:clade] .== clade)
# x=synrefB[ii,:V6V8]
# idx = unique(z -> x[z], 1:length(x))

## per clade, find the unique sequences, including short ones:
clade_list=unique(synref[!,:clade]) #find unqiue sequences and reorganize!
# clade_list=sort(setdiff(clade_list,["-";"?"]))

#indexin(["5.3"],clade_list0)
## for this, it's just easier to do by hand!
clade_list=["5.3";
"5.2MV1";
 "5.2MV2";
 "5.2MV3";
 "CB4"   ;
 "CB5"   ;
 "I"    ;
 "II"    ;
 "III"   ;
 "IV"    ;
 "V"     ;
 "VI"    ;
 "VII"   ;
 "VIII";
  "IX" ;
 "CRD1" ;
 "WPC1"  ;
 "WPC2"  ;
 "XV"    ;
 "XVI"]
##
unqrefs=DataFrame(ID=String[],strain_reps=Set[],clade=String[],V6V8=BioSequence[])

for clade in clade_list

    println(string("analyzing sequences within clade: ",clade))
    jj=findall((synref[!,:clade] .== clade) .& (synref[!,:V6V8_length] .>300))
    withinclade_array=deepcopy(synref[jj,[:strain]])
    scoremodel = AffineGapScoreModel(match=5, mismatch=-4, gap_open=-10,gap_extend=-10)

    for q=1:size(jj,1)
        #println(q)
        temp_df=align_to_database2!(synref[jj[q],:V6V8],synref[jj,:],scoremodel,:V6V8,:strain) #align to unique sequences only here
        #synrefD[j,:alignments]=temp_df #record alignments
        tempID=synref[jj[q],:strain] #column name
        withinclade_array[Symbol(tempID)] = temp_df[!,:mismatch]
    end

    # now, sort through this identity matrix!
    global counter=1
    for q=1:size(jj,1)

        #first, check if this strain has already been accounted for:
        if !any(map(x -> withinclade_array[q,:strain] in x,unqrefs[!,:strain_reps])) #if not, then proceed!

            qq=findall(convert(Array{Any},withinclade_array[q,:]) .== 0)
            qq=qq .- 1 #to match columns

            if length(qq) == 1 #then just this strain has this sequence and it should be unique!
                push!(unqrefs,[string(clade,"-",counter),Set([synref[jj[q],:strain]]),clade,synref[jj[q],:V6V8]])
                global counter += 1
            else
                strains=Set(withinclade_array[qq,:strain])

                #indexes back into synref!
                qqind=indexin(withinclade_array[qq,:strain],synref[!,:strain])
                y=argmax(synref[qqind,:V6V8_length]) #index of longest v6v8 read in qqind

                push!(unqrefs,[string(clade,"-",counter),strains,clade,synref[qqind[y],:V6V8]])
                global counter += 1
            end
        end
    end
end

## sweet! and we can double check this:

findall(map(x -> "HB1133" in x,unqrefs[!,:strain_reps])) #this should just be one!
findall(map(x -> "RCC1026" in x,unqrefs[!,:strain_reps])) #this should just be one!

## add the lengths!
unqrefs.V6V8_length  = length.(unqrefs[!,:V6V8])

## and now, add in the oligotypes!

oligo_file = open(FASTA.Reader, "/home/kristen/Documents/V6V8_analysis/analysis_products/syn_oligotyping/mvco_syn/mvco_syn_timeseries_for_oligotyping_pynast_aligned-c19-s1-a0.0-A0-M200/oligo_repseqs.fasta") #test-c25-s1-a0.0-A100-M40/oligo_representative_no_gaps.fasta") #open the fasta file

for record in oligo_file
    println(string(FASTA.identifier(record)))
    push!(unqrefs, [string(FASTA.identifier(record)), Set([]),"",FASTA.sequence(record),length(FASTA.sequence(record))])
end
close(oligo_file)

## cool, now we can do alignments for each sequence to see how they match up!

clade_array=deepcopy(unqrefs[!,[:clade,:ID]])

scoremodel = AffineGapScoreModel(match=5, mismatch=-4, gap_open=-10,gap_extend=-10)
#scoremodel = AffineGapScoreModel(match=5, mismatch=-4, gap_open=-4,gap_extend=-1)

for q=1:size(unqrefs,1)

        println(q)
        temp_df=align_to_database2!(unqrefs[q,:V6V8],unqrefs,scoremodel,:V6V8,:ID) #align to unique sequences only here
        #synrefD[j,:alignments]=temp_df #record alignments
        tempID=unqrefs[q,:ID]
        clade_array[Symbol(tempID)] = temp_df[!,:mismatch]

end

## a double check on the oligotypes:

k=findall(clade_array[!,:ID] .== "O3")
kk=findall(convert(Array{Any},clade_array[k,:]) .== 0)
qq= [ i[2] for i in kk ]
qq=qq .- 2
println(clade_array[qq,:ID])
## sweet! and the plotting!

spy(convert(Array,clade_array[:,3:end]))

xlabel=unqrefs[!,:ID]

xticks = collect(1:size(unqrefs,1))
labels2 = Dict(zip(xticks, xlabel))

spy(convert(Array,clade_array[:,3:end]),
    Scale.x_discrete(labels = x -> labels2[x]),
    Scale.y_discrete(labels = x -> labels2[x]))

## save these matches

@save "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_clade_match_dataframes.jld2" unqrefs clade_array

## and export for fig in matlab, if needed:

CSV.write("/home/kristen/Documents/V6V8_analysis/scripts/rel_abnd_analysis_scripts/matlab_scripts/V6V8_matrix_match.csv", clade_array, delim = ',')

# notes; ah-ya! O3 has additional matches to M11.1 and M11.2, but clade designation is ambiguous from full length 16S
#O14 also matches M12.1, which is clade II from 16S...

## SI table to show sequences used:

synref_for_table=deepcopy(synref[!,[:clade,:strain,:accession,:Source,:V6V8_length]])

## some housekeeping:

#exclude the strains that don't have clade designation:
ii=findall((synref[!,:clade] .== "?") .| (synref[!,:clade] .== "-"))
jj=setdiff(collect(1:size(synref,1)),ii)
synref_for_table=synref_for_table[jj,:]

#replace!(synref_for_table[!,:Source],"Hunter-Cevera et al. 2016" => "Hunter-Cevera 2014")

## now match each strain to a unqiue V6V8 sequence yo!

v6v8_match=Array{Union{Missing, String}}(missing, size(synref_for_table,1))
for q=1:size(synref_for_table,1)
    println(q)
    qq=findall(map(x -> synref_for_table[q,:strain] in x,unqrefs[!,:strain_reps]))
    v6v8_match[q] = unqrefs[qq,:ID][1]
end

synref_for_table.V6V8_unq_match = v6v8_match

##

sort!(synref_for_table,:V6V8_unq_match)
CSV.write("/home/kristen/Documents/V6V8_analysis/paper_V6V8/PNAS_style/syn_ref_table.csv", synref_for_table, delim = ",")

## If need to write out these sequences in a fasta format:

synref_for_table_seqs=deepcopy(synref[!,[:clade,:strain,:accession,:Source,:sequence]])

#exclude the strains that don't have clade designation:
ii=findall((synref[!,:clade] .== "?") .| (synref[!,:clade] .== "-"))
jj=setdiff(collect(1:size(synref,1)),ii)
synref_for_table_seqs=synref_for_table_seqs[jj,:]

fastafile="/home/kristen/Documents/V6V8_analysis/analysis_products/marine_syn_db.fasta"

w = open(FASTA.Writer, fastafile)
for q=1:size(synref_for_table_seqs,1)
    rec = FASTA.Record(string(synref_for_table_seqs[q,:strain]),string(synref_for_table_seqs[q,:accession]),synref_for_table_seqs[q,:sequence])
    write(w, rec)
end
close(w)

## or if want to put this table in latex:
t=Array{Union{Missing, String}}(missing, size(synref_for_table,1))
t .= "&"

n=Array{Union{Missing, String}}(missing, size(synref_for_table,1))
n .= "\\"

insert!(synref_for_table, 2, t, :L1)
insert!(synref_for_table, 4, t, :L2)
insert!(synref_for_table, 6, t, :L3)
insert!(synref_for_table, 8, t, :L4)
insert!(synref_for_table, 10, t, :L5)
insert!(synref_for_table, 12, n, :E)

##

CSV.write("/home/kristen/Documents/V6V8_analysis/scripts/rel_abnd_analysis_scripts/matlab_scripts/syn_ref_table.csv", synref_for_table, delim = "   ")

## if need to double check alignments directly:

scoremodel = AffineGapScoreModel(match=5, mismatch=-4, gap_open=-10,gap_extend=-10)

# s1 = synrefB[ii[idx[1]],:V6V8];
# s2 = synrefB[ii[idx[2]],:V6V8];
s1 = unqrefs[12,:V6V8]#unqseq[6,:V6V8];
s2 = unqrefs[18,:V6V8];
res = pairalign(LocalAlignment(), s1, s2, scoremodel)
aln=alignment(res)

max_stars_repo_name: hsosik/NES-LTER

max_issues_repo_name: hsosik/NES-LTER

max_forks_repo_name: hsosik/NES-LTER


text:
using PolyChaos, Test
import LinearAlgebra: norm
myfile = open("dataQuadratureRules/config.txt")
ns = parse.(Int,readlines(myfile))
close(myfile)
myfile = open("dataQuadratureRules/endPts.txt")
endPts = parse.(Float64,readlines(myfile))
close(myfile)
nodes = ns[1]:ns[2]:ns[3]
high = ns[3]
tol = 1e-7

@time @testset "Fejer" begin
    for n in nodes
            myfile = open("dataQuadratureRules/fejer$n.txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = fejer(n)
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            close(myfile)
    end
end

@time @testset "Fejer2" begin
    for n in nodes
            myfile = open("dataQuadratureRules/fejer2_$n.txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = fejer2(n)
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            close(myfile)
    end
end

@time @testset "Clenshaw Curtis" begin
    for n in nodes
            myfile = open("dataQuadratureRules/cc$n.txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = clenshaw_curtis(n)
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            close(myfile)
    end
end

@time @testset "Gauss" begin
        data1 = open("dataQuadratureRules/logHigh.txt")
        d1 = reshape(parse.(Float64,readlines(data1)),:,2)
        data2 = open("dataQuadratureRules/hermiteHigh.txt")
        d2 = reshape(parse.(Float64,readlines(data2)),:,2)
        close(data1)
        close(data2)
    for n in ns
            myfile = open("dataQuadratureRules/gaussLog$n.txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = gauss(n,d1[:,1],d1[:,2])
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            q = Quad(n,d1[:,1],d1[:,2],Measure("gaussian"))
            @test isapprox(norm(αβref - vcat(nw(q)...),Inf),0.;atol=tol)
            close(myfile)
            myfile = open("dataQuadratureRules/gaussHerm$n.txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = gauss(n,d2[:,1],d2[:,2])
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            q = Quad(n,d2[:,1],d2[:,2],Measure("hermite"))
            @test isapprox(norm(αβref - vcat(nw(q)...),Inf),0.;atol=tol)
            close(myfile)
    end
end

@time @testset "Radau with Logistic" begin
        data1 = open("dataQuadratureRules/logHigh.txt")
        d1 = reshape(parse.(Float64,readlines(data1)),:,2)
        close(data1)
    for n in ns
        for i in [endPts;-endPts]
            myfile = open("dataQuadratureRules/radauLog$(n)pt$i.txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = radau(n,d1[:,1],d1[:,2],i)
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            close(myfile)

        end
    end
end

@time @testset "Radau with Hermite" begin
        data2 = open("dataQuadratureRules/hermiteHigh.txt")
        d2 = reshape(parse.(Float64,readlines(data2)),:,2)
        close(data2)
    for n in ns
        for i in [endPts;-endPts]
            myfile = open("dataQuadratureRules/radauHerm$(n)pt$i.txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = radau(n,d2[:,1],d2[:,2],i)
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            close(myfile)
        end
    end
end

function testInt(lb,ub)
        # analytic integral of f(x) = exp(-|x|)
        @assert lb < ub "lower bound has to be smaller than upper bound"
        lb > 0 && return -exp(-ub) + exp(-lb)
        ub < 0 && return exp(ub) - exp(lb)
        return 2 - exp(lb) - exp(-ub)
end

ns[3]=700
intervals = [(-5.,5.),(-Inf,Inf),(-Inf,-1.),(1.,Inf),(-3.,5.),(-10.,1.)]
@time @testset "gpquad" begin
        for n in intervals
                nod1,wei1 = quadgp(x->exp.(-abs.(x)),n[1],n[2],ns[3],quadrature=fejer)
                #print("Nodes: ",nod1,"\n")
                #print("Weights: ",wei1,"\n")
                analytic = testInt(n[1],n[2])
            @test isapprox(sum(wei1)-analytic,0;atol=1e-4)
        end
        for n in intervals
                nod2,wei2 = quadgp(x->1.,n[1],n[2],ns[3],quadrature=fejer)
                analytic = testInt(n[1],n[2])
            @test isapprox(exp.(-abs.(nod2))'*wei2-analytic,0;atol=1e-4)
        end
end


text:
#Shows what a typical ISD and Viral Infections looks like
using RCall
include("ProblemGenerator.jl")
include("VirusCallBacks.jl")
########################################################################
#Heatmaps of IFN dynamics (Figure 2)
########################################################################

#-----------ISD, Not Stochastic, Homo-----------
probISD_nS_Ho = ModelSetup(:ISD,:notStochastic,:Homo)
    solISD_nS_Ho = @time solve(probISD_nS_Ho,CVODE_BDF(linear_solver=:GMRES),saveat=0.1)
ISD_nS_Ho = DataFrame()
ISD_nS_Ho.x = repeat(1:N,N)
ISD_nS_Ho.y = repeat(1:N,inner=N)
ISD_nS_Ho.IFN = vec(solISD_nS_Ho(10.0)[:,:,7])

#-----------Virus, Not Stochastic, Homo-----------
probVirus_nS_Ho = ModelSetup(:Virus,:notStochastic,:Homo)
    solVirus_nS_Ho = @time solve(probVirus_nS_Ho,CVODE_BDF(linear_solver=:GMRES),saveat=0.1,callback=cb)
Virus_nS_Ho = DataFrame()
Virus_nS_Ho.x = repeat(1:N,N)
Virus_nS_Ho.y = repeat(1:N,inner=N)
Virus_nS_Ho.IFN = vec(solVirus_nS_Ho(10.0)[:,:,7])

#-----------Virus, Not Stochastic, Homo, Cell States-----------
function CellStateGrid(θ::ParContainer,t::Float64)
  grid = zeros(Int64,size(θ.cellsInfected))
  grid[θ.cellsInfected .<= t] .= 1
  grid[θ.cellsDead .<= t] .= 2

  return grid
end

CellStateVirus_nS_Ho_raw = CellStateGrid(probVirus_nS_Ho.p,10.0)
CellStateVirus_nS_Ho = DataFrame()
CellStateVirus_nS_Ho.x = repeat(1:N,N)
CellStateVirus_nS_Ho.y = repeat(1:N,inner=N)
CellStateVirus_nS_Ho.Cell = vec(CellStateVirus_nS_Ho_raw)

#-----------ISD, Stochastic, Hetero-----------
probISD_S_He = ModelSetup(:ISD,:Stochastic,:Hetero)
    solISD_S_He = @time solve(probISD_S_He,CVODE_BDF(linear_solver=:GMRES),saveat=0.1)
ISD_S_He = DataFrame()
ISD_S_He.x = repeat(1:N,N)
ISD_S_He.y = repeat(1:N,inner=N)
ISD_S_He.IFN = vec(solISD_S_He(10.0)[:,:,7])

#-----------Virus, Stochastic, Hetero-----------
probVirus_S_He = ModelSetup(:Virus,:Stochastic,:Hetero)
    solVirus_S_He = @time solve(probVirus_S_He,CVODE_BDF(linear_solver=:GMRES),saveat=0.1,callback=cb)
Virus_S_He = DataFrame()
Virus_S_He.x = repeat(1:N,N)
Virus_S_He.y = repeat(1:N,inner=N)
Virus_S_He.IFN = vec(solVirus_S_He(10.0)[:,:,7])

#-----------Virus, Not Stochastic, Homo, Cell States-----------

CellStateVirus_S_He_raw = CellStateGrid(probVirus_S_He.p,10.0)
CellStateVirus_S_He = DataFrame()
CellStateVirus_S_He.x = repeat(1:N,N)
CellStateVirus_S_He.y = repeat(1:N,inner=N)
CellStateVirus_S_He.Cell = vec(CellStateVirus_S_He_raw)

#-----------R Plot-----------
@rput ISD_nS_Ho
@rput Virus_nS_Ho
@rput CellStateVirus_nS_Ho
@rput ISD_S_He
@rput Virus_S_He
@rput CellStateVirus_S_He
R"""
  library(ggplot2)
  library(ggpubr)

  #To draw an arc on the heatmap to show initial condition
  arc <- data.frame(
  x0 = seq(0, 200, length.out = 1000),
  y0 = sqrt(200^2 - (seq(0, 200, length.out = 1000)^2))
  )

  #Color Range
  clim = c(min(ISD_S_He$IFN), max(ISD_nS_Ho$IFN))
  commonFigureOptions <- list(scale_fill_distiller(palette = "Spectral",guide = guide_colorbar(frame.colour = "black", ticks.colour = "black"),limits = clim),
    scale_x_continuous(expand=c(0,0)),
    scale_y_continuous(expand=c(0,0)),
    theme_bw(base_size = 12),
    ylab("Cell"),
    xlab("Cell"),
    labs(fill="IFN (nM)"),
    theme(plot.title = element_text(hjust = 0.5),aspect.ratio = 1))

  p1 <- ggplot() +
    geom_raster(ISD_nS_Ho, mapping =aes(x, y, fill=IFN)) +
    geom_line(data=arc, aes(x=x0, y=y0),size=0.2,linetype="dashed") +
    labs(fill="IFN (nM)") +
    ggtitle("ISD \n Deterministic \n Homogeneous") +
    commonFigureOptions

  p2 <- ggplot(Virus_nS_Ho, aes(x, y, fill=IFN)) +
    geom_raster(aes(fill=IFN)) +
    ggtitle("Virus \n Deterministic \n Homogeneous") +
    commonFigureOptions

  p3 <- ggplot(CellStateVirus_nS_Ho, aes(x, y, fill=factor(Cell))) +
    geom_raster(aes(fill=factor(Cell))) +
    scale_fill_manual(values=c("#377EB8","#4DAF4A","#E41A1C"),labels = c("Healthy", "Infected","Dead")) +
    scale_x_continuous(expand=c(0,0)) +
    scale_y_continuous(expand=c(0,0)) +
    theme_bw() +
    labs(title ="Viral \n Infection \n Cell States", x="Cell", y="Cell", fill = "Cell State") +
    theme(plot.title = element_text(hjust = 0.5),aspect.ratio = 1)

  p4 <- ggplot() +
    geom_raster(ISD_S_He, mapping =aes(x, y, fill=IFN)) +
    geom_line(data=arc, aes(x=x0, y=y0),size=0.2,linetype="dashed") +
    ggtitle("ISD \n Stochastic \n Heterogeneous") +
    commonFigureOptions

  p5 <- ggplot(Virus_S_He, aes(x, y, fill=IFN)) +
    geom_raster(aes(fill=IFN)) +
    ggtitle("Virus \n Stochastic \n Heterogeneous") +
    commonFigureOptions

  p6 <- ggplot(CellStateVirus_S_He, aes(x, y, fill=factor(Cell))) +
    geom_raster(aes(fill=factor(Cell))) +
    scale_fill_manual(values=c("#377EB8","#4DAF4A","#E41A1C"),labels = c("Healthy", "Infected","Dead")) +
    scale_x_continuous(expand=c(0,0)) +
    scale_y_continuous(expand=c(0,0)) +
    theme_bw() +
    labs(title ="Viral \n Infection \n Cell States",x="Cell", y="Cell", fill = "Cell State") +
    theme(plot.title = element_text(hjust = 0.5),aspect.ratio = 1)

  figure1 <- ggarrange(p1, p2, p3, p4, p5, p6,
                      labels = "AUTO",
                      align = "h",
                      ncol = 3, nrow = 2)

  ggsave("./Figures/Figure2New.pdf",width=12,height=8,units="in")
"""


text:
using SafeTestsets
@safetestset "aI Test aITaskDispatcher" begin include("aITestTaskDispatcher.jl") end


text:
using Interpolations
using HDF5
using StaticArrays

try
    global const DIRECTORY_SL = dirname(dirname(@__FILE__)) * "/Data/" # Directory in which the sl are stored
    namefile = DIRECTORY_SL * "data_sl_l_" * string(2) * ".hf5"
    file = h5open(namefile,"r")
    close(file)
catch
    println("File Coupling_coefficients_interpolated - ERROR - could not load the sl coefficients from the file")
end

const ALPHA_MIN = 10^(-3) # Minimum value of α over which we interpolated the sₗ
const ALPHA_MAX = 10^(0) # Maximum value of α over which we interpolated the sₗ
const ECC_MIN = 0.0 # Minimum value of e over which we interpolated the sₗ
const ECC_MAX = 0.99 # Maximum value of e over which we interpolated the sₗ
const LMAX = 50 # largest mode we can compute from the interpolated sₗ
const LMAX_HALF = div(50, 2)

function Interpolating_function_sl(l::Int64)
    #### Returns one interpolating function sₗ ####
    #Getting the name of the file
    namefile = DIRECTORY_SL*"data_sl_l_"*string(l)*".hf5" # Name of the file

    # Reading the file
    file = h5open(namefile,"r") # Opening the file in read-only mode
    tab_logalpha = read(file,"tab_logalpha") # Range of log(alpha)
    tab_ein  = read(file,"tab_ein")  # Range of ein
    tab_eout = read(file,"tab_eout") # Range of eout
    tab_sl = read(file,"tab_sl") # Table of sl
    close(file) # Closing the file

    # getting the bounds and steps
    logalpha_min, delta_logalpha, logalpha_max = tab_logalpha # Values for the range in log(alpha)
    ein_min,  delta_ein,  ein_max  = tab_ein  # Values for the range in ein
    eout_min, delta_eout, eout_max = tab_eout # Values for the range in eout

    # getting the lists of points for which we have exact values
    range_logalpha = logalpha_min:delta_logalpha:logalpha_max # Range in log(alpha)
    range_ein = ein_min:delta_ein:ein_max # Range in ein
    range_eout = eout_min:delta_eout:eout_max # Range in eout

    # getting the interpolation function -- 229 allocations / 100 MiB
    INT_sl = scale(interpolate(tab_sl, BSpline(Cubic(Line(OnGrid())))),
                   range_logalpha,range_ein,range_eout) # Making the 3D gridded interpolation

    # Returning the interpolation function
    return INT_sl
end

const TAB_INT_SL = SVector{LMAX_HALF}([Interpolating_function_sl(l) for l = 2:2:LMAX]) # This is a list of the interpolating function. It is constant that can be accessed by all scripts

function Coupling_coefficients_interpolated(K::Array{Tuple{Float64,Float64,Float64},1}, lmax::Int64; G::Float64 = 1.0)
    #### returns an array of size n_pop * n_pop * lmax_half. Each cell contains an (even axisymmetric) coupling coefficient ####
    # getting all the constants we will need
    n_pop = length(K)
    lmax_half = div(lmax, 2)

    # preallocating all the arrays we will need
    J = zeros(n_pop, n_pop, lmax_half) # will store all the coupling coefficients

    # Let's filling up J
    for star_1 = 1:n_pop
        m, a, ecc = K[star_1] # (m,a,e) of the unprimed particle

        for star_2 = 1:star_1
            mp, ap, eccp = K[star_2]  # (m,a,e) of the primed particle

            if (a <= ap) # Unprimed particle is inner, and primed is outer
                ain, ein = a, ecc # Inner particle
                aout, eout = ap, eccp # Outer particle
            else # Primed particle is inner, and unprimed is outer
                ain, ein = ap, eccp # Inner particle
                aout, eout = a, ecc # Outer particle
            end

            alpha = ain / aout # Ratio of the semi-major axes, with alpha <= 1.0

            # Checking that (alpha,ein,eout) are within the interpolated ranges
            @assert ALPHA_MIN <= alpha <= ALPHA_MAX "get_Jl: alpha not in the interpolated range"
            @assert ECC_MIN   <= ein   <= ECC_MAX   "get_Jl: ein not in the interpolated range"
            @assert ECC_MIN   <= eout  <= ECC_MAX   "get_Jl: eout not in the interpolated range"

            for l_half = 1:lmax_half
                l = 2* l_half

                # Checking that l is within the interpolated range
                @assert 2 <= l <= LMAX "get_Jl: l not in the interpolation range"

                # Reading sl
                logalpha = log(alpha) # Value of log(alpha), i.e. the variable w.r.t. which the interpolation was performed
                l_half = div(l, 2) # index of mode l
                sl = TAB_INT_SL[l_half](logalpha, ein, eout) # Getting the value from the interpolated functions

                # Computing the coupling coefficient
                Jl = (sl * G * m * mp) / aout # Jₗ[K₁, K₂] = G m₁ m₂ sₗ[K₁, K₂] / a_out by definition

                # Storing the coupling coefficient -- We use the symmetry
                J[star_1, star_2, l_half] = Jl
                J[star_2, star_1, l_half] = Jl
            end
        end
    end

    return(J)
end

function Coupling_coefficient_interpolated_v2(l::Int64,
                                                K::Tuple{Float64,Float64,Float64}, Kp::Tuple{Float64,Float64,Float64};
                                                G::Float64 = 1.0)
    #### Returns a single (even axisymmetric) coupling coefficient ####

    # Giving a name to the annuli parameters
    m,  a,  ecc  = K  # (m,a,e) of the unprimed particle
    mp, ap, eccp = Kp # (m,a,e) of the primed particle

    if (a <= ap) # Unprimed particle is inner, and primed is outer
        ain, ein = a, ecc # Inner particle
        aout, eout = ap, eccp # Outer particle
    else # Primed particle is inner, and unprimed is outer
        ain, ein = ap, eccp # Inner particle
        aout, eout = a, ecc # Outer particle
    end

    alpha = ain / aout # Ratio of the semi-major axes, with alpha <= 1.0

    ##### Checking that (alpha,ein,eout) are within the interpolated ranges
    @assert ALPHA_MIN <= alpha <= ALPHA_MAX "get_Jl: alpha not in the interpolated range"
    @assert ECC_MIN   <= ein   <= ECC_MAX   "get_Jl: ein not in the interpolated range"
    @assert ECC_MIN   <= eout  <= ECC_MAX   "get_Jl: eout not in the interpolated range"
    @assert 2 <= l <= LMAX "get_Jl: l not in the interpolation range"

    # Reading sl
    logalpha = log(alpha) # Value of log(alpha), i.e. the variable w.r.t. which the interpolation was performed
    l_half = div(l, 2) # index of mode l
    sl = TAB_INT_SL[l_half](logalpha, ein, eout) # Getting the value from the interpolated functions

    # Computing the coupling coefficient
    Jl = (sl * G * m * mp) / aout # Jₗ[K₁, K₂] = G m₁ m₂ sₗ[K₁, K₂] / a_out by definition

    # Returning the coupling coefficient
    return Jl # Output
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe


text:
#Pkg.add("ExcelFiles")
# Info: Using the Python distribution in the Conda package by default.
# To use a different Python version, set ENV["PYTHON"]="pythoncommand" and re-run Pkg.build("PyCall").

using ExcelFiles, DataFrames

years = [Symbol("$i") for i in 2006:2015];

demand = DataFrame(load("Data.xlsx", "Load!A1:K8761"));

head(demand)

demand = demand[2:end];

demand.colindex

size(demand)

windSouthEast = DataFrame(load("Data.xlsx", "SouthEast!A1:K8761"));

head(windSouthEast)

windNorth = DataFrame(load("Data.xlsx", "North!A1:K8761"));

windNorth =windNorth[2:end];

windSouthWest = DataFrame(load("Data.xlsx", "SouthWest!A1:K8761"));

windSouthWest = windSouthWest[2:end];

#SouthsolarInput rng=PINCHER-CREEK
solarSouth = DataFrame(load("Data.xlsx", "PINCHER-CREEK"));

solarSouth = solarSouth[2:end-1];

# par= MiddleSolarInput rng=EDMONTON-STONY-PLAIN-CS!A1:K8761 Rdim=1 Cdim=1
# par= NorthSolarInput rng=GRANDE-PRAIRIE-A!A1:K8761 Rdim=1 Cdim=1
solarMiddle = DataFrame(load("Data.xlsx", "EDMONTON-STONY-PLAIN-CS!A1:K8761"));

#solarMiddle

solarMiddle = solarMiddle[2:end];

solarNorth = DataFrame(load("Data.xlsx", "GRANDE-PRAIRIE-A!A1:K8761"));

solarNorth = solarNorth[2:end];

#solarNorth[Symbol("2006")]

size(solarSouth)

#solarSouth

#solarSouth

typeof(solarSouth)

#solarSouth1 = similar(solarSouth)

#* the unit of output in data table is w. transfor to MW. A solar farm with 50000 PV panel? may cause floating point problem, number is too small

for (i,j) in eachcol(solarSouth) # a tuple with symbol and array
  solarSouth[i] = j*50000/1000000;
end

for (i,j) in eachcol(solarMiddle) # a tuple with symbol and array
  solarMiddle[i] = j*50000/1000000;
end

for (i,j) in eachcol(solarNorth) # a tuple with symbol and array
  solarNorth[i] = j*50000/1000000;
end

#Para.csv
using CSV
para =  CSV.read("Para.csv")

para

para[1].== "wind"

turbine = parse(para[ para[1].== "wind",:UnitCap][1]);

turbine

panel =  parse(para[ para[1].== "solar",:UnitCap][1]);
panel

generation_assets = [Symbol(i) for i in para[:,1][2:end]]

fossil_fuel = [Symbol(i) for i in para[:,1][end-1:end]]

renewables = [Symbol(i) for i in para[:,1][2:3]]

windLoc = [:se, :sw, :n]
solarLoc = [:s, :m, :n]

size(para)

#para.columns[1]

para.colindex

FixOMCosts = Dict( zip( generation_assets    , ([parse(i) for i in para[:FixOMCost][2:end]]) )   )

VarOMCosts = Dict( zip( generation_assets    , ([parse(i) for i in para[:VarOMCost][2:end]]) )   )

ONights= Dict( zip( generation_assets    , ([parse(i) for i in para[:ONight][2:end]]) )   )

UnitCaps  = Dict( zip( generation_assets    , ([parse(i) for i in para[:UnitCap][2:end]]) )   )

Maxis  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Max][2:end]]) )   )

Emits  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Emit][2:end]]) )   )

Ramps  = Dict( zip( generation_assets    , ([parse(i) for i in para[:Ramp][2:end]]) )   )

LTimes= Dict( zip( generation_assets    , ([parse(i) for i in para[:LTime][2:end]]) )   )

# # parameters
# ''''''
# round = Para("battery", "Round") ;
# duration = Para("battery", "duration") ;
# taxEm = Para("CT","CTax");
# turbine = Para("wind","UnitCap");
# maxturbine = Para("wind","Max");
# panel = Para("solar","UnitCap");
# maxPanel = Para("solar","Max");
# rate = Para("CT","Discount");
# ''''''

maximums = Dict( zip( generation_assets    , ([parse(i) for i in para[:Max][2:end]]) )   )

roundEfficiency = Dict( zip( generation_assets    , ([parse(i) for i in para[:Round][2:end]]) )   )

duration= Dict( zip( generation_assets    , ([i for i in para[:duration][2:end]]) )   )

taxEm = Dict( zip( generation_assets    , ([parse(i) for i in para[:Ctax][2:end]]) )   )

rate = Dict( zip( generation_assets    , ([parse(i) for i in para[:Discount][2:end]]) )   )

Vcost = Dict( zip( generation_assets    , ( VarOMCosts[i] + taxEm[i]* Emits[i]  for i in generation_assets) )   ) ;

Vcost

Fcost =   Dict( zip(generation_assets, (ONights[i]*rate[i]*((1+rate[i])^LTimes[i])/(((1+rate[i])^LTimes[i])-1) +
                                            FixOMCosts[i] for i in generation_assets)
                    )
                );

Fcost

using JuMP
using Gurobi

m = Model(solver=GurobiSolver())

T = 8760

@variable(m, batteryCap>=0)

# energy in battery
@variable(m, v[1:T]>=0)

@variable(m, charge[1:T]>=0)
@variable(m, discharge[1:T]>=0)

@variable(m, fuelCap[fossil_fuel] >=0)

@variable(m, gen[ 1:T,fossil_fuel] >=0)

@variable(m, nSolarPlant[solarLoc] >=0)
@variable(m, nWindTurbine[windLoc] >=0)

wind = DataFrame(s = 1.0:T,m = 1.0:T, n = 1.0:T);

solar = DataFrame(s = 1.0:T,m = 1.0:T, n = 1.0:T);

year = 2006

Symbol("$year")

wind[:se] = windSouthEast[Symbol("$year")];
wind[:sw] = windSouthWest[Symbol("$year")];
wind[:n] = windNorth[ Symbol("$year")];

#* solar unit output per hour  (MW), number of panel is in 50000, so total solar ouput is MW
solar[:s] = solarSouth[Symbol("$year")];
solar[:m] = solarMiddle[Symbol("$year")];
solar[:n] = solarNorth[Symbol("$year")];

demand0 = demand[Symbol("$year")];

@objective(m, Min, sum(sum(wind[t,L]*nWindTurbine[L]*Vcost[:wind] for L in windLoc) +
                            (sum(solar[t,Ls]*nSolarPlant[Ls]*Vcost[:solar] for Ls in solarLoc))  +
                            (sum(gen[t,f] *Vcost[f] for f in fossil_fuel )) +
                            discharge[t]*Vcost[:battery]
                            for t in 1:T
                            ) +
                    sum(Fcost[:wind]*nWindTurbine[L]*turbine for L in windLoc) +
                    sum(Fcost[:solar]*nSolarPlant[Ls]*turbine for Ls in solarLoc) +
                    sum(Fcost[f]*fuelCap[f] for f in fossil_fuel) +
                    Fcost[:battery]*batteryCap
            );

@constraint(m,  sum(wind[1,L]*nWindTurbine[L] for L in windLoc) +
                                sum(solar[1,Ls]*nSolarPlant[Ls] for Ls in solarLoc) +
                                  sum(gen[1,f] for f in fossil_fuel) +
                                   discharge[1] -charge[1] >= demand0[1]);

# for t in 1:T
#     @constraint(m, Balance[t], sum(wind[t,L].*nWindTurbine[L] for L in windLoc) +
#                                 sum(solar[t,Ls].*nSolarPlant[Ls] for Ls in solarLoc) +
#                                   sum(gen[t,f] for f in fossil_fuel) +
#                                    discharge[t] -charge[t] >= demand0[t]);
# end

@constraint(m, Balance0[t=1:T], sum(wind[t,L].*nWindTurbine[L] for L in windLoc) +
                            sum(solar[t,Ls].*nSolarPlant[Ls] for Ls in solarLoc) +
                              sum(gen[t,f] for f in fossil_fuel) +
                               discharge[t] -charge[t] >= demand0[t]);

 @constraint(m, capacityFuel[t=1:T, f = fossil_fuel], gen[t,f] <= fuelCap[f] );

 @constraint(m, capacityWind[L=windLoc], nWindTurbine[L] <= maximums[:wind] );

 @constraint(m, capacitySolar[Ls=solarLoc], nSolarPlant[Ls] <= maximums[:solar] );

@constraint(m, RampUp[t=2:T, f = fossil_fuel], gen[t, f] <= gen[t-1,f]+Ramps[f]*fuelCap[f] );

@constraint(m, RampDn[t=2:T, f = fossil_fuel], gen[t, f] >= gen[t-1,f]-Ramps[f]*fuelCap[f] );

@constraint(m, Storage[t=2:T], v[t] ==   v[t-1]+(roundEfficiency[:battery]*charge[t-1]-discharge[t-1]));

@constraint(m, MaxStorage[t=1:T], v[t] <= duration[:battery] * batteryCap);

@constraint(m, MaxCharge[t=1:T], charge[t] <= duration[:battery] * batteryCap);

@constraint(m, MaxDischarge[t=1:T], discharge[t] <= batteryCap);

solve(m)

# # include("examples_julia/2_uc.jl");

# n_machines = 3;
# time_steps = 12;

# fixed_costs = [5, 10, 1000];
# variable_costs = [800, 700, 20];
# max_power = [500, 700, 1500];
# startup_costs = [10, 5, 850];
# demands = [1500, 1200, 1800, 2000, 2200, 1750, 1600, 1500, 1200, 1500, 1600, 1800];

# using JuMP
# using Gurobi

# m = Model(solver=GurobiSolver())

# @variable(m, on[1:time_steps, 1:n_machines], Bin)
# @variable(m, power[1:time_steps, 1:n_machines] >= 0)
# @variable(m, start[1:time_steps, 1:n_machines], Bin)

# @objective(m, Min, sum(dot(fixed_costs, vec(on[t, :])) +
#         dot(variable_costs, vec(power[t, :])) + dot(startup_costs, vec(start[t, :])) for t in 1:time_steps))

# for t in 1:time_steps
#     for i in 1:n_machines
#         @constraint(m, power[t, i] <= max_power[i] * on[t, i])
#     end
#     @constraint(m, sum(power[t, :]) == demands[t])
# end

# for i in 1:n_machines
#     @constraint(m, start[1, i] == on[1, i])
# end

# for t in 2:time_steps
#     for i in 1:n_machines
#         @constraint(m, start[t, i] >= on[t, i] - on[t - 1, i])
#     end
# end

# for t in 1:time_steps
#     for i in 1:n_machines
#         @constraint(m, on[t, i] >= start[t, i])
#         if t <= time_steps - 1
#             @constraint(m, on[t + 1, i] >= start[t, i])
#         end
#         if t <= time_steps - 2
#             @constraint(m, on[t + 2, i] >= start[t, i])
#         end
#     end
# end

# print(m)
# solve(m)
# getvalue(on)

Gurobi.Env()


text:
abstract type ActivationFunction end
abstract type NaiveActivationFunction <: ActivationFunction end
abstract type TrainedActivationFunction <: ActivationFunction end

struct Cube <: NaiveActivationFunction end
(::Cube)(x::T) where {T <: Number} = x^3.0

struct Linear <: NaiveActivationFunction end
(::Linear)(x::T) where {T <: Number} = x

struct ReLU <: NaiveActivationFunction end
(::ReLU)(x::T) where {T <: Number} = max(zero(x), x)

struct Sigmoid <: NaiveActivationFunction end
(::Sigmoid)(x::T) where {T <: Number} = 1.0 / (1.0 + exp(-x))

struct Square <: NaiveActivationFunction end
(::Square)(x::T) where {T <: Number} = x^2.0

struct Tanh <: NaiveActivationFunction end
(::Tanh)(x::T) where {T <: Number} = tanh(x)

struct KDEComparator{T₁ <: KernelDensity.UnivariateKDE,
                     T₂ <: KernelDensity.UnivariateKDE} <: TrainedActivationFunction
    pdf₁::KernelDensity.InterpKDE{T₁}
    pdf₂::KernelDensity.InterpKDE{T₂}

    function KDEComparator(samples::T₁,
                           targets::T₂,
                           weights::T₃;
                           boundary_offset::Int = 1) where {T₁ <: AbstractMatrix,
                                                            T₂ <: AbstractVector,
                                                            T₃ <: AbstractVector}
        X = samples
        t = targets
        w = weights

        h = vec(w' * X)
        h₀ = h[t .!= 1]
        h₁ = h[t .== 1]

        boundary = (min(minimum(h₀), minimum(h₁)) - boundary_offset,
                    max(maximum(h₀), maximum(h₁)) + boundary_offset)

        kde₀ = KernelDensity.kde_lscv(h₀, npoints = 100, boundary = boundary)
        kde₁ = KernelDensity.kde_lscv(h₁, npoints = 100, boundary = boundary)

        pdf₀ = KernelDensity.InterpKDE(kde₀)
        pdf₁ = KernelDensity.InterpKDE(kde₁)

        param(::KernelDensity.InterpKDE{T}) where {T} = T
        new{param(pdf₀),param(pdf₁)}(pdf₀, pdf₁)
    end
end

function (c::KDEComparator)(x::T) where {T <: Number}
    Distributions.pdf(c.pdf₁, x) - Distributions.pdf(c.pdf₂, x)
end

function train_kde_comparators(::Type{T₁},
                               samples::T₂,
                               targets::T₃,
                               n_neurons::Int;
                               boundary_offset::Int = 1) where {T₁ <: Number,
                                                                T₂ <: AbstractMatrix,
                                                                T₃ <: AbstractVector}
    X = samples
    t = targets
    L = n_neurons
    D = first(size(X))
    W = gaussian_projection_matrix(T₁, L, D)
    fs = [KDEComparator(X, t, W[l,:], boundary_offset = boundary_offset) for l in 1:L]
    W, fs
end


text:
using ..ITensorAutoHOOT
using ..ITensorAutoHOOT: generate_optimal_tree

# contract into one ITensor
ITensors.ITensor(t::TreeTensor; kwargs...) = contract(collect(t.tensors)...; kwargs...)

ITensors.contract(t1::TreeTensor; kwargs...) = t1

function ITensors.contract(t1::TreeTensor, t2::TreeTensor...; kwargs...)
  ts = mapreduce(t -> t.tensors, vcat, t2)
  return contract(t1, TreeTensor(ts...); kwargs...)
end

ITensors.contract(t_list::Vector{TreeTensor}; kwargs...) = contract(t_list...; kwargs...)

function ITensors.contract(t1::TreeTensor, t2::TreeTensor; cutoff, maxdim)
  # print("\n inputs are ", t1, t2, "\n")
  connect_inds = intersect(inds(t1), inds(t2))
  if length(connect_inds) <= 1
    return TreeTensor(t1.tensors..., t2.tensors...)
  end
  network = [t1.tensors..., t2.tensors...]
  if length(noncommoninds(network...)) <= 1
    return TreeTensor(contract(network...))
  end
  # TODO: add caching here
  uncontract_inds = noncommoninds(network...)
  inds_btree = inds_binary_tree(network, uncontract_inds; algorithm="mincut")
  tree = tree_approximation(network, inds_btree; cutoff=cutoff, maxdim=maxdim)
  out = TreeTensor(tree...)
  # print("output is ", out, "\n")
  return out
end

function uncontract_inds_binary_tree(path::Vector, uncontract_inds::Vector)
  inds1 = uncontract_inds_binary_tree(path[1], uncontract_inds)
  inds2 = uncontract_inds_binary_tree(path[2], uncontract_inds)
  if inds1 == []
    return inds2
  elseif inds2 == []
    return inds1
  else
    return [inds1, inds2]
  end
end

function uncontract_inds_binary_tree(tensor::ITensor, uncontract_inds::Vector)
  return intersect(inds(tensor), uncontract_inds)
end

function tree_approximation(
  network::Vector{ITensor}, inds_btree::Vector; cutoff=1e-15, maxdim=10000
)
  # inds_dict map each inds_btree node to generated indices
  # tnets_dict map each inds_btree node to a tensor network
  inds_dict = Dict()
  tnets_dict = Dict()
  projectors = []

  function insert_projectors(tree::Vector)
    insert_projectors(tree[1])
    insert_projectors(tree[2])
    # get the tensor network of the current tree node
    net1, net2 = tnets_dict[tree[1]], tnets_dict[tree[2]]
    net_bra = vcat(network, symdiff(net1, net2))
    # get the two indices to merge at this step
    ind1, ind2 = inds_dict[tree[1]], inds_dict[tree[2]]
    # form normal equations to the system to factorize
    net_ket = sim(linkinds, net_bra)
    sim_dict = Dict()
    for ind in [ind1..., ind2...]
      sim_dict[ind] = sim(ind)
    end
    net_ket = sim([ind1..., ind2...], net_ket, sim_dict)
    net_normal = vcat(net_bra, net_ket)
    tensor_normal = contract(generate_optimal_tree(net_normal))
    # use eig to factorize and get the projector
    rinds = (ind1..., ind2...)
    linds = map(i -> sim_dict[i], rinds)
    diag, U = eigen(
      tensor_normal, linds, rinds; cutoff=cutoff, maxdim=maxdim, ishermitian=true
    )
    dr = commonind(diag, U)
    # add the projector to the list projectors
    projectors = vcat(projectors, [U])
    tnets_dict[tree] = vcat(net_bra, [U])
    return inds_dict[tree] = (dr,)
  end

  function insert_projectors(tree::Union{Vector{<:Index},<:Index})
    inds_dict[tree] = Tuple(tree)
    return tnets_dict[tree] = network
  end

  @assert (length(inds_btree) >= 2)
  insert_projectors(inds_btree[1])
  insert_projectors(inds_btree[2])
  last_network = Vector{ITensor}(vcat(network, projectors))
  last_tensor = contract(generate_optimal_tree(last_network))
  return vcat(projectors, [last_tensor])
end


text:
using NeuralNets
using HDF5, JLD
using ProgressMeter

image_names = readdir("../images")

nTrain = 20 # 2000
N_stim = 100 # input dim
T = 300 # msec
dt = NeuralNets.dt
imgs = image_names[rand(1:length(image_names),nTrain)] # chosen images

net = buildnet()

progress = Progress(nTrain,1)
println("Starting Training Procedure with $nTrain iterations...")
for itr=1:nTrain
	img_name = string("../images/", imgs[itr])
	spike_train = image_stim(img_name, int(sqrt(N_stim)), T, dt)

	stim_spikes = cell(length(net.l))
	for li=1:length(net.l)
		if li in 1:3
			stim_spikes[li] = spike_train[:,:,li]
		else
			stim_spikes[li] = nothing
		end
	end
	rs = RunSettings(T,stim_spikes)
	resetLayers!(net)
	result = runSpikingNet!(net,rs)
	normalizeWeights!(net)

	if itr%5==0 || itr == 1
		fname = string("../data/t2_", itr, ".jld")
		println("Saving Data...")
		save(fname,"pre", net.pre, "runresult", result)
	end
	next!(progress)
end