text:
#-----------------------------------------------------------#
# Build optimization model for circuit decomsposition here  #
#-----------------------------------------------------------#
import JuMP: MOI

function build_QCModel(data::Dict{String, Any}; options = nothing)

    m_qc = QCO.QuantumCircuitModel(data)

    # Update defaults to user-defined options
    if options !== nothing
        for i in keys(options)
            QCO.set_option(m_qc, i, options[i])
        end
    end

    QCO._catch_options_errors(m_qc)

    # convex-hull formulation per depth, but larger number of variables and constraints
    if m_qc.options.model_type == "balas_formulation"

        QCO.variable_QCModel_balas(m_qc)
        QCO.constraint_QCModel_balas(m_qc)

    # minimal variables and constraints, but not a convex-hull formulation per depth
    elseif m_qc.options.model_type == "compact_formulation"

        QCO.variable_QCModel_compact(m_qc)
        QCO.constraint_QCModel_compact(m_qc)

    end

    QCO.objective_QCModel(m_qc)

    return m_qc
end

function variable_QCModel_balas(qcm::QuantumCircuitModel)
    QCO.variable_gates_per_depth(qcm)
    QCO.variable_gates_onoff(qcm)
    QCO.variable_sequential_gate_products(qcm)
    QCO.variable_gate_products_copy(qcm)
    QCO.variable_gate_products_linearization(qcm)

    if qcm.data["decomposition_type"] == "approximate"
        QCO.variable_slack_for_feasibility(qcm)
    end

    return
end

function constraint_QCModel_balas(qcm::QuantumCircuitModel)

    QCO.constraint_single_gate_per_depth(qcm)
    QCO.constraint_gates_onoff_per_depth(qcm)
    QCO.constraint_gate_initial_condition(qcm)
    QCO.constraint_gate_intermediate_products(qcm)
    QCO.constraint_gate_product_linearization(qcm)
    QCO.constraint_gate_target_condition(qcm)
    QCO.constraint_cnot_gate_bounds(qcm)
    (!qcm.data["are_gates_real"]) && (QCO.constraint_complex_to_real_symmetry(qcm))

    QCO.constraint_QCModel_valid(qcm)

    return
end

function variable_QCModel_compact(qcm::QuantumCircuitModel)

    QCO.variable_gates_onoff(qcm)
    QCO.variable_sequential_gate_products(qcm)
    QCO.variable_gate_products_linearization(qcm)

    if qcm.data["decomposition_type"] == "approximate"
        QCO.variable_slack_for_feasibility(qcm)
    end

    return
end

function constraint_QCModel_compact(qcm::QuantumCircuitModel)

    QCO.constraint_single_gate_per_depth(qcm)
    QCO.constraint_gate_initial_condition_compact(qcm)
    QCO.constraint_gate_intermediate_products_compact(qcm)
    QCO.constraint_gate_product_linearization(qcm)
    QCO.constraint_gate_target_condition_compact(qcm)
    QCO.constraint_cnot_gate_bounds(qcm)
    (!qcm.data["are_gates_real"]) && (QCO.constraint_complex_to_real_symmetry(qcm))

    QCO.constraint_QCModel_valid(qcm)

    return
end

function constraint_QCModel_valid(qcm::QuantumCircuitModel)

    if qcm.options.all_valid_constraints != -1

        if qcm.options.all_valid_constraints == 1
            QCO.constraint_commutative_gate_pairs(qcm)
            QCO.constraint_involutory_gates(qcm)
            QCO.constraint_redundant_gate_product_pairs(qcm)
            QCO.constraint_idempotent_gates(qcm)
            QCO.constraint_identity_gate_symmetry(qcm)
            QCO.constraint_convex_hull_complex_gates(qcm)
            QCO.constraint_complex_unit_magnitude(qcm)

        elseif qcm.options.all_valid_constraints == 0
            qcm.options.commute_gate_constraints            && QCO.constraint_commutative_gate_pairs(qcm)
            qcm.options.involutory_gate_constraints         && QCO.constraint_involutory_gates(qcm)
            qcm.options.redundant_gate_pair_constraints     && QCO.constraint_redundant_gate_product_pairs(qcm)
            qcm.options.idempotent_gate_constraints         && QCO.constraint_idempotent_gates(qcm)
            qcm.options.identity_gate_symmetry_constraints  && QCO.constraint_identity_gate_symmetry(qcm)
            qcm.options.convex_hull_gate_constraints        && QCO.constraint_convex_hull_complex_gates(qcm)
            qcm.options.unit_magnitude_constraints          && QCO.constraint_complex_unit_magnitude(qcm)
        end

    end

    return
end

function objective_QCModel(qcm::QuantumCircuitModel)

    if qcm.data["objective"] == "minimize_depth"
        QCO.objective_minimize_total_depth(qcm)

    elseif qcm.data["objective"] == "minimize_cnot"
        QCO.objective_minimize_cnot_gates(qcm)
    end

    return
end

""
function optimize_QCModel!(qcm::QuantumCircuitModel; optimizer=nothing)

    if qcm.options.relax_integrality
        JuMP.relax_integrality(qcm.model)
    end

    if JuMP.mode(qcm.model) != JuMP.DIRECT && optimizer !== nothing
        if JuMP.backend(qcm.model).optimizer === nothing
            JuMP.set_optimizer(qcm.model, optimizer)
        else
            Memento.warn(_LOGGER, "Model already contains optimizer, cannot use optimizer specified in `optimize_QCModel!`")
        end
    end

    JuMP.set_time_limit_sec(qcm.model, qcm.options.time_limit)

    if !qcm.options.optimizer_log
        JuMP.set_silent(qcm.model)
    end

    if JuMP.mode(qcm.model) != JuMP.DIRECT && JuMP.backend(qcm.model).optimizer === nothing
        Memento.error(_LOGGER, "No optimizer specified in `optimize_QCModel!` or the given JuMP model.")
    end

    start_time = time()

    _, solve_time, solve_bytes_alloc, sec_in_gc = @timed JuMP.optimize!(qcm.model)

    try
        solve_time = JuMP.solve_time(qcm.model)
    catch
        Memento.warn(_LOGGER, "The given optimizer does not provide the SolveTime() attribute, falling back on @timed.  This is not a rigorous timing value.");
    end

    Memento.debug(_LOGGER, "JuMP model optimize time: $(time() - start_time)")

    qcm.result = QCO.build_QCModel_result(qcm, solve_time)

    return qcm.result
end

function run_QCModel(params::Dict{String, Any},
                     qcm_optimizer::MOI.OptimizerWithAttributes;
                     options = nothing)

    data = QCO.get_data(params)

    model_qc = QCO.build_QCModel(data, options = options)

    result_qc = QCO.optimize_QCModel!(model_qc, optimizer = qcm_optimizer)

    if model_qc.options.relax_integrality
        if result_qc["primal_status"] == MOI.FEASIBLE_POINT
            Memento.info(_LOGGER, "Integrality-relaxed solutions can be found in the results dictionary")
        else
            Memento.info(_LOGGER, "Infeasible primal status for the integrality-relaxed problem")
        end
    else
        if model_qc.options.visualize_solution
            QCO.visualize_solution(result_qc, data)
        end
    end

    return result_qc
end

function set_option(qcm::QuantumCircuitModel, s::Symbol, val)
    Base.setproperty!(qcm.options, s, val)
end

function _catch_options_errors(qcm::QuantumCircuitModel)

    if !(qcm.options.model_type in ["compact_formulation", "balas_formulation"])
        Memento.warn(_LOGGER, "Invalid model_type. Setting it to default value (compact_formulation).")
        QCO.set_option(qcm, :model_type, "compact_formulation")
    end

    if !(qcm.options.all_valid_constraints in [-1,0,1])
        Memento.warn(_LOGGER, "Invalid all_valid_constraints; choose a value ∈ [-1,0,1]. Setting it to default value of 0.")
        QCO.set_option(qcm, :all_valid_constraints, 0)
    end

end


text:
####################################
#
# options
#
####################################

struct CDOptions
  maxIter::Int64
  optTol::Float64
  randomize::Bool           # coordinate are visitied in a randomized order or not
  warmStart::Bool           # when running CD, should we do path following or not
  numSteps::Int64           # when pathFollowing, how many points are there on the path
end

CDOptions(;
  maxIter::Int64=2000,
  optTol::Float64=1e-7,
  randomize::Bool=true,
  warmStart::Bool=true,
  numSteps::Int=50) = CDOptions(maxIter, optTol, randomize, warmStart, numSteps)

# this is used for ScaledLasso and FeasibleLasso
struct IterLassoOptions
  maxIter::Int64
  optTol::Float64
  initProcedure::Symbol    # :Screening, :InitStd, :WarmStart
  sinit::Int64             # how many columns of X to use to estimate initial variance or obtain initial residuals
  σinit::Float64
  optionsCD::CDOptions
end

IterLassoOptions(;
  maxIter::Int64=20,
  optTol::Float64=1e-2,
  initProcedure::Symbol=:Screening,
  sinit::Int64=5,
  σinit::Float64=1.,
  optionsCD::CDOptions=CDOptions()) = IterLassoOptions(maxIter, optTol, initProcedure, sinit, σinit, optionsCD)

####################################
#
# helper functions
#
####################################

"""
Helper function that finds an initial estimate for σ that is needed
for Lasso and ScaledLasso procedures.

The procedure works as follows:

* s input variables that are most correlated with the response y are found
* y is regressed on those s features
* σ is estimated based on the residuals, which gives an upper bound on the true sigma
"""
_findInitSigma!(
  X::AbstractMatrix{T},
  y::AbstractVector{T},
  s::Int,
  storage::Vector{T}) where {T <: AbstractFloat} =  std(_findInitResiduals!(X, y, s, storage))

function _findInitResiduals!(
  X::AbstractMatrix{T},
  y::AbstractVector{T},
  s::Int,
  storage::Vector{T}) where {T <: AbstractFloat}

  S = _findLargestCorrelations(X, y, s)
  Xs = view(X, :, S)
  mul!(storage, Xs, Xs \ y)
  @. storage = y - storage
  return storage
end

function _findInitResiduals!(
  w::AbstractVector{T},
  X::AbstractMatrix{T},
  y::AbstractVector{T},
  s::Int,
  storage::Vector{T}) where {T <: AbstractFloat}

  S = _findLargestCorrelations(w, X, y, s)

  Xs = view(X, :, S)
  mul!(storage, Xs, (Xs' * Diagonal(w) * Xs) \ (Xs' * Diagonal(w) * y))
  @. storage = y - storage
  return storage
end

# return a bit array containing indices of columns
function _findLargestCorrelations(
  X::AbstractMatrix{T},
  y::AbstractVector{T},
  s::Int) where {T <: AbstractFloat}

  p = size(X, 2)
  storage = Array{T}(undef, p)
  mul!(storage, transpose(X), y)
  @. storage = abs(storage)
  S = storage .>= nlargest(s, storage)[end]
end

function _findLargestCorrelations(
  w::AbstractVector{T},
  X::AbstractMatrix{T},
  y::AbstractVector{T},
  s::Int) where {T <: AbstractFloat}

  n, p = size(X)
  storage = Array{T}(undef, p)
  @inbounds for j=1:p
    val = zero(T)
    @simd for i=1:n
      val += X[i,j] * w[i] * y[i]
    end
    storage[j] = abs(val)
  end
  S = storage .>= nlargest(s, storage)[end]
end

function _stdX!(out::Vector{T}, X::AbstractMatrix{T}) where {T <: AbstractFloat}
  n, p = size(X)

  @inbounds for j=1:p
    v = zero(T)
    @simd for i=1:n
      v += X[i, j]^2.
    end
    out[j] = sqrt(v / n)
  end
  out
end

function _stdX!(out::Vector{T}, w::AbstractVector{T}, X::AbstractMatrix{T}) where {T <: AbstractFloat}
  n, p = size(X)

  @inbounds for j=1:p
    v = zero(T)
    @simd for i=1:n
      v += w[i] * X[i, j]^2.
    end
    out[j] = sqrt(v / n)
  end
  out
end

function _getLoadings!(out::Vector{T}, X::AbstractMatrix{T}, e::AbstractVector{T}) where {T <: AbstractFloat}
  n, p = size(X)

  @inbounds for j=1:p
    v = zero(T)
    @simd for i=1:n
      v += (X[i, j]*e[i])^2.
    end
    out[j] = sqrt(v / n)
  end
  out
end

function _getSigma(w::AbstractVector{T}, r::AbstractVector{T}) where {T <: AbstractFloat}
    n = length(w)
    σ = zero(T)
    for ii = 1:n
        σ += r[ii]^2 * w[ii]
    end
    σ /= sum(w)
    sqrt(σ)
end

ext: jl


text:
using VaxData, Test, InteractiveUtils

include("vaxints.jl")
include("vaxfloatf.jl")
include("vaxfloatd.jl")
include("vaxfloatg.jl")


text:
# this file is to run JWAS
using DataFrames,CSV,Random,Distributions,JWAS,DelimitedFiles, LinearAlgebra
Random.seed!(123)

method            = ARGS[1]
rep               = parse(Int, ARGS[2])  #rep1,...,rep20
pct               = parse(Float64, ARGS[3])

@show method,rep,pct

mainpath="/group/qtlchenggrp/tianjing/omics_fix_all/missing_pattern_legarra/"
cd(mainpath*"$pct/$method/rep$rep/")

############ READ DATA ##########
#read missing index
#some individuals have no omics data
missingIndex=Int.(vec(readdlm(mainpath*"missing_index/missingIndex.pct$pct.rep$rep.txt")))
@show missingIndex[1:3]
@show length(missingIndex)

#read phenotype
phenofile="/group/qtlchenggrp/tianjing/omics_data/andrew1055/data1/y_no_ar.csv"
phenotypes = CSV.read(phenofile,DataFrame,delim = ',',header=true,missingstrings=["NA"])

#read omics
omicsfile  = "/group/qtlchenggrp/tianjing/omics_data/andrew1055/data1/omics.csv"
omics  = CSV.read(omicsfile,DataFrame,delim = ',',header=true,missingstrings=["NA"]);
insertcols!(omics,2,:y => phenotypes[:,:y]); #insert phenotype into 2nd column of omics
for i in 2:size(omics,2) #skip ID
    omics[!,i]= convert(Vector{Union{Missing,Float64}}, omics[!,i]); #convert the data type to allow missing data
end
omics[956:1055,:y] .= missing  #missing phenotype for individuals in testing datasett
omics[missingIndex,3:end].=missing; #some individuals have no omics data
@show [count(ismissing,col) for col in eachcol(omics)]

omics_name=["gene$i" for i in 1:1200]

@show size(omics)
@show omics[1:5,1:5];
@show omics[end-4:end,1:5];

#fixed variances components required in Christensenet al.(2021)
if occursin("gblup",method) #NN-GBLUP
	G=Matrix(Diagonal(repeat([1.22],1200)))  #genetic variance for each omics is 1.22 in Christensenet al.(2021)
	R=Matrix(Diagonal(repeat([0.78],1200)))  #residual variance for omics is 0.78 in Christensenet al.(2021)
	genofile   = "/group/qtlchenggrp/tianjing/omics_data/andrew1055/data1/GRM_geno.csv"
	genotypes  = get_genotypes(genofile,G;separator=',',method="GBLUP",header=false,G_is_marker_variance = false,estimateVariance=false);
end

# build model
model_equation = "y = intercept + genotypes";
#fixed variances components required in Christensenet al.(2021)
user_σ2_yobs=19.37777777777778 #residual variance for observed phenotype
user_σ2_weightsNN=0.01         #variance for neural network weight between omics and observed phenotype

if occursin("gblup",method)   #NN-GBLUP
	model = build_model(model_equation,R;
                        num_hidden_nodes=1200,
                        latent_traits=omics_name,
                        nonlinear_function="linear");
	out   = runMCMC(model,omics,chain_length=5000,printout_model_info=false,user_σ2_yobs=user_σ2_yobs,user_σ2_weightsNN=user_σ2_weightsNN,estimate_variance=false);
end

# Step 5: Check Accuruacy
#read true bv
bv=vec(readdlm("/group/qtlchenggrp/tianjing/omics_data/andrew1055/data1/tbv_m.txt"));
accuruacy_all  = cor(out["EBV_NonLinear"][!,:EBV],bv)
accuruacy_test  = cor(out["EBV_NonLinear"][!,:EBV][956:1055],bv[956:1055])

@show pct,method,rep
@show accuruacy_all
@show accuruacy_test

res=[rep,accuruacy_test, accuruacy_all]
writedlm("res$rep.txt",res)

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
module SnoopCompile

using SnoopCompileCore
export @snoopc
isdefined(SnoopCompileCore, Symbol("@snoopi")) && export @snoopi
if isdefined(SnoopCompileCore, Symbol("@snoopr"))
    export @snoopr, uinvalidated, invalidation_trees, filtermod, findcaller, ascend
end

using Core: MethodInstance, CodeInfo
using Serialization, OrderedCollections

# Parcel Regex
const anonrex = r"#{1,2}\d+#{1,2}\d+"         # detect anonymous functions
const kwrex = r"^#kw##(.*)$|^#([^#]*)##kw$"   # detect keyword-supplying functions
const kwbodyrex = r"^##(\w[^#]*)#\d+"         # detect keyword body methods
const genrex = r"^##s\d+#\d+$"                # detect generators for @generated functions
const innerrex = r"^#[^#]+#\d+"               # detect inner functions

# Parcel
include("parcel_snoopc.jl")

if VERSION >= v"1.2.0-DEV.573"
    include("parcel_snoopi.jl")
end

if isdefined(SnoopCompileCore, Symbol("@snoopr"))
    include("invalidations.jl")
end

# Write
include("write.jl")

end # module


text:
# include("../../distributions/src/distributions_module.jl")
include("../ode_test_functions/modelling_module.jl")
include("../src/ode_module.jl")
include("../ode_test_functions/ode_test_functions.jl")

module odetest
using ..modelling
using ..ode
using ..odesample
using ..numder
# using ..distributions
using PyPlot
PyPlot.pygui(true)
using BenchmarkTools

function force_fcn(t)
    t_start = 1.5
    t_end  = 3.5 # 0.8
    F_amplitude = 30
    return F_amplitude * (modelling.step(t - t_start) - modelling.step(t - t_end))
end

# implement cash karp
# look for rosenbrock23

function test_rk23()

    omega = 30
    D = 0.1
    m = 1.0
    c = m * omega^2
    d = 2 * D * omega
    params = odesample.HarmonicOscillatorOptions(;
        m = m,
        d = d,
        c = c,
        f = force_fcn)
    # ode_fcn(der_q, t, q) = odesample.harmonic_oscillator!(der_q, t, q, params)

    params = odesample.VanDerPolOptions(1_000.0)
    ode_fcn(der_q, t, q) = odesample.van_der_pol!(der_q, t, q, params)

    # ode_fcn = odesample.stiff_test_v1

    time_interval = ode.TimeInterval(5.0) # 50

    # time_interval = ode.TimeInterval(50.0) # 50
    time_interval = ode.TimeInterval(3_000.0) # 50
    init_cond     = [0.0, 1.0]
    # init_cond     = [1.0, 1.0, 0.0]

    ode_problem = ode.OdeProblem(
        ode_fcn,
        time_interval,
        init_cond)

    # ode_solver = ode.BogaczkiShampine23()
    # ode_solver = ode.DormandPrince45()
    ode_solver = ode.Tsitouras45()
    # ode_solver = ode.Ros23()

    basic_options = ode.BasicOdeOptions(
        ode_problem;
        abs_tol = 1e-9,
        rel_tol = 1e-9,
        max_iter = 10_000_000
    )

    step_size_controller = ode.ClassicalStepSizeControl()
    # step_size_controller = ode.PIStepSizeControl()

    stepper_options = ode.StepperOptions(;
        max_step_size = 5.0)

    @time (ode_res, step_stat) = ode.solve(
        ode_problem,
        ode_solver;
        basic_options = basic_options,
        step_size_controller = step_size_controller,
        stepper_options = stepper_options
        ) #

    t = ode.time(ode_res)
    q = ode.variables(ode_res)

# #=
    PyPlot.figure()

    a1 = PyPlot.subplot(2,1,1)
    PyPlot.grid()
    PyPlot.plot( t, q[1, :] )

    PyPlot.subplot(2,1,2, sharex = a1)
    PyPlot.grid()
    PyPlot.plot( t, q[2, :] )

    dt = diff(t)
    tmids = 0.5 * (t[1:end-1] + t[2:end])

    PyPlot.figure()
    PyPlot.plot(tmids, dt)
    PyPlot.grid()
# =#

    return step_stat
end

# test_rk23()

function test()
    x  = [1.0, 0.0]
    k2 = zeros(2)
    f_tmp = zeros(2)
    h_tmp = zeros(2)
    t = 10.0

    params = odesample.VanDerPolOptions(10.0)
    ode_fcn(der_q, t, q) = odesample.van_der_pol!(der_q, t, q, params)

    f_q(q) = ode_fcn(k2, t, q) # k2 is just a dummy

    L = length(x)
    J = zeros(L, L)
    # J      = numder.jacobian_fw!(J, f_q, x)
    J      = ode.jacobian_fw!(f_q, x, J, h_tmp, f_tmp)
    f_t(t) = ode_fcn(k2, t, x) # k2 is just a dummy
    # dfdt   = numder.numdiff_fw(f_t, t)

end

# test()

# x  = [1.0, 0.0]
# k2 = zeros(2)
# t = 10.0

# params = odesample.VanDerPolOptions(10.0)
# ode_fcn(der_q, t, q) = odesample.van_der_pol!(der_q, t, q, params)

# f_q(q) = ode_fcn(k2, t, q) # k2 is just a dummy

end # odetest

#=
module old_odetest

function force_fcn(t)
    t_start = 1.5
    t_end  = 3.5 # 0.8
    F_amplitude = 30
    return F_amplitude * (modelling.step(t - t_start) - modelling.step(t - t_end))
end

struct HarmonicalOscillatorOptions
    m
    D
    omega
    F
end
function HarmonicalOscillatorOptions(; m = 1, D = 0.1, omega = 3, F = force_fcn)
    HarmonicalOscillatorOptions(
        m,
        D,
        omega,
        F)
end

function simple_ode_fcn(der_q, t, q, options::HarmonicalOscillatorOptions = HarmonicalOscillatorOptions() )

    der_x = q[2]
    der_v =  options.F(t) / options.m  - options.omega^2 * q[1] - 2 * options.D * options.omega * q[2]

    der_q[1] = der_x
    der_q[2] = der_v

    return der_q

end

function test_exp_euler()

    time_interval   = ode.TimeInterval(5.0)
    initial_values  = [0.0, 1.0]
    options         = ode.OdeOptions(step_size = 1e-3)

    params = HarmonicalOscillatorOptions(; omega = 30, D = 0.1, F = force_fcn)
    fcn(t, q) = simple_ode_fcn(t, q, params)

    ode_res = ode.ode_solver_explicit_euler(fcn, time_interval, initial_values, options)

    PyPlot.figure()

    PyPlot.subplot(2,1,1)
    PyPlot.grid()
    PyPlot.plot( ode_res.t, ode_res.q[:, 1] )

    PyPlot.subplot(2,1,2)
    PyPlot.grid()
    PyPlot.plot( ode_res.t, ode_res.q[:, 2] )

end

function test_semiimp_euler()

    time_interval   = ode.TimeInterval(5.0)
    initial_values  = [0.0, 1.0]
    options         = ode.OdeOptions(step_size = 1e-3)

    params = HarmonicalOscillatorOptions(; omega = 50, D = 0.1, F = force_fcn)
    fcn(t, q) = simple_ode_fcn(t, q, params)

    ode_res = ode.ode_solver_semiimplicit_euler(fcn, time_interval, initial_values, options)

    PyPlot.figure()

    PyPlot.subplot(2,1,1)
    PyPlot.grid()
    PyPlot.plot( ode_res.t, ode_res.q[:, 1] )

    PyPlot.subplot(2,1,2)
    PyPlot.grid()
    PyPlot.plot( ode_res.t, ode_res.q[:, 2] )

end

function test_imp_euler()

    time_interval   = ode.TimeInterval(5.0)
    initial_values  = [0.0, 1.0]
    options         = ode.OdeOptions(step_size = 1e-3)

    params = HarmonicalOscillatorOptions(; omega = 30, D = 0.1, F = force_fcn)
    fcn(t, q) = simple_ode_fcn(t, q, params)

    ode_res = ode.ode_solver_implicit_euler(fcn, time_interval, initial_values, options)

    PyPlot.figure()

    a1 = PyPlot.subplot(2,1,1)
    PyPlot.grid()
    PyPlot.plot( ode_res.t, ode_res.q[:, 1] )

    PyPlot.subplot(2,1,2, sharex = a1)
    PyPlot.grid()
    PyPlot.plot( ode_res.t, ode_res.q[:, 2] )

end

function add_noise()

    no_force_fcn(t) = 0.0

    time_interval   = ode.TimeInterval(5.0)
    initial_values  = [0.0, 1.0]
    options         = ode.OdeOptions(step_size = 1e-3)

    params = HarmonicalOscillatorOptions(; omega = 30, D = 0.1, F = no_force_fcn)
    fcn(t, q) = simple_ode_fcn(t, q, params)

    ode_res = ode.ode_solver_implicit_euler(fcn, time_interval, initial_values, options)

    # add noise
    n_points = length(ode_res.t)

    (min_x, max_x) = extrema(ode_res.q[:,1])
    (min_v, max_v) = extrema(ode_res.q[:,2])

    x_noise_distr = distributions.NormalDistribution(0.0, 1.0)
    v_noise_distr = distributions.NormalDistribution(0.0, 1.0)

    x_mult_noise_level = 0.05
    v_mult_noise_level = 0.05

    x_add_noise_level = 0.05 * (max_x - min_x)
    v_add_noise_level = 0.05 * (max_v - min_v)

    x_mult_noise = distributions.random(x_noise_distr, n_points)
    v_mult_noise = distributions.random(v_noise_distr, n_points)
    x_add_noise  = distributions.random(x_noise_distr, n_points)
    v_add_noise  = distributions.random(v_noise_distr, n_points)

    t_orig = copy(ode_res.t)
    q_mod = similar(ode_res.q)
    q_mod[:,1] = ode_res.q[:,1] .* (1 .+ x_mult_noise_level .* x_mult_noise) .+ x_add_noise_level .* x_add_noise
    q_mod[:,2] = ode_res.q[:,2] .* (1 .+ v_mult_noise_level .* v_mult_noise) .+ v_add_noise_level .* v_add_noise
    mod_ode_res = ode.OdeResult(t_orig, q_mod)

    PyPlot.figure()

    a1 = PyPlot.subplot(2,1,1)
    PyPlot.grid()
    PyPlot.plot( mod_ode_res.t, mod_ode_res.q[:, 1] )
    PyPlot.plot( ode_res.t, ode_res.q[:, 1], linewidth = 2 )

    PyPlot.subplot(2,1,2, sharex = a1)
    PyPlot.grid()
    PyPlot.plot( mod_ode_res.t, mod_ode_res.q[:, 2] )
    PyPlot.plot( ode_res.t, ode_res.q[:, 2], linewidth = 2 )

end

end

=#

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using OrdinaryDiffEq, Flux, Optim, Random, Plots
using Zygote
using ForwardDiff
using LinearAlgebra, Statistics
using ProgressBars, Printf
using Flux.Optimise: update!, ExpDecay
using Flux.Losses: mae, mse
using BSON: @save, @load
using Catalyst

Random.seed!(1234);

# TODO: use YAML file for configuration
# Argments
is_restart = false;
p_cutoff = 0.0;
n_epoch = 1000000;
n_plot = 100;
opt = ADAMW(0.0001, (0.9, 0.999), 0.0);
datasize = 100;
tstep = 0.1;
n_exp_train = 20;
n_exp_test = 10;
n_exp = n_exp_train + n_exp_test;
noise = 1.f-3;
ns = 5;
nr = 10;
k = [];
alg = Tsit5();
atol = 1e-5;
rtol = 1e-2;

maxiters = 10000;

lb = 1.f-5;

rn = @reaction_network begin
    (1.0, 1.0), A ↔ B
    (1.0, 1.0), B ↔ C
    (1.0, 1.0), C ↔ D
    # (1.0, 1.0), 2A ↔ B + C
    # (1.0, 1.0), 2B ↔ C + D
    (1.0, 1.0), 2C ↔ D + E
end

# Generate data sets
u0_list = rand(Float32, (n_exp, ns));
u0_list[:, 1:2] .+= 0.2;
# u0_list[:, 3:end] .= 0.0;
tspan = Float32[0.0, datasize * tstep];
tsteps = range(tspan[1], tspan[2], length=datasize);
ode_data_list = zeros(Float32, (n_exp, ns, datasize));
std_list = [];

function max_min(ode_data)
    return maximum(ode_data, dims=2) .- minimum(ode_data, dims=2) .+ lb
end

for i in 1:n_exp
    u0 = u0_list[i, :];
    prob_trueode = ODEProblem(rn, u0, tspan, k);
    ode_data = Array(solve(prob_trueode, alg, saveat=tsteps));
    ode_data += randn(size(ode_data)) .* ode_data .* noise
    ode_data_list[i, :, :] = ode_data
    push!(std_list, max_min(ode_data));
end
y_std = maximum(hcat(std_list...), dims=2);

p = randn(Float32, nr * (ns + 1)) .* 0.5;

function p2vec(p)
    w_kf = p[1:nr];
    w_kb = w_kf;  # assume Kc = 1
    w_out = reshape(p[nr + 1:end], ns, nr);
    w_out = clamp.(w_out, -2.5, 2.5)
    return w_kf, w_kb, w_out
end

function crnn!(du, u, p, t)
    w_in_f = clamp.(-w_out, 0, 2.5);
    w_in_b = clamp.(w_out, 0, 2.5);

    u_in = @. log(clamp(u, lb, Inf))
    w_in_x_f = w_in_f' * u_in;
    w_in_x_b = w_in_b' * u_in;

    du .= w_out * @. (exp(w_in_x_f + w_kf) - exp(w_in_x_b + w_kb));
end

u0 = u0_list[1, :]
prob = ODEProblem(crnn!, u0, tspan, saveat=tsteps,
                  atol=atol, rtol=rtol)

function pred_ode(u0, p)
    global w_kf, w_kb, w_out
    w_kf, w_kb, w_out = p2vec(p)
    pred = Array(solve(prob, alg, u0=u0, p=p; maxiters=maxiters))
    return pred
end
pred_ode(u0, p);

function display_p(p)
    local w_kf, w_kb, w_out
    w_kf, w_kb, w_out = p2vec(p);

    println("species (column) reaction (row)")

    println("\n w_out | w_kf | w_kb")
    println(rn.states)
    show(stdout, "text/plain", round.(vcat(w_out, w_kf', w_kb')', digits=3))

    println("\n")
end
display_p(p)

function loss_ode(p, i_exp)
    pred = pred_ode(u0_list[i_exp, :], p)
    loss = mae(ode_data_list[i_exp, :, :] ./ y_std, pred ./ y_std)
    return loss
end
loss_ode(p, 1)

# Callback function to observe training

species = rn.states;
cbi = function (p, i_exp)
    ode_data = ode_data_list[i_exp, :, :]
    pred = pred_ode(u0_list[i_exp, :], p)
    l_plt = []
    for i in 1:ns
        plt = Plots.scatter(tsteps, ode_data[i,:],
                      markercolor=:transparent,
                      label="Exp",
                      framestyle=:box)
        plot!(plt, tsteps, pred[i,:], lw=3, label="CRNN")
        plot!(xlabel="Time", ylabel=species[i])

        if i == 1
            plot!(plt, legend=true, framealpha=0)
        else
        plot!(plt, legend=false)
        end

        push!(l_plt, plt)
    end
    plt_all = plot(l_plt...)
    png(plt_all, string("figs/i_exp_", i_exp))

    println("plot for exp $i_exp is updated \n")
    return false
end

l_loss_train = []
l_loss_val = []
iter = 1
cb = function (p, loss_train, loss_val)

    global l_loss_train, l_loss_val, iter
    push!(l_loss_train, loss_train)
    push!(l_loss_val, loss_val)

    if iter % n_plot == 0
        display_p(p)

        @printf("min loss train %.4e val %.4e\n", minimum(l_loss_train), minimum(l_loss_val))

        l_exp = randperm(n_exp)[1:1];
        println("update plot for ", l_exp)
        for i_exp in l_exp
            cbi(p, i_exp)
        end

        plt_loss = plot(l_loss_train, xscale=:log10, yscale=:log10, label="train");
        plot!(plt_loss, l_loss_val, label="val");

        png(plt_loss, "figs/loss");

        @save "./checkpoint/mymodel.bson" p opt l_loss_train l_loss_val iter
    end

iter += 1;
end

if is_restart
    @load "./checkpoint/mymodel.bson" p opt l_loss_train l_loss_val iter;
    iter += 1;
end

i_exp = 1
epochs = ProgressBar(iter:n_epoch)
loss_epoch = zeros(Float32, n_exp);
for epoch in epochs
    global p
    for i_exp in randperm(n_exp_train)
        grad = ForwardDiff.gradient(x -> loss_ode(x, i_exp), p)
        update!(opt, p, grad)
    end
    for i_exp in 1:n_exp
        loss_epoch[i_exp] = loss_ode(p, i_exp)
    end
    loss_train = mean(loss_epoch[1:n_exp_train]);
    loss_val = mean(loss_epoch[n_exp_train + 1:end]);
    set_description(epochs, string(@sprintf("Loss train %.4e val %.4e", loss_train, loss_val)))
    cb(p, loss_train, loss_val);
end

max_stars_repo_path: case1 rev/case1.jl

max_stars_repo_name: DENG-MIT/CRNN

max_issues_repo_path: case1 rev/case1.jl

max_issues_repo_name: DENG-MIT/CRNN

max_forks_repo_path: case1 rev/case1.jl

max_forks_repo_name: DENG-MIT/CRNN


text:
"""
Provides functionality for working with HTTP headers in Genie.
"""
module Headers

import Revise, HTTP
import Genie

"""
    set_headers!(req::HTTP.Request, res::HTTP.Response, app_response::HTTP.Response) :: HTTP.Response

Configures the response headers.
"""
function set_headers!(req::HTTP.Request, res::HTTP.Response, app_response::HTTP.Response) :: HTTP.Response
  if req.method == Genie.Router.OPTIONS || req.method == Genie.Router.GET || req.method == Genie.Router.POST
    request_origin = get(Dict(req.headers), "Origin", "")

    allowed_origin_dict = Dict("Access-Control-Allow-Origin" =>
      in(request_origin, Genie.config.cors_allowed_origins) || in("*", Genie.config.cors_allowed_origins)
      ? request_origin
      : strip(Genie.config.cors_headers["Access-Control-Allow-Origin"])
    )

    app_response.headers = [d for d in merge(Genie.config.cors_headers, allowed_origin_dict, Dict(res.headers), Dict(app_response.headers))]
  end

  app_response.headers = [d for d in
                            merge(
                              Dict(res.headers),
                              Dict(app_response.headers),
                              Dict("Server" => Genie.config.server_signature)
                            )
                          ]

  app_response
end

"""
    normalize_headers(req::HTTP.Request)

Makes request headers case insensitive.
"""
function normalize_headers(req::HTTP.Request) :: HTTP.Request
  headers = Dict(req.headers)
  normalized_headers = Dict{String,String}()

  for (k,v) in headers
    normalized_headers[normalize_header_key(string(k))] = string(v)
  end

  req.headers = [k for k in normalized_headers]

  req
end

"""
    normalize_header_key(key::String) :: String

Brings header keys to standard casing.
"""
function normalize_header_key(key::String) :: String
  join(map(x -> uppercasefirst(lowercase(x)), split(key, '-')), '-')
end

end


text:
const TREES = readlines("day03.txt")

function numhit(slope::Tuple{Int64,Int64}, trees::Array{<:AbstractString,1})
    total = 0
    height = length(trees)
    width = length(trees[1])
    for i = 2:slope[1]:height
        j = mod1(slope[2] * (i - 1) + 1, width)
        if trees[i][j] == '#'
            total += 1
        end
    end
    return total
end

println("Part 1: $(numhit((1, 3), TREES))")

const SLOPES = [(1, 1), (1, 3), (1, 5), (1, 7), (2, 1)]

println("Part 2: $(prod(numhit(slope, TREES) for slope in SLOPES))")


text:
using PartonDensity
using Random, Distributions
using Test

@testset "Valence PDF parametrisation" begin

    local val_pdf_params

    for i in 1:100

        λ_u = rand(Uniform(0, 1))
        K_u = rand(Uniform(2, 10))
        λ_d = rand(Uniform(0, 1))
        K_d = rand(Uniform(2, 10))
        λ_g1 = rand(Uniform(0, 1))
        λ_g2 = rand(Uniform(-1, 0))
        K_g = rand(Uniform(2, 10))
        λ_q = rand(Uniform(-1, 0))
        weights = ones(7)

        val_pdf_params = ValencePDFParams(λ_u=λ_u, K_u=K_u,
                                          λ_d=λ_d, K_d=K_d, λ_g1=λ_g2,
                                          λ_g2=λ_g2, K_g=K_g, λ_q=λ_q,
                                          weights=weights)

        @test int_xtotx(val_pdf_params) ≈ 1.0

    end

    @test typeof(val_pdf_params) == ValencePDFParams

    @test val_pdf_params.param_type == VALENCE_TYPE

end

@testset "Dirichlet PDF parametrisation" begin

    local dir_pdf_params

    for i in 1:100

        K_u = rand(Uniform(2, 10))
        K_d = rand(Uniform(2, 10))
        λ_g1 = rand(Uniform(0, 1))
        λ_g2 = rand(Uniform(-1, 0))
        K_g = rand(Uniform(2, 10))
        λ_q = rand(Uniform(-1, 0))
        weights = ones(9)

        dir_pdf_params = DirichletPDFParams(K_u=K_u, K_d=K_d, λ_g1=λ_g2,
                                            λ_g2=λ_g2, K_g=K_g, λ_q=λ_q,
                                            weights=weights)

        @test int_xtotx(dir_pdf_params) ≈ 1.0

    end

    @test typeof(dir_pdf_params) == DirichletPDFParams

    @test dir_pdf_params.param_type == DIRICHLET_TYPE

end

@testset "Forward model" begin

    # Define different parametriations for testing
    val_pdf_params = ValencePDFParams(λ_u=0.6, K_u=3.4,
                                      λ_d=0.7, K_d=4.7,
                                      λ_g1=0.4, λ_g2=-0.6,
                                      K_g=4.2, λ_q=-0.2,
                                      weights=[5., 5., 1., 1., 1., 0.5, 0.5])

    dir_pdf_params = DirichletPDFParams(K_u=3.4, K_d=4.7,
                                        λ_g1=0.4, λ_g2=-0.6,
                                        K_g=4.2, λ_q=-0.2,
                                        weights=[3., 1., 5., 5., 1., 1., 1., 0.5, 0.5])

    pdf_params_list = [val_pdf_params, dir_pdf_params]

    # Initialise
    qcdnum_grid = QCDNUMGrid(x_min=[1.0e-3], x_weights=[1], nx=100,
                             qq_bounds=[1.0e2, 3.0e4], qq_weights=[1.0, 1.0],
                             nq=50, spline_interp=3)

    qcdnum_params = QCDNUMParameters(order=2, α_S=0.118, q0=100.0,
                                     grid=qcdnum_grid, n_fixed_flav=5,
                                     iqc=1, iqb=1, iqt=1, weight_type=1)

    splint_params = SPLINTParameters(nuser=1000)
    quark_coeffs = QuarkCoefficients()

    forward_model_init(qcdnum_grid, qcdnum_params, splint_params)

    # Run forward model
    for pdf_params in pdf_params_list

        counts_pred_ep, counts_pred_em = forward_model(pdf_params, qcdnum_params,
                                                       splint_params, quark_coeffs)

        @test all(counts_pred_ep .>= 0.0)
        @test all(counts_pred_ep .<= 1.0e3)

        @test all(counts_pred_em .> 0.0)
        @test all(counts_pred_em .<= 1.0e3)

        nbins = size(counts_pred_ep)[1]
        counts_obs_ep = zeros(UInt64, nbins)
        counts_obs_em = zeros(UInt64, nbins)

        for i in 1:nbins
            counts_obs_ep[i] = rand(Poisson(counts_pred_ep[i]))
            counts_obs_em[i] = rand(Poisson(counts_pred_em[i]))
        end

        sim_data = Dict{String, Any}()
        sim_data["nbins"] = nbins
        sim_data["counts_obs_ep"] = counts_obs_ep
        sim_data["counts_obs_em"] = counts_obs_em

        mktempdir() do tmp_dir

            output_file = joinpath(tmp_dir, "test_sim.h5")

            pd_write_sim(output_file, pdf_params, sim_data)

            new_pdf_params, new_sim_data = pd_read_sim(output_file)

            @test typeof(new_pdf_params) == typeof(pdf_params)

            @test new_sim_data == new_sim_data

        end

    end

end


text:
"
Implementation of the QuadEdge-Datastructure. [Guibas85]_

.. [Guibas85] Leonidas Guibas and Jorge Stolfi, Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams, *ACM Transactions on Graphics*, 4(2), 1985, 75-123.
"
@enum ConstraintType NoConstraint OnLeftSide OnRightSide

type Edge
  origin::Int
  next::Int
  constraintType::ConstraintType
end

type SubDivision{T}
  edges::Vector{Edge}
  vertexCache::Vector{Int}
  vertices::Vector{T}
  faces::Vector{T}
end

typealias VertexIndex Int
typealias EdgeIndex Int

"""
    makeedge!(sd::SubDivision)

Creates an empty edge and adds it to given SubDivision.

#returns
* Number of newly-added edge.
"""
function makeedge!(sd::SubDivision)
  local newEdgeNo = length(sd.edges) + 0x01
  push!(sd.edges, Edge(0, newEdgeNo, NoConstraint))
  push!(sd.edges, Edge(0, newEdgeNo+3, NoConstraint))
  push!(sd.edges, Edge(0, newEdgeNo+2, NoConstraint))
  push!(sd.edges, Edge(0, newEdgeNo+1, NoConstraint))
  return newEdgeNo
end

##
# basic edge index operations
"""
    base(ei::Int)

Computes base index of index `ei`.
"""
function base(ei::Int)
  return ei - ((ei-1) % 4)
end

"""
    sym(ei::Int)

Computes index of symmetric index `ei`.
"""
function sym(ei::Int)
  local qi = (ei-1) % 4
  if (qi < 2)
     return ei + 2
  else
     return ei - 2
   end
 end

"""
    rot(ei::Int)

Computes the index of the dual edge to index `ei`.

# Remarks
* The dual edge is obtained by rotating `ei` counter-clockwise by 90 degrees.
* The start(end) points of the dual edge can be interpreted as representing
  the right(left) face.
"""
function rot(ei::Int)
  local qi = (ei-1) % 4
  if (qi <3)
    return ei + 1
  else
    return ei - 3
  end
end

"""
    invrot(ei::Int)

Inverse operation of `rot(ei)`.
"""
function invrot(ei::Int)
  local qi = (ei-1) % 4
  if (qi == 0)
    return ei + 3
  else
    return ei - 1
  end
end

##
# Basic edge operations
"""
    org(sd::SubDivision, ei::Int)

Gets the index of the origin vertex of edge index `ei`.
"""
function org(sd::SubDivision, ei::Int)
  return sd.edges[ei].origin
end

"""
    dest(sd::SubDivision, ei::Int)

Gets the index of the destination vertex of edge index `ei`.
"""
function dest(sd::SubDivision, ei::Int)
  local di = sym(ei)
  return sd.edges[di].origin
end

"""
    onext(sd::SubDivision, ei::Int)

Gets the next edge starting at the origin of `ei` in counter-clockwise
direction.
"""
function onext(sd::SubDivision, ei::Int)
  return sd.edges[ei].next
end

"""
    oprev(sd::SubDivision, ei::Int)

Gets the next edge starting at the origin of `ei` in clockwise direction.
"""
function oprev(sd::SubDivision, ei::Int)
  return rot(onext(sd, rot(ei)))
end

"""
    dprev(sd::SubDivision, ei::Int)

Gets the next edge ending at the destination of `ei` in clockwise direction.
"""
function dprev(sd::SubDivision, ei::Int)
  return invrot(onext(sd, invrot(ei)))
end

"""
    lnext(sd::SubDivision, ei::Int)

Gets the next edge around the left face of `ei` and starting at destination
of `ei`.
"""
function lnext(sd::SubDivision, ei::Int)
  return rot(onext(sd, invrot(ei)))
end

"""
    lprev(sd::SubDivision, ei::Int)

Gets the previous edge around the left face of `ei` and ending at source
of `ei`.
"""
function lprev(sd::SubDivision, ei::Int)
  return sym(onext(sd, ei))
end

"""
    lface(sd::SubDivision, ei::Int)

Returns the index to face left of given edge `ei`.
"""
function lface(sd::SubDivision, ei::Int)
  return org(sd, sym(rot(ei)))
end

"""
    lface(sd::SubDivision, ei::Int)

Returns the index to face left of given edge `ei`.
"""
function rface(sd::SubDivision, ei::Int)
  return org(sd, rot(ei))
end

function lface_edges(sd::SubDivision, ei::Int)
  faceEdges = Vector{Int}()
  push!(faceEdges, ei)
  local currentEdge = lnext(sd, ei)
  while (currentEdge != ei)
    push!(faceEdges, currentEdge)
    currentEdge = lnext(sd, currentEdge)
  end

  return faceEdges
end

function lface_vertices(sd::SubDivision, ei::Int)
  local faceVertices = Vector{Int}()
  push!(faceVertices, org(sd, ei))

  local currentEdge = lnext(sd, ei)
  while (currentEdge != ei)
    push!(faceVertices, org(sd, currentEdge))
    currentEdge = lnext(sd, currentEdge)
  end

  return faceVertices
end

function ostar(sd::SubDivision, vi::Int)
  local ostarEdges = Vector{Int}()
  local startEdge = sd.vertexCache[vi]
  push!(ostarEdges, startEdge)
  local curEdge = onext(sd, startEdge)
  while (curEdge != startEdge)
    push!(ostarEdges, curEdge)
    curEdge = onext(sd, curEdge)
  end

  return ostarEdges
end

function isconstraint(sd::SubDivision, ei::Int)
  return !(sd.edges[ei].constraintType == NoConstraint)
end

function setconstraint!(sd::SubDivision, ei::Int, ct::ConstraintType)
  sd.edges[ei].constraintType = ct
  local sct = ct == OnRightSide ? OnLeftSide : OnRightSide
  sd.edges[sym(ei)].constraintType = sct
end

"""
    endpoints!(sd::SubDivision, ei::Int, vio::Int, vid::Int, lfi::Int, rfi::Int)

Assigns the start(end) point `vio`(`vid`) and the left(right) faces `lfi`(`rfi`)
to given edge `ei`.

# Remarks
* The start(end) points of the symmetric and symmetric dual edges are assigned
  correctly.
"""
function endpoints!(sd::SubDivision, ei::Int, vio::Int, vid::Int, lfi::Int, rfi::Int)
  local sei = sym(ei)

  # update end points
  sd.edges[ei].origin = vio
  sd.edges[sei].origin = vid

  # update end points of dual edges
  local re = rot(ei)
  local sre = sym(re)
  sd.edges[re].origin = rfi
  sd.edges[sre].origin = lfi

  # update vertex cache
  sd.vertexCache[vio] = ei
  sd.vertexCache[vid] = sei
end

"""
    splice!(sd::SubDivision, ai::Int, bi::Int)

Combines or splices two given edge rings `ai` and `bi`.

# Remarks
* If the origins of `ai` and `bi` are identical, it cuts the two edge rings
* If the origins differ, it joins the two edge rings
* After the operation, the new edge rings satisfy the identities
  `onext(ai) = onext(bi)` and `onext(bi) = onext(ai)` and the corresponding
  identities for the dual rings.
* The origins of the edges are neither evaluated nor are theyaffected by this
  operation and need to be updated manually.
"""
function splice!(sd::SubDivision, ai::Int, bi::Int)
  local ain = sd.edges[ai].next
  local bin = sd.edges[bi].next
  local alphai = rot(ain)
  local betai  = rot(bin)
  local alphain = sd.edges[alphai].next
  local betain  = sd.edges[betai].next

  sd.edges[ai].next = bin
  sd.edges[bi].next = ain
  sd.edges[alphai].next = betain
  sd.edges[betai].next = alphain
end

function create_subdivision{T}(v1::T, v2::T, v3::T, innerFace::T, outerFace::T)
  local sd = SubDivision(Array{Edge}(0), [-1,-1,-1], [v1, v2, v3],[innerFace, outerFace])

  # assign edges
  local ea = makeedge!(sd)
  endpoints!(sd, ea, 1, 2, 1, 2)

  local eb = makeedge!(sd)
  endpoints!(sd, eb, 2, 3, 1, 2)
  splice!(sd, sym(ea), eb)

  local ec = makeedge!(sd)
  endpoints!(sd, ec, 3, 1, 1, 2)
  splice!(sd, sym(eb), ec)
  splice!(sd, sym(ec), ea)

  return sd
end

###
# High level functions
"""
    connect!{T}(sd::SubDivision, ai::Int, bi::Int, newFace::T)

Connects the end vertex of edge `ai` to the start vertex `bi`.

# Remarks
* This operation adds another face to the subdivision. The new face will be
  to the right of the new edge and will be correctly assigned to all adjacent
  edges.
"""
function connect!{T}(sd::SubDivision{T}, ai::Int, bi::Int, newFace::T)
  local newEdge = makeedge!(sd)
  local newOrigin = dest(sd, ai)
  local newdestination = org(sd, bi)
  local oldFaceIndex = dest(sd, rot(ai))
  local newFaceIndex = length(sd.faces) + 1
  endpoints!(sd, newEdge, newOrigin, newdestination, oldFaceIndex, newFaceIndex)

  splice!(sd, newEdge, lnext(sd, ai))
  splice!(sd, sym(newEdge), bi)

  # update vertex cache
  sd.vertexCache[newOrigin] = newEdge
  sd.vertexCache[newdestination] = sym(newEdge)

  # add new face
  push!(sd.faces, newFace)

  local symNewEdge = sym(newEdge)
  local newFaceEdgeIt = lnext(sd, symNewEdge)
  while (newFaceEdgeIt != symNewEdge)
    sd.edges[sym(rot(newFaceEdgeIt))].origin = newFaceIndex
    newFaceEdgeIt = lnext(sd, newFaceEdgeIt)
  end

  return newEdge
end

function deleteedge!(sd::SubDivision, ei::Int)
  local op = oprev(sd, ei)
  local es = sym(ei)
  local esp = oprev(sd, es)

  local ov = org(sd, ei)
  local dv = org(sd, es)

  # remember faces
  local keepFaceIndex = lface(sd, ei)
  local deleteFaceIndex = rface(sd, ei)

  splice!(sd, ei, op)
  splice!(sd, es, esp)

  # mark edge as invalid
  local baseIndex = base(ei)
  sd.edges[baseIndex] = Edge(-1, -1, NoConstraint)
  sd.edges[baseIndex + 1] = Edge(-1, -1, NoConstraint)
  sd.edges[baseIndex + 2] = Edge(-1, -1, NoConstraint)
  sd.edges[baseIndex + 3] = Edge(-1, -1, NoConstraint)

  # update vertex cache
  sd.vertexCache[ov] = op
  sd.vertexCache[dv] = esp

  # update faces
  local keepFaceEdges = lface_edges(sd, op)
  for curEdge in keepFaceEdges
    sd.edges[sym(rot(curEdge))].origin = keepFaceIndex
  end
end

"""
    removedeletededges!(sd::SubDivision)

Removes all edges that were marked as deleted in DelaunayMeshes.
"""
function removedeletededges!(sd::SubDivision)

  # compute partial sum of offsets
  local delEdgeIndex = find((x::Edge) -> x.origin == -1, sd.edges[1:4:end])
  local offsets = zeros(Int, length(sd.edges[1:4:end]))
  offsets[delEdgeIndex] = 4
  cumsum!(offsets, offsets)

  # substract offset of edges
  for (i, e) in enumerate(sd.edges)
    local baseOffset = Int((base(e.next) - 1)/4 + 1)
    e.next -= offsets[baseOffset]
  end

  # remove all offending edges
  filter!((x::Edge) -> x.origin != -1, sd.edges)

  # rebuild vertex cache
  fill!(sd.vertexCache, -1)
  for i in 1:2:length(sd.edges)
    local org = sd.edges[i].origin
    if (sd.vertexCache[org] == -1)
      sd.vertexCache[org] = i
    end
  end
end


text:
struct S
 a end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Documenter, CSV

makedocs(modules = [CSV], sitename = "CSV.jl", format = Documenter.HTML(edit_link = "main", canonical = "https://juliadata.github.io/CSV.jl/stable/"))

deploydocs(repo = "github.com/JuliaData/CSV.jl.git", devbranch = "main")


text:
using Gaston, TimeSeries

module PriceCharts

using Gaston, TimeSeries

export plotreturn,
       plotseries,
## testing
       @pricecharts

################## include files #####################

include("plot.jl")
include("testpricecharts.jl")

end #module


text:
#
# Defines AstrodynamicalCalculations methods for
# OrbitalStates types
#

export massparameter, massparameters, normalized_massparameter
export primary_massparameter, secondary_massparameter
export RAAN, argument_of_periapsis, inclination
export periapsis_velocity, apoapsis_velocity
export mean_motion_vector
export distance_to_primary, distance_to_secondary
export massparameter, massparameters
export primary_massparameter, secondary_massparameter

include("States.jl")
include("Systems.jl")


text:
using YearMonths
using CVMData
using Test

ym = YearMonths.YearMonth(2020,1)
quotes = CVMData.Funds.get_month_quotes(ym)
@test length(quotes) > 0

ym = YearMonths.YearMonth(2020,4)
quotes = CVMData.Funds.get_month_quotes(ym)
@test length(quotes) > 0

ym = YearMonths.YearMonth(2021,1)
quotes = CVMData.Funds.get_month_quotes(ym)
@test length(quotes) > 0

ym = YearMonths.YearMonth(2999,1)
@test_throws BoundsError CVMData.Funds.get_month_quotes(ym)


text:
@testset "Catchability Deviation scaling" begin
    run_sims(QDevScalingSpec, 1;
             prep_file = "prep.h5",
             csv = true,
             feather = true)

    # Construct all of the expected output file names
    output_files = joinpath.("qdevscaling",
                             "repl_01",
                             "qdevscale_" .*
                             string.(1:length(sim_values(QDevScalingSpec)),
                                     pad = 2) .*
                             ["_popstate.h5" "_popstate.csv" "_catch.csv" "_popstate.feather" "_catch.feather"])

    @test isdir("qdevscaling")
    @test isdir("qdevscaling", "repl_01")
    @test all(isfile.(output_files))
end


text:
@testset "Subsets" begin

  # need expanding
  subset_result = Subset("Subset 1")
  @test subset_result.Name == "Subset 1"
  @test name(subset_result) == "Subset 1"

  subset_2 = JSON.parse("""
  {
    "Name": "Subset 2"
  }
  """)

  @test name(Subset(subset_2)) == "Subset 2"
  @test Subset(subset_2) != subset_result

end


text:
using Combinatorics
using GraphModularDecomposition
using GraphModularDecomposition: overlap
using IterTools
using LinearAlgebra
using LinearAlgebra: checksquare
using Random
using SparseArrays
using Test

const \ = setdiff

function is_module(G::AbstractMatrix, S::Vector{Int})
    n = checksquare(G)
    all(G[i,k] == G[j,k] && G[k,i] == G[k,j]
        for i in S for j in S for k in 1:n if k ∉ S)
end

function all_modules(G::AbstractMatrix)
    n = checksquare(G)
    filter!(collect(subsets(1:n))) do S
        1 < length(S) && is_module(G, S)
    end
end

function maximal_independent_sets(G::AbstractMatrix)
    n = checksquare(G)
    C(i, j) = G[i,j] == G[j,i] == 0
    filter!(collect(subsets(1:n))) do S
        all(C(i, j) for i in S for j in S) &&
        !any(all(C(i, j) for i in S) for j in 1:n if j ∉ S)
    end
end

function all_strong_modules(G::AbstractMatrix)
    modules = all_modules(G)
    filter(modules) do A
        all(modules) do B
            !overlap(A, B)
        end
    end
end

function is_modular_permutation(G::AbstractMatrix, p::Vector{Int};
    modules = all_strong_modules(G))
    isempty(modules) && return true
    diffs = map(modules) do M
        diff(findall(in(M), p))
    end
    maximum(maximum, diffs) == 1
end

findin_partition(P, S) = sort!(map(x->findfirst(X->x in X, P), S))

function is_modular_partition(G::AbstractMatrix, P::Vector{Vector{Int}};
    modules = all_strong_modules(G))
    sort!(vcat(P...)) == collect(1:size(G,2)) || error("not a partition")
    diffs = map(modules) do M
        diff(findin_partition(P, M))
    end
    maximum(maximum, diffs) == 1
end

is_tournament(G::AbstractMatrix) = all(G + G' + I .== 1)


text:
abstract type AbstractAllocationRecord end
struct AllocationRecord <: AbstractAllocationRecord
    allocations::Vector{Array}
    initial_sizes::Vector{Any}  # collection of differently sizes tuples of integers
end
AllocationRecord() = AllocationRecord(Vector{Array}(), Vector{Any}())

struct FrozenAllocationRecord{A,S} <: AbstractAllocationRecord
    allocations::A
    initial_sizes::S
end

function FrozenAllocationRecord(record)
    return FrozenAllocationRecord(Tuple(record.allocations), Tuple(record.initial_sizes))
end

function Base.copy(record::AllocationRecord)
    return AllocationRecord(copy.(record.allocations), record.initial_sizes)
end

function Base.copy(record::FrozenAllocationRecord)
    return FrozenAllocationRecord(copy.(record.allocations), record.initial_sizes)
end

"""
    freeze(record)

Freezes an allocation `record`, to not allow new allocations to be added.
This converts the memory of the what happenned to be stored as a `Tuple`.
This could give better performance during replay; or it could give worse.
"""
freeze(record::AbstractAllocationRecord) = FrozenAllocationRecord(record)

"""
    reinitialize!(record)

Reset all arrays in the allocation `record` back to their initial sizes.
"""
function reinitialize!(record)
    for ii in eachindex(record.allocations)
        alloc = record.allocations[ii]
        sz = record.initial_sizes[ii]
        if size(alloc) !== sz
            # fix any vectors that were e.g. `push!`ed to.
            resize!(alloc, sz...)
        end
    end
    return record
end

# This can be use to write the record to file so that one can reuse it.
function Base.show(io::IO, record::AllocationRecord)
    allocs = join(("$(typeof(x))(undef, $(size(x)))" for x in record.allocations), ", ")
    sizes = join(record.initial_sizes, ", ")
    println(io, "$(@__MODULE__).AllocationRecord(")
    println(io, "    [$allocs],")
    println(io, "    [$sizes]")
    print(io, ")")
end

Base.hash(record::AbstractAllocationRecord, k::UInt) = hash(record.initial_sizes)
function Base.:(==)(r1::AbstractAllocationRecord, r2::AbstractAllocationRecord)
    return (
        r1.initial_sizes == r2.initial_sizes
        && eltype.(r1.allocations) == eltype.(r2.allocations)
    )
end


text:
using TestPackage
using Test

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


text:
export NBCopulaVCObs, NBCopulaVCModel
### first structures
mutable struct NBCopulaVCObs{T <: BlasReal, D, Link} # d changes, so must be mutable
    # data
    y::Vector{T}
    X::Matrix{T}
    V::Vector{Matrix{T}}
    n::Int
    p::Int
    m::Int # number of variance components
    # working arrays
    ∇β::Vector{T}   # gradient wrt β
    ∇resβ::Matrix{T}# residual gradient matrix d/dβ_p res_ij (each observation has a gradient of residual is px1)
    ∇τ::Vector{T}   # gradient wrt τ
    ∇θ::Vector{T}   # gradient wrt θ2
    ∇r::Vector{T}   # gradient wrt r (NB)
    Hβ::Matrix{T}   # Hessian wrt β
    Hθ::Matrix{T}   # Hessian wrt variance components θ
    Hr::Matrix{T}   # Hessian wrt r (NB)
    Hτ::Matrix{T}   # Hessian wrt τ
    res::Vector{T}  # residual vector res_i
    t::Vector{T}    # t[k] = tr(V_i[k]) / 2
    q::Vector{T}    # q[k] = res_i' * V_i[k] * res_i / 2
    xtx::Matrix{T}  # Xi'Xi
    storage_n::Vector{T}
    storage_n2::Vector{T}
    storage_n3::Vector{T}
    storage_p1::Vector{T}
    storage_p2::Vector{T}
    storage_nn::Matrix{T}
    storage_np::Matrix{T}
    storage_pp::Matrix{T}
    added_term_numerator::Matrix{T}
    added_term2::Matrix{T}
    η::Vector{T}    # η = Xβ systematic component
    μ::Vector{T}    # μ(β) = ginv(Xβ) # inverse link of the systematic component
    varμ::Vector{T} # v(μ_i) # variance as a function of the mean
    dμ::Vector{T}   # derivative of μ
    d::D            # distribution()
    link::Link      # link function ()
    wt::Vector{T}   # weights wt for GLM.jl
    w1::Vector{T}   # working weights in the gradient = dμ/v(μ)
    w2::Vector{T}   # working weights in the information matrix = dμ^2/v(μ)
end

function NBCopulaVCObs(
    y::Vector{T},
    X::Matrix{T},
    V::Vector{Matrix{T}},
    d::D,
    link::Link) where {T <: BlasReal, D, Link}
    n, p, m = size(X, 1), size(X, 2), length(V)
    @assert length(y) == n "length(y) should be equal to size(X, 1)"
    # working arrays
    ∇β  = Vector{T}(undef, p)
    ∇resβ  = Matrix{T}(undef, n, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇θ  = Vector{T}(undef, m)
    ∇r  = Vector{T}(undef, 1)
    Hβ  = Matrix{T}(undef, p, p)
    Hθ  = Matrix{T}(undef, m, m)
    Hr  = Matrix{T}(undef, 1, 1)
    Hτ  = Matrix{T}(undef, 1, 1)
    res = Vector{T}(undef, n)
    t   = [tr(V[k])/2 for k in 1:m]
    q   = Vector{T}(undef, m)
    xtx = transpose(X) * X
    storage_n = Vector{T}(undef, n)
    storage_n2 = Vector{T}(undef, n)
    storage_n3 = Vector{T}(undef, n)
    storage_p1 = Vector{T}(undef, p)
    storage_p2 = Vector{T}(undef, p)
    storage_nn = Matrix{T}(undef, n, n) # stores Γ = a1*V1 + ... am*Vm
    storage_np = Matrix{T}(undef, n, p)
    storage_pp = Matrix{T}(undef, p, p)
    added_term_numerator = Matrix{T}(undef, n, p)
    added_term2 = Matrix{T}(undef, p, p)
    η = Vector{T}(undef, n)
    μ = Vector{T}(undef, n)
    varμ = Vector{T}(undef, n)
    dμ = Vector{T}(undef, n)
    wt = Vector{T}(undef, n)
    fill!(wt, one(T))
    w1 = Vector{T}(undef, n)
    w2 = Vector{T}(undef, n)
    # constructor
    NBCopulaVCObs{T, D, Link}(y, X, V, n, p, m, ∇β, ∇resβ, ∇τ, ∇θ, ∇r, Hβ, Hθ, Hr,
        Hτ, res, t, q, xtx, storage_n, storage_n2, storage_n3, storage_p1,
        storage_p2, storage_nn, storage_np, storage_pp, added_term_numerator, added_term2,
        η, μ, varμ, dμ, d, link, wt, w1, w2)
end

"""
NBCopulaVCModel
NBCopulaVCModel(gcs)
Negative Binomial copula variance component model, which contains a vector of
`NBCopulaVCObs` as data, model parameters, and working arrays.
"""
struct NBCopulaVCModel{T <: BlasReal, D, Link} <: MathProgBase.AbstractNLPEvaluator
    # data
    data::Vector{NBCopulaVCObs{T, D, Link}}
    Ytotal::T
    ntotal::Int     # total number of singleton observations
    n::Int          # number of random vectors
    p::Int          # number of mean parameters in linear regression
    m::Int          # number of variance components
    # parameters
    β::Vector{T}    # p-vector of mean regression coefficients
    τ::Vector{T}    # inverse of linear regression variance parameter
    θ::Vector{T}    # m-vector: [θ12, ..., θm2]
    r::Vector{T}    # r parameter for negative binomial
    # working arrays
    ∇β::Vector{T}   # gradient from all observations
    ∇τ::Vector{T}
    ∇θ::Vector{T}
    ∇r::Vector{T}
    XtX::Matrix{T}  # X'X = sum_i Xi'Xi
    Hβ::Matrix{T}    # Hessian from all observations
    Hθ::Matrix{T}
    Hr::Matrix{T}
    Hτ::Matrix{T}
    Ainv::Matrix{T}
    Aevec::Matrix{T}
    M::Matrix{T}
    vcov::Matrix{T}
    ψ::Vector{T}
    TR::Matrix{T}         # n-by-m matrix with tik = tr(Vi[k]) / 2
    QF::Matrix{T}         # n-by-m matrix with qik = res_i' Vi[k] res_i
    hess1::Matrix{T}      # holds transpose(gcm.QF) * Diagonal(gcm.storage_n) required for outer product in hessian term 1
    hess2::Matrix{T}      # holds transpose(gcm.TR) * Diagonal(gcm.storage_n2) required for outer product in hessian term 2
    storage_n::Vector{T}
    storage_n2::Vector{T}
    storage_m::Vector{T}
    storage_θ::Vector{T}
    d::Vector{D}
    link::Vector{Link}
    penalized::Bool
end

function NBCopulaVCModel(gcs::Vector{NBCopulaVCObs{T, D, Link}}; penalized::Bool = false) where {T <: BlasReal, D, Link}
    n, p, m = length(gcs), size(gcs[1].X, 2), length(gcs[1].V)
    β   = Vector{T}(undef, p)
    τ   = [1.0]
    θ   = Vector{T}(undef, m)
    r   = [1.0]
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇θ  = Vector{T}(undef, m)
    ∇r  = Vector{T}(undef, 1)
    XtX = zeros(T, p, p) # sum_i xi'xi
    Hβ  = Matrix{T}(undef, p, p)
    Hθ  = Matrix{T}(undef, m, m)
    Hr  = Matrix{T}(undef, 1, 1)
    Hτ  = Matrix{T}(undef, 1, 1)
    Ainv    = zeros(T, p + m + 1, p + m + 1)
    Aevec   = zeros(T, p + m + 1, p + m + 1)
    M       = zeros(T, p + m + 1, p + m + 1)
    vcov    = zeros(T, p + m + 1, p + m + 1)
    ψ       = Vector{T}(undef, p + m + 1)
    TR  = Matrix{T}(undef, n, m) # collect trace terms
    Ytotal = 0.0
    ntotal = 0.0
    d = Vector{D}(undef, n)
    link = Vector{Link}(undef, n)
    for i in eachindex(gcs)
        ntotal  += length(gcs[i].y)
        Ytotal  += sum(gcs[i].y)
        BLAS.axpy!(one(T), gcs[i].xtx, XtX)
        TR[i, :] = gcs[i].t
        d[i] = gcs[i].d
        link[i] = gcs[i].link
    end
    QF        = Matrix{T}(undef, n, m)
    hess1     = Matrix{T}(undef, m, n)
    hess2     = Matrix{T}(undef, m, n)
    storage_n = Vector{T}(undef, n)
    storage_n2 = Vector{T}(undef, n)
    storage_m = Vector{T}(undef, m)
    storage_θ = Vector{T}(undef, m)
    NBCopulaVCModel{T, D, Link}(gcs, Ytotal, ntotal, n, p, m, β, τ, θ, r,
        ∇β, ∇τ, ∇θ, ∇r, XtX, Hβ, Hθ, Hr, Hτ, Ainv, Aevec, M, vcov, ψ, TR, QF, hess1, hess2,
        storage_n, storage_n2, storage_m, storage_θ, d, link, penalized)
end

"""
    loglikelihood!(gc::NBCopulaVCObs{T, D, Link}, β, τ, θ, r)
Calculates the loglikelihood of observing `y` given parameters for `β`, `τ = 1`, `θ`, and `r`for Negative Binomial distribution.
"""
function loglikelihood!(
  gc::NBCopulaVCObs{T, D, Link},
  β::Vector{T},
  τ::T, # inverse of linear regression variance
  θ::Vector{T},
  r::T,
  needgrad::Bool = false,
  needhess::Bool = false;
  penalized::Bool = false
  ) where {T <: BlasReal, D, Link}
  # n, p, m = size(gc.X, 1), size(gc.X, 2), length(gc.V)
  needgrad = needgrad || needhess
  if needgrad
      fill!(gc.∇β, 0)
      fill!(gc.∇τ, 0)
      fill!(gc.∇θ, 0)
  end
  needhess && fill!(gc.Hβ, 0)
  fill!(gc.∇β, 0.0)
  update_res!(gc, β)
  standardize_res!(gc)
  fill!(gc.∇resβ, 0.0) # fill gradient of residual vector with 0
  std_res_differential!(gc) # this will compute ∇resβ

  # evaluate copula loglikelihood
  @inbounds for k in 1:gc.m
      mul!(gc.storage_n, gc.V[k], gc.res) # storage_n = V[k] * res
      if needgrad
          BLAS.gemv!('T', θ[k], gc.∇resβ, gc.storage_n, 1.0, gc.∇β) # stores ∇resβ*Γ*res (standardized residual) SHOULDNT THIS BE res^t*Γ*res? NOT ∇resβ*Γ*res
      end
      gc.q[k] = dot(gc.res, gc.storage_n) / 2
  end
  # 2nd term of loglikelihood
  logl = QuasiCopula.component_loglikelihood(gc, r)
  # 1st term of loglikelihood
  tsum = dot(θ, gc.t)
  logl += -log(1 + tsum)
  # 3rd term of loglikelihood
  qsum  = dot(θ, gc.q)
  logl += log(1 + qsum)

  # add L2 ridge penalty
  if penalized
      logl -= 0.5 * dot(θ, θ)
  end
  if needgrad
      inv1pq = inv(1 + qsum)

        # gradient with respect to r
        gc.∇r .= nb_first_derivative(gc, θ, r)
        if penalized
            gc.∇θ .-= θ
        end
      if needhess
          BLAS.syrk!('L', 'N', -abs2(inv1pq), gc.∇β, 1.0, gc.Hβ) # only lower triangular
          # does adding this term to the approximation of the hessian violate negative semidefinite properties?
          fill!(gc.added_term_numerator, 0.0) # fill gradient with 0
          fill!(gc.added_term2, 0.0) # fill hessian with 0
          @inbounds for k in 1:gc.m
              mul!(gc.added_term_numerator, gc.V[k], gc.∇resβ) # storage_n = V[k] * res
              BLAS.gemm!('T', 'N', θ[k], gc.∇resβ, gc.added_term_numerator, one(T), gc.added_term2)
          end
          gc.added_term2 .*= inv1pq
          gc.Hβ .+= gc.added_term2
          gc.Hβ .+= QuasiCopula.glm_hessian(gc)

          # hessian for r
          gc.Hr .= nb_second_derivative(gc, θ, r)
      end
      gc.storage_p2 .= gc.∇β .* inv1pq
      gc.res .= gc.y .- gc.μ
      gc.∇β .= QuasiCopula.glm_gradient(gc)
      gc.∇β .+= gc.storage_p2
  end
  logl
end

function loglikelihood!(
    gcm::NBCopulaVCModel{T, D, Link},
    needgrad::Bool = false,
    needhess::Bool = false
    ) where {T <: BlasReal, D, Link}
    logl = zero(T)
    if needgrad
        fill!(gcm.∇β, 0.0)
        fill!(gcm.∇θ, 0.0)
        fill!(gcm.∇r, 0.0)
    end
    if needhess
        fill!(gcm.Hβ, 0.0)
        fill!(gcm.Hθ, 0.0)
        fill!(gcm.Hr, 0.0)
    end
    logl = zeros(Threads.nthreads())
    Threads.@threads for i in eachindex(gcm.data)
        @inbounds logl[Threads.threadid()] += loglikelihood!(gcm.data[i], gcm.β,
            gcm.τ[1], gcm.θ, gcm.r[1], needgrad, needhess; penalized = gcm.penalized)
    end
    @inbounds for i in eachindex(gcm.data)
        if needgrad
            gcm.∇β .+= gcm.data[i].∇β
            gcm.∇r .+= gcm.data[i].∇r
        end
        if needhess
            gcm.Hβ .+= gcm.data[i].Hβ
            gcm.Hr .+= gcm.data[i].Hr
        end
    end
    if needgrad
        gcm.∇θ .= update_∇θ!(gcm)
    end
    if needhess
        gcm.Hθ .= update_Hθ!(gcm)
    end
    # logl
    return sum(logl)
end

"""
1st derivative of loglikelihood of a sample with θ being the variance components
"""
function nb_first_derivative(gc::NBCopulaVCObs, θ::Vector{T}, r::Number) where T <: BlasReal
    s = zero(T)
    # 2nd term of logl
    @inbounds for j in eachindex(gc.y)
        yi, μi = gc.y[j], gc.μ[j]
        s += -(yi+r)/(μi+r) - log(μi+r) + 1 + log(r) + digamma(r+yi) - digamma(r)
    end
    # 3rd term of logl
    # Γ = θ' * gc.V # Γ = a1*V1 + ... + am*Vm
    # η = gc.η
    # D = Diagonal([sqrt(exp(η[j])*(exp(η[j])+r) / r) for j in 1:length(η)])
    # dD = Diagonal([-exp(2η[i]) / (2r^1.5 * sqrt(exp(η[i])*(exp(η[i])+r))) for i in 1:length(η)])
    # dresid = -inv(D)*dD*resid
    # s += resid'*Γ*dresid / (1 + 0.5resid'*Γ*resid)
    resid = gc.res # res = inv(D)(y - μ)
    Γ = gc.storage_nn # Γ = a1*V1 + ... + am*Vm
    η = gc.η
    for j in 1:length(η)
        gc.storage_n[j] = inv(sqrt(exp(η[j])*(exp(η[j])+r) / r)) # storage_n = inv(Di) = 1 / sqrt(var(yi))
    end
    for j in 1:length(η)
        # storage_n2 = -inv(D) * d2D
        gc.storage_n2[j] = -gc.storage_n[j] *
            ((exp(3η[j]) / (4r^1.5 * (exp(η[j])*(exp(η[j])+r))^(1.5))) +
            (3exp(2η[j]) / (4r^(2.5)*sqrt(exp(η[j])*(exp(η[j])+r)))))
    end
    for j in 1:length(η)
        # storage_n = inv(D) * dD
        gc.storage_n[j] *= -exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r)))
    end
    for j in 1:length(η)
        # storage_n2 = 2inv(D)*dD*inv(D)*dD -inv(D)*d2D
        gc.storage_n2[j] += 2 * abs2(gc.storage_n[j])
    end
    gc.storage_n .*= -1.0 .* resid # storage_n = dr(β) = derivative of residuals
    gc.storage_n2 .*= resid # storage_n2 = dr²(β) = 2nd derivative of residuals
    mul!(gc.storage_n3, Γ, resid) # storage_n3 = Γ * resid
    denom = 1 + 0.5 * dot(resid, gc.storage_n3)
    mul!(gc.storage_n3, Γ, gc.storage_n) # storage_n3 = Γ * dresid
    term1 = (dot(resid, gc.storage_n3) / denom)^2 # (resid'*Γ*dresid / denom)^2
    term2 = dot(gc.storage_n, gc.storage_n3) / denom # term2 = dresid'*Γ*dresid / denom
    mul!(gc.storage_n3, Γ, gc.storage_n2) # storage_n3 = Γ * d2resid
    term3 = dot(resid, gc.storage_n3) / denom # term3 = resid'*Γ*d2resid / denom
    s += -term1 + term2 + term3
    return s
end

"""
2nd derivative of loglikelihood of a sample with θ being the variance components
"""
function nb_second_derivative(gc::NBCopulaVCObs, θ::Vector{T}, r::Number) where T <: BlasReal
    s = zero(T)
    # 2nd term of logl
    @inbounds for j in eachindex(gc.y)
        yi, μi = gc.y[j], gc.μ[j]
        s += (yi+r)/(μi+r)^2 - 2/(μi+r) + 1/r + trigamma(r+yi) - trigamma(r)
    end
    # 3rd term of logl
    # Γ = θ' * gc.V # Γ = a1*V1 + ... + am*Vm
    # η = gc.η
    # D = Diagonal([sqrt(exp(η[j])*(exp(η[j])+r) / r) for j in 1:length(η)])
    # dD = Diagonal([-exp(2η[i]) / (2r^1.5 * sqrt(exp(η[i])*(exp(η[i])+r))) for i in 1:length(η)])
    # d2D = Diagonal([(exp(3η[i]) / (4r^1.5 * (exp(η[i])*(exp(η[i])+r))^(1.5))) +
    #     (3exp(2η[i]) / (4r^(2.5)*sqrt(exp(η[i])*(exp(η[i])+r)))) for i in 1:length(η)])
    # resid = gc.res
    # dresid = -inv(D)*dD*resid
    # d2resid = (2inv(D)*dD*inv(D)*dD - inv(D)*d2D)*resid
    # denom = 1 + 0.5resid'*Γ*resid
    # term1 = (resid'*Γ*dresid / denom)^2
    # term2 = dresid'*Γ*dresid / denom
    # term3 = resid'*Γ*d2resid / denom
    # s += -term1 + term2 + term3
    resid = gc.res # res = inv(D)(y - μ)
    Γ = gc.storage_nn # Γ = a1*V1 + ... + am*Vm
    η = gc.η
    for j in 1:length(η)
        gc.storage_n[j] = inv(sqrt(exp(η[j])*(exp(η[j])+r) / r)) # storage_n = inv(Di) = 1 / sqrt(var(yi))
    end
    for j in 1:length(η)
        # storage_n2 = -inv(D) * d2D
        gc.storage_n2[j] = -gc.storage_n[j] *
            ((exp(3η[j]) / (4r^1.5 * (exp(η[j])*(exp(η[j])+r))^(1.5))) +
            (3exp(2η[j]) / (4r^(2.5)*sqrt(exp(η[j])*(exp(η[j])+r)))))
    end
    for j in 1:length(η)
        # storage_n = inv(D) * dD
        gc.storage_n[j] *= -exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r)))
    end
    for j in 1:length(η)
        # storage_n2 = 2inv(D)*dD*inv(D)*dD -inv(D)*d2D
        gc.storage_n2[j] += 2 * abs2(gc.storage_n[j])
    end
    gc.storage_n .*= -1.0 .* resid # storage_n = dr(β) = derivative of residuals
    gc.storage_n2 .*= resid # storage_n2 = dr²(β) = 2nd derivative of residuals
    mul!(gc.storage_n3, Γ, resid) # storage_n3 = Γ * resid
    denom = 1 + 0.5 * dot(resid, gc.storage_n3)
    mul!(gc.storage_n3, Γ, gc.storage_n) # storage_n3 = Γ * dresid
    term1 = (dot(resid, gc.storage_n3) / denom)^2 # (resid'*Γ*dresid / denom)^2
    term2 = dot(gc.storage_n, gc.storage_n3) / denom # term2 = dresid'*Γ*dresid / denom
    mul!(gc.storage_n3, Γ, gc.storage_n2) # storage_n3 = Γ * d2resid
    term3 = dot(resid, gc.storage_n3) / denom # term3 = resid'*Γ*d2resid / denom
    s += -term1 + term2 + term3

    return s
end


text:
_tokenize(text::String) = map(lowercase, split(text))
_all() = reduce(union, values(_INDEX));


text:
# Measures of anisotropy

"""
    Au(C) -> au

Return the Universal Elastic Anisotropy Index, `au`, of the tensor `C`.

See: Ranganathan & Ostoja-Starzewksi, Universal elastic anisotropy index,
     Phys Rev Lett (2008) vol. 101 (5) pp. 055504
"""
function Au(C)
    Kv, Gv, Kr, Gr = VoigtK(C), VoigtG(C), ReussK(C), ReussG(C)
    5*(Gv/Gr) + Kv/Kr - 6
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE
using Behavior

using TypeDBClient

# include(joinpath(@__DIR__,"test/behaviour/config/ConfigEnvironment.jl"))

g = TypeDBClient

@given("connection has been opened") do context
    client = g.CoreClient("127.0.0.1",1729)
    context[:client] = client
end

@given("connection does not have any database") do context
    all_dbs = g.get_all_databases(context[:client])
    if !isempty(all_dbs)
        delete_all_databases(context[:client])
        all_dbs = g.get_all_databases(context[:client])
    end
    @expect length(all_dbs) == 0
end

sessions(context) = collect(values(context[:client].sessions))

transactions(context::Behavior.StepDefinitionContext) = ordered_transaction_result(context[:session])
transactions(session::g.CoreSession) = ordered_transaction_result(session)
function ordered_transaction_result(session::g.CoreSession)
    result = g.CoreTransaction[]
    for (_, transaction) in session.transactions
        push!(result, transaction)
    end
    return result
end

trans_read = g.Proto.Transaction_Type.READ
trans_write = g.Proto.Transaction_Type.WRITE

################# Utility functions #######################

function group_count_items(items::Vector)
    unique_vector = unique(items)
    result = Dict{Any, Number}()
    for search in unique_vector
        count_nr = 0
        for item in items
            search == item && (count_nr +=1)
        end
        result[search] = count_nr
    end
    return result
end


text:
struct bootstrapResult
        fullDataEstimate
        replicateEstimate
        resampleIndexBB
        resampleDataBB
end
function bootstrapStatistic(data::Array{<:Real,1}, statisticHandle::Function, sampler::samplingAlgorithm, blockLength::Integer, NbootstrapReplicates::Integer)
        # Bootstrap a statistic
        fullDataEstimate = statisticHandle(data)
        estimateType = typeof(fullDataEstimate)

        # Default number of blocks
        Ndata = length(data)
        Nblocks = ceil(Int, Ndata/blockLength)

        # Sampling
        bootstrapSampleHandle = sampler.handle
        resampleIndexBB, resampleDataBB = bootstrapSampleHandle(data, blockLength, Nblocks, NbootstrapReplicates)

        # Statistic on samples
        replicateEstimate = Vector{estimateType}(undef, NbootstrapReplicates)
        for ii in 1:NbootstrapReplicates
                replicateEstimate[ii] = statisticHandle(resampleDataBB[ii])
        end

        bootstrapResultOut = bootstrapResult(fullDataEstimate, replicateEstimate, resampleIndexBB, resampleDataBB)

        return bootstrapResultOut
end


text:
module PSU
using Libdl # for NI dlls
include("../INSTR_TYPES.jl")
include("../lib_locations.jl")
include("../visa/VISA.jl")
## include all PSU files
include("./AgilentE36xx/AgilentE364x.jl")

end #endmodule


text:
### Source excitation generation ###

function excite(pitch::StridedVector{Cdouble}, hopsize=80;
                interp_period::Int=1,
                gaussian::Bool=false,
                seed::Int=1)
    expectedlen = trunc(hopsize * (length(pitch) - 1))
    excitation = Array{Cdouble}(expectedlen)

    ccall((:excite, libSPTK), Void,
          (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cint, Bool, Cint),
          pitch, length(pitch), excitation, hopsize, interp_period,
          gaussian, seed)

    excitation
end


text:
nprocs()==1 && addprocs()

@everywhere logis(x) = 1 ./ (1 .+ exp(-x))
@everywhere dlogis(x) = logis(x) .* (1 .- logis(x))

@everywhere function dfa!(w1,w2,w3,B1,B2,inp,out)
    a1 = w1 * inp
    z1 = logis(a1)
    a2 = w2 * z1
    z2 = logis(a2)
    ay = w3 * z2
    y = logis(ay)
    err = y - out
    d_a1 = (B1.*err) .* dlogis(a1)
    d_a2 = (B2.*err) .* dlogis(a2)
    w1 .-= d_a1 * inp'
    w2 .-= d_a2 * z1'
    w3 .-= err * z2'
    return Dict("output"=>y,"err"=>err)
end

@everywhere function bp!(bp_w1,bp_w2,bp_w3,inp,out)
    bp_a1 = bp_w1*inp
    bp_z1 = logis(bp_a1)
    bp_a2 = bp_w2 * bp_z1
    bp_z2 = logis(bp_a2)
    bp_ay = bp_w3 * bp_z2
    bp_y = logis(bp_ay)
    bp_err = bp_y - out
    bp_d_a3 = (bp_err) .* dlogis(bp_ay)
    bp_d_a2 = (bp_w3' * bp_d_a3) .* dlogis(bp_a2)
    bp_d_a1 = (bp_w2' * bp_d_a2) .* dlogis(bp_a1)
    bp_w1 .-= bp_d_a1 * inp'
    bp_w2 .-= bp_d_a2 * bp_z1'
    bp_w3 .-= bp_d_a3 * bp_z2'
    return Dict("output"=>bp_y,"err"=>bp_err)
end

inp = [1 1 ; 0 1 ; 1 0 ; 0 0]'
out = [0 1 1 0]
y=zeros(size(out))
n_hidden = 15; # Number of hidden units
num_iterations = 1000; #Number of learning steps
trials = 50;

e_store = SharedArray(Float64,(num_iterations,trials))
bp_e_store = SharedArray(Float64,(num_iterations,trials))
y_store = SharedArray(Float64,(size(out,2),num_iterations))
bp_y_store = SharedArray(Float64,(size(out,2),num_iterations))

@elapsed @sync @parallel for jj = 1:trials
    w1 = randn(n_hidden,size(inp,1))
    w2 = randn(n_hidden,n_hidden)
    w3 = randn(size(out,1),n_hidden)

    B1=rand(n_hidden,1)
    B2=rand(n_hidden,1)

    bp_w1 = randn(n_hidden,size(inp,1))
    bp_w2 = randn(n_hidden,n_hidden)
    bp_w3 = randn(size(out,1),n_hidden)

    for ii = 1:num_iterations
        dfa_res = dfa!(w1,w2,w3,B1,B2,inp,out)
        bp_res = bp!(bp_w1,bp_w2,bp_w3,inp,out)
        if jj == 1
            y_store[:,ii] = dfa_res["output"]'
            bp_y_store[:,ii] = bp_res["output"]'
        end
        e_store[ii,jj] = sum(abs(dfa_res["err"]))
        bp_e_store[ii,jj] = sum(abs(bp_res["err"]))
    end
end

using Plots
l = @layout [
    a b
    c d
]

gr()
p1=plot(y_store',label=["1,1" "0,1" "1,0" "0,0"],title="DFA")
p2=plot(bp_y_store',label=["1,1" "0,1" "1,0" "0,0"],title="BP")
p3=plot([e_store[:,1] bp_e_store[:,1]],label=["DFA" "BP"],title="First Trial")
p4=plot([mean(e_store,2) mean(bp_e_store,2)],label=["DFA" "BP"],title="Mean Error")
plot(p1,p2,p3,p4,layout=l)


text:
using Test
using Peridynamics.Discretization

@testset "tensor product grid" begin
    @test create_grid([0; 2], [2]) == [0.0, 2.0]
    @test create_grid([0; 2], [2]) == [0.0, 2.0]
    @test create_grid([0; 2], [3]) == [0.0, 1.0, 2.0]
    @test create_grid([0; 2], [2.0]) == [0.0, 2.0]
    @test create_grid([0; 2], [1.0]) == [0.0, 1.0, 2.0]
    @test create_grid([0; 2]) == [0.0, 1.0, 2.0]
    @test create_grid([0 0; 2 2], [2 2]) == [0.0 0.0; 0.0 2.0; 2.0 0.0; 2.0 2.0]
    @test create_grid([0 0; 1 1]) == [0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0]
    @test create_grid([0 0; 2 2], [2 3]) == [0.0 0.0; 0.0 1.0; 0.0 2.0; 2.0 0.0; 2.0 1.0; 2.0 2.0]
    @test create_grid([0 0; 2 2], [2.0 2.0]) == [0.0 0.0; 0.0 2.0; 2.0 0.0; 2.0 2.0]
    @test create_grid([0 0; 2 2], [2.0 1.0]) == [0.0 0.0; 0.0 1.0; 0.0 2.0; 2.0 0.0; 2.0 1.0; 2.0 2.0]
    @test create_grid([0 0 0; 2 2 2], [2 2 2]) == [0.0 0.0 0.0; 0.0 0.0 2.0; 0.0 2.0 0.0; 0.0 2.0 2.0; 2.0 0.0 0.0; 2.0 0.0 2.0; 2.0 2.0 0.0; 2.0 2.0 2.0]
    @test create_grid([0 0 0; 1 1 1]) == [0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 0.0; 0.0 1.0 1.0; 1.0 0.0 0.0; 1.0 0.0 1.0; 1.0 1.0 0.0; 1.0 1.0 1.0]
    create_grid([0 0 0; 2 2 2], [2 2 2])
    @test create_grid([0 0 0; 2 2 2], [2 2 3]) == [0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 0.0 2.0; 0.0 2.0 0.0; 0.0 2.0 1.0; 0.0 2.0 2.0; 2.0 0.0 0.0; 2.0 0.0 1.0; 2.0 0.0 2.0; 2.0 2.0 0.0; 2.0 2.0 1.0; 2.0 2.0 2.0]
    @test create_grid([0 0 0; 2 2 2], [2.0 2.0 2.0]) == [0.0 0.0 0.0; 0.0 0.0 2.0; 0.0 2.0 0.0; 0.0 2.0 2.0; 2.0 0.0 0.0; 2.0 0.0 2.0; 2.0 2.0 0.0; 2.0 2.0 2.0]
    @test create_grid([0 0 0; 2 2 2], [2.0 2.0 1.0]) == [0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 0.0 2.0; 0.0 2.0 0.0; 0.0 2.0 1.0; 0.0 2.0 2.0; 2.0 0.0 0.0; 2.0 0.0 1.0; 2.0 0.0 2.0; 2.0 2.0 0.0; 2.0 2.0 1.0; 2.0 2.0 2.0]
    @test_throws ArgumentError create_grid([0; 2], [2 2])
    @test_throws ArgumentError create_grid([0 0 0 0; 2 2 2 2])
    @test_throws ArgumentError create_grid([0 0 0; 2 2 2], [2 2])
    @test_throws ArgumentError create_grid([0 0; 2 2], [2 2 2])
    @test_throws ArgumentError create_grid([0 0 0 0; 2 2 2 2], [2 2 2 2])
end