text:
mutable struct NfLattice{T}
  basis::Vector{T}
  discriminant::fmpz
  is_minkowski_exact::Bool
  minkowski_gram_exact::fmpz_mat
  minkowski_gram_scaled::Tuple{Int, fmpz_mat}
  function NfLattice{T}(v::Vector{T}, discriminant::fmpz) where {T <: NumFieldElem}
    n = new{T}()
    n.basis = v
    n.discriminant = discriminant
    is_minkowski_exact = false
    return n
  end
end

function lattice(v::Vector{T}, disc::fmpz; isexact::Bool = false) where T <: NumFieldElem
  @assert !isempty(v)
  L = NfLattice{T}(v, abs(disc))
  L.is_minkowski_exact = isexact
  return L
end

function dim(L::NfLattice)
  return length(L.basis)
end

function basis(L::NfLattice)
  return L.basis
end

function discriminant(L::NfLattice)
  return L.discriminant
end

function nf(L::NfLattice{T}) where T <: NumFieldElem
  return parent(basis(L)[1])
end

function minkowski_matrix(L::NfLattice, p::Int)
  return minkowski_matrix(basis(L), p)
end

#apply the change of basis given by M, creating a new lattice.
function apply(L::NfLattice{T}, t::fmpz_mat) where T <: NumFieldElem
  K = nf(L)
  B = basis(L)
  new_basis = Vector{T}(undef, dim(L))
  tmp = K()
  for i = 1:length(new_basis)
    a = K()
    for j = 1:ncols(t)
      tmp = B[j]*t[i, j]
      a = add!(a, a, tmp)
    end
    new_basis[i] = a
  end
  return lattice(new_basis, discriminant(L), isexact = L.is_minkowski_exact)
end

function minkowski_gram_mat_scaled(L::NfLattice, p::Int)
  if L.is_minkowski_exact
    L.minkowski_gram_exact = _exact_minkowski_matrix(basis(L))
    return L.minkowski_gram_exact
  end
  K = nf(L)
  if isdefined(L, :minkowski_gram_scaled) && L.minkowski_gram_scaled[1] >= p
    A = deepcopy(L.minkowski_gram_mat_scaled[2])
    shift!(A, p - L.minkowski_gram_mat_scaled[1])
  else
    c = minkowski_matrix(L, p)
    B = basis(L)
    d = zero_matrix(FlintZZ, length(B), absolute_degree(K))
    A = zero_matrix(FlintZZ, length(B), length(B))
    round_scale!(d, c, p)
    ccall((:fmpz_mat_gram, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}), A, d)
    shift!(A, -p)
    L.minkowski_gram_scaled = (p, deepcopy(A))
  end
  for i=1:absolute_degree(K)
    fmpz_mat_entry_add_ui!(A, i, i, UInt(nrows(A)))
  end
  return A
end

function weighted_minkowski_gram_scaled(L::NfLattice, v::fmpz_mat, prec::Int)
  c = deepcopy(minkowski_matrix(L, prec))
  mult_by_2pow_diag!(c, v)
  d = zero_matrix(FlintZZ, nrows(c), ncols(c))
  round_scale!(d, c, prec)
  g = zero_matrix(FlintZZ, nrows(c), nrows(c))
  ccall((:fmpz_mat_gram, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}), g, d)
  shift!(g, -prec)
  for i=1:n
    fmpz_mat_entry_add_ui!(g, i, i, UInt(nrows(d)))
  end
  return g
end

function lll(L::NfLattice, weights::fmpz_mat = zero_matrix(FlintZZ, 1, 1); starting_prec::Int = 100 + 25*div(dim(L), 3) + Int(round(log(abs(discriminant(L))))))
  if L.is_minkowski_exact
    M = _exact_minkowski_matrix(basis(L))
    l, v = lll_gram_with_transform(M)
    return FakeFmpqMat(l), v
  end
  prec = starting_prec
  local l1::FakeFmpqMat
  local v::fmpz_mat
  n = dim(L)
  starting_profile = sum(nbits(Hecke.upper_bound(t2(x), fmpz)) for x in basis(L))
  @vprint :LLL 1 "Starting profile: $(starting_profile) \n"
  @vprint :LLL 1 "Target profile: $(nbits(discriminant(L)^2)+divexact(n*(n-1), 2)) \n"
  @vprint :LLL 1 "Starting precision: $(prec) \n"
  while true
    if prec > 2^30
      error("Precision too large!")
    end
    fl, l1, v = _lll(L, weights, prec)
    if fl
      break
    end
    Lnew = apply(L, v)
    new_profile = sum(nbits(Hecke.upper_bound(t2(x), fmpz)) for x in basis(Lnew))
    @vprint :LLL 1 "New profile: $(new_profile) \n"
    if new_profile < starting_profile
      @vprint :LLL 1 "Using transformation!\n"
      l2, v2 = lll(Lnew, weights, starting_prec = ceil(Int, prec*1.5))
      return l2, v2*v
    end
    prec *= 2
    @vprint :LLL 1 "Increasing precision to $(prec) \n"
  end
  return l1, v
end

function _lll(L::NfLattice, weights::fmpz_mat, prec::Int)
  @vprint :LLL 1 "Computing Minkowski Gram matrix with precision $(prec) \n"
  local d::fmpz_mat
  local sv::fmpz
  while true
    try
      if iszero(weights)
        @vtime :LLL 1 d = minkowski_gram_mat_scaled(L, prec)
        sv = fmpz(0)
      else
        @vtime :LLL 1 d = weighted_minkowski_gram_scaled(L, weights, prec)
        sv = max(fmpz(0), sum(weights[1, i] for i=1:ncols(weights)))
      end
      break
    catch e
      prec += 100
    end
  end
  n = dim(L)
  g = identity_matrix(FlintZZ, n)
  g1 = identity_matrix(FlintZZ, n)
  ctx1 = Nemo.lll_ctx(0.4, 0.51, :gram)
  ctx2 = Nemo.lll_ctx(0.99, 0.51, :gram)
  @vtime :LLL 1 ccall((:fmpz_lll, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{Nemo.lll_ctx}), d, g, ctx1)
  @vtime :LLL 1 ccall((:fmpz_lll, libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{Nemo.lll_ctx}), d, g1, ctx2)
  mul!(g, g1, g)
  fl = true
  if nbits(maximum(abs, g)) >  div(prec, 2)
    fl = false
  end
  if fl
    n = nrows(d)
    disc = discriminant(L) * fmpz(2)^(2*sv)
    di = root(disc, n) + 1
    di *= fmpz(2)^(div(n+1,2)) * fmpz(2)^prec
    if compare_index(d, 1, 1, di) > 0
      fl = false
    end
    pr = prod_diagonal(d)
    if pr > fmpz(2)^(div(n*(n-1), 2)) * disc * fmpz(2)^(n*prec)
      fl = false
    end
  end
  return fl, FakeFmpqMat(d, fmpz(2)^prec), g
end

function lll_basis(L::NfLattice{T}) where T
  K = nf(L)
  l, t = lll(L)
  L1 = apply(L, t)
  return basis(L)
end

_abs_disc(O::NfRelOrd) = absolute_discriminant(O)
_abs_disc(I::NfRelOrdIdl) = absolute_discriminant(order(I))*absolute_norm(I)

function _get_nice_basis(OL::T) where T <: Union{NfRelOrdIdl, NfRelOrd}
  L = nf(OL)
  B = pseudo_basis(OL, copy = false)
  ideals = Dict{typeof(B[1][2]), Vector{elem_type(base_field(L))}}()
  for i = 1:length(B)
    if haskey(ideals, B[i][2])
      continue
    end
    ideals[B[i][2]] = lll_basis(B[i][2])
  end
  abs_bas = Vector{elem_type(L)}(undef, absolute_degree(L))
  ind = 1
  for i = 1:length(B)
    J = ideals[B[i][2]]
    for j = 1:length(J)
      abs_bas[ind] = J[j]*B[i][1]
      ind += 1
    end
  end
  return abs_bas
end

function lll_basis(OL::T) where T <: Union{NfRelOrdIdl, NfRelOrd}
  L = nf(OL)
  B = _get_nice_basis(OL)
  isexact = false
  if istotally_real(L)
    isexact = true
  end
  V = lattice(B, _abs_disc(OL), isexact = isexact)
  return lll_basis(V)
end


text:
struct Scene{T<:Tuple}
    shapes::T
end

shapes(scene::Scene) = scene.shapes

min_dist_index(scene::Scene, point::Point{3, T}) where {T<:Real} = findmin([abs(shape(point)) for shape in shapes(scene)])

ext: jl


text:
#!/usr/bin/env julia

home = homedir()
namespace = "pandola"
target = pwd()

# Main call.
main() = begin
    file_name = "index.css"
    pcs_dist_path = "$home/$namespace/pcs/dist/$file_name"
    target_file = "$target/$file_name"

    if target_file |> isfile |> !
        println("\nOK, Checked done.\n")
        println("Make Link:\n$pcs_dist_path ->\n$target_file")
        symlink(pcs_dist_path, target_file)
        println("\nOK!\n")
    else
        println("\nSorry, Files existed. Please use pandola-css-clear remove files.\n")
    end
end

main()

max_stars_repo_name: yuffiy/pandola

max_issues_repo_name: HairyRabbit/pandola

max_forks_repo_name: HairyRabbit/pandola


text:
using GRAMPC
using GRAMPC_jll
using Test

####################################################################
# Introduce an error and see if it is reported by the library
####################################################################
@test_throws ErrorException ccall( (:printError, libgrampc_julia), Cvoid, (Cstring,), "Test error" )
@test_throws ErrorException ccall( (:printErrorAddString, libgrampc_julia), Cvoid, (Cstring, Cstring), "Error", "Message" )

####################################################################
# Test warning reporting
####################################################################
warnCString = Vector{UInt8}(undef, 1024)

# Ensure our clear function works
GRAMPC._clear_warning()
@test ccall( (:grampcjl_get_last_warning, libgrampc_julia), Cint, (Ptr{UInt8},), warnCString ) == 0

# Introduce a warning and see if it is reported by the library (these are only done by polling)
ccall( (:printWarningAddString, libgrampc_julia), Cvoid, (Cstring, Cstring), "Warning", "Message" )
@test ccall( (:grampcjl_get_last_warning, libgrampc_julia), Cint, (Ptr{UInt8},), warnCString ) == 1
@test unsafe_string( pointer( warnCString ) ) == "Warning: Message"

# Make sure the warning was cleared after we read it
@test ccall( (:grampcjl_get_last_warning, libgrampc_julia), Cint, (Ptr{UInt8},), warnCString ) == 0

# Test the warning macro
ccall( (:printWarningAddString, libgrampc_julia), Cvoid, (Cstring, Cstring), "Warning", "Message" )
@test_logs (:warn, "Warning: Message") GRAMPC.@check_warning

# Make sure the warning was cleared after we read it
@test ccall( (:grampcjl_get_last_warning, libgrampc_julia), Cint, (Ptr{UInt8},), warnCString ) == 0


text:
######## Internal Simulation Object Structs ########
abstract type AbstractStage end

mutable struct _Stage{M<:AbstractOperationsProblem} <: AbstractStage
    key::Int64
    reference::OperationsTemplate
    op_problem::Type{M}
    sys::PSY.System
    canonical::Canonical
    optimizer::JuMP.OptimizerFactory
    executions::Int64
    execution_count::Int64
    interval::Dates.Period
    chronology_ref::Dict{Int64, <:Chronology}
    ini_cond_chron::Union{<:Chronology, Nothing}
    cache::Dict{Type{<:AbstractCache}, AbstractCache}

    function _Stage(key::Int64,
                    reference::OperationsTemplate,
                    op_problem::Type{M},
                    sys::PSY.System,
                    canonical::Canonical,
                    optimizer::JuMP.OptimizerFactory,
                    executions::Int64,
                    interval::Dates.Period,
                    chronology_ref::Dict{Int64, <:Chronology},
                    cache::Vector{<:AbstractCache}) where M <: AbstractOperationsProblem

    ini_cond_chron = get(chronology_ref, 0, nothing)
    if !isempty(get_initial_conditions(canonical))
        if isnothing(ini_cond_chron)
            @warn("Initial Conditions chronology set for Stage $(key) which contains Initial conditions")
        end
    end

    pop!(chronology_ref, 0, nothing)

    cache_dict = Dict{Type{<:AbstractCache}, AbstractCache}()
    for c in cache
        cache_dict[typeof(c)] = c
    end

    new{M}(key,
           reference,
           op_problem,
           sys,
           canonical,
           optimizer,
           executions,
           0,
           interval,
           chronology_ref,
           ini_cond_chron,
           cache_dict)

    end

end

######## Exposed Structs to define a Simulation Object ########

mutable struct Stage <: AbstractStage
    op_problem::Type{<:AbstractOperationsProblem}
    model::OperationsTemplate
    initial_time::Dates.DateTime
    horizon::Int64
    interval::Dates.Period
    execution_count::Int64
    sys::PSY.System
    optimizer::JuMP.OptimizerFactory
    chronology_ref::Dict{Int64, <:Chronology}
    cache::Vector{<:AbstractCache}

    function Stage(::Type{M},
                   model::OperationsTemplate,
                   initial_time::Dates.DateTime,
                   horizon::Int64,
                   interval::Dates.Period,
                   execution_count::Int64,
                   sys::PSY.System,
                   optimizer::JuMP.OptimizerFactory,
                   chronology_ref=Dict{Int, <:Chronology}(),
                   cache::Vector{<:AbstractCache}=Vector{AbstractCache}()) where M<:AbstractOperationsProblem

        new(M,
            model,
            initial_time,
            horizon,
            interval,
            execution_count,
            sys,
            optimizer,
            chronology_ref,
            cache)
    end

end

function Stage(::Type{M},
               model::OperationsTemplate,
               horizon::Int64,
               interval::Dates.Period,
               execution_count::Int64,
               sys::PSY.System,
               optimizer::JuMP.OptimizerFactory,
               chronology_ref::Dict{Int64, <:Chronology},
               cache::Union{Nothing, AbstractCache}=nothing;
               kwargs...) where M<:AbstractOperationsProblem

    initial_time = get(kwargs, :initial_time, PSY.get_forecast_initial_times(sys)[1])
    cacheinput = isnothing(cache) ? Vector{AbstractCache}() : [cache]
    return Stage(M,
                model,
                initial_time,
                horizon,
                interval,
                execution_count,
                sys,
                optimizer,
                chronology_ref,
                cacheinput)

end

function Stage(model::OperationsTemplate,
               horizon::Int64,
               interval::Dates.Period,
               execution_count::Int64,
               sys::PSY.System,
               optimizer::JuMP.OptimizerFactory,
               chronology_ref::Dict{Int64, <:Chronology},
               cache::Union{Nothing, AbstractCache}=nothing;
               kwargs...)

    initial_time = get(kwargs, :initial_time, PSY.get_forecast_initial_times(sys)[1])
    return Stage(GenericOpProblem,
                model,
                horizon,
                interval,
                execution_count,
                sys,
                optimizer,
                chronology_ref,
                cache; initial_time = initial_time, kwargs...)

end

get_execution_count(s::S) where S <: AbstractStage = s.execution_count
get_sys(s::S) where S <: AbstractStage = s.sys
get_chronology_ref(s::S) where S <: AbstractStage = s.chronology_ref

get_template(s::Stage) = s.model


text:
# AUTO GENERATED FILE - DO NOT EDIT

export daq_thermometer

"""
    daq_thermometer(;kwargs...)

A Thermometer component.
A thermometer component that
fills to a value between some
range
Keyword arguments:
- `id` (String; optional): The ID used to identify this compnent in Dash callbacks
- `base` (Real; optional): Base to be used in logarithmic scale.
- `className` (String; optional): Class to apply to the root component element.
- `color` (String; optional): The color of the thermometer fill/current value text
- `height` (Real; optional): The height of the thermometer in pixels
- `label` (optional): Description to be displayed alongside the control. To control styling, pass an object with label and style properties.. label has the following type: String | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `label` (String; optional)
- `labelPosition` (a value equal to: 'top', 'bottom'; optional): Where the component label is positioned.
- `logarithmic` (Bool; optional): If set to true, a logarithmic scale will be
used.
- `max` (Real; optional): The maximum value of the thermometer. If logarithmic,
represents the maximum exponent.
- `min` (Real; optional): The minimum value of the thermometer. If logarithmic,
represents the minimum exponent.
- `scale` (optional): Configuration for the component scale.. scale has the following type: lists containing elements 'start', 'interval', 'labelInterval', 'custom'.
Those elements have the following types:
  - `start` (Real; optional): Value to start the scale from. Defaults
to min.
  - `interval` (Real; optional): Interval by which the scale goes up. Attempts
to dynamically divide min-max range by
default.
  - `labelInterval` (Real; optional): Interval by which labels are added to
scale marks. Defaults to 2 (every other
mark has a label).
  - `custom` (optional): Custom scale marks. The key determines the position
and the value determines what will show. If you want
to set the style of a specific mark point, the value
should be an object which contains style and label
properties. custom has the following type: Real | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (String; optional)
  - `label` (String; optional)
- `showCurrentValue` (Bool; optional): If true, the current value of the
thermometer will be displayed
- `style` (Dict; optional): Style to apply to the root component element.
- `theme` (Dict; optional): Theme configuration to be set by a ThemeProvider
- `units` (String; optional): Label for the current value
- `value` (Real; optional): The value of thermometer. If logarthmic, the value
displayed will be the logarithm of the inputted value.
- `width` (Real; optional): The width of the thermometer in pixels
"""
function daq_thermometer(; kwargs...)
        available_props = Symbol[:id, :base, :className, :color, :height, :label, :labelPosition, :logarithmic, :max, :min, :scale, :showCurrentValue, :style, :theme, :units, :value, :width]
        wild_props = Symbol[]
        return Component("daq_thermometer", "Thermometer", "dash_daq", available_props, wild_props; kwargs...)
end


text:
module Stride

export stride_run,
       stride_update_path,
       stride_parse,
       exists_in_path

type Stride_guts
    stride_path :: AbstractString
    delete_tmp_files :: Bool
end

"""
SSDataResidue holds the information related to the secondary structure
associated to a residue.

The fields are as follows:
| Field Name   | Field Type     | Description                        |
| ---          | ---            | ---                                |
| name         | AbstractString | The 3 letter residue name          |
| ss_type      | Char           | The one letter secondary structure |
| ss_type_full | AbstractString | Full SS name                       |
| phi          | Float64        | The phi                            |
| psi          | Float64        | and psi dihedral angle             |
| area         | Float64        | Solvent accessible area            |

The following Secondary Structures are defined using the [DSSP8 standard](https://zhanglab.ccmb.med.umich.edu/literature/dssp.pdf) and defined as follows:
| --- | ---            |
| H   | Alpha Helix    |
| B   | Beta Bridge    |
| E   | Strand         |
| T   | H Bounded Turn |
| I   | Pi Helix       |
| G   | 3/10 Helix     |
| S   | Bend           |
| C   | Random Coil    |
"""
type SSDataResidue
    name         :: AbstractString
    ss_type      :: Char
    ss_type_full :: AbstractString
    phi          :: Float64
    psi          :: Float64
    area         :: Float64
end

"""
SSData holds the information related to a given protein
"""
type SSData
    chain :: Array{SSDataResidue, 1}
end

stride_state = Stride_guts("", true)

"""
    stride_run(in_path :: AbstractString; out_path::AbstractString="")

Runs Stride for a given pdb file. if out_path is set then the temporary
Stride file generated will be kept. The function returns a SSData object.
"""
function stride_run(in_path :: AbstractString; out_path::AbstractString="")
    if length(stride_state.stride_path) == 0
        if exists_in_path("stride")
            stride_path = "stride"
        else
            error("Stride not found in PATH!\nPlease add stride to PATH" *
                  " or run stride_update_path(\"pathtostride\")")
        end
    else
        stride_path = stride_state.stride_path
    end

    if out_path == ""
        out_path = tempname()
    else
        stride_state.delete_tmp_files = false
    end

    args = [in_path, "-f$(out_path)"]
    cmd_string = `$(stride_path) $args`
    run(cmd_string)

    ss = stride_parse(out_path)

    if stride_state.delete_tmp_files
        rm(out_path)
    end

    return ss
end

"""
    stride_update_path(path :: AbstractString)

stride_update_path can be used to change the path where the package will
look for the Stride binary.
"""
function stride_update_path(path :: AbstractString)
    stride_state.stride_path = path
end

"""
    stride_delete_tmp_files(delete :: Bool)

stride_delete_tmp_files recieves a boolean value that determines if the tempfiles
will be deleted after use. Default option is to delete it.
"""
function stride_delete_tmp_files(delete :: Bool)
    stride_state.delete_tmp_files = delete
end

"""
    stride_parse is reponsible for parsing the output of Stride into a SSData object
"""
function stride_parse(path :: AbstractString)
    data = readdlm(path)
    nlines, ncolumns = size(data)

    ss = SSData([])

    for i in 1:nlines
        if data[i] == "REM"
            continue
        elseif data[i] == "LOC"
            continue
        elseif data[i] == "ASG"
            push!(ss.chain, SSDataResidue(data[i, 2], data[i, 6][1],
                                            data[i, 7], data[i, 8],
                                            data[i, 9], data[i, 10]
                                            ))
        end
    end

    return ss
end

"""
    exists_in_path(prog::String)

Determine whether the given program name or path is executable using the current user's
permissions. This is roughly equivalent to querying `which program` at the command line
and checking that a result is found, but no shelling out occurs.
"""
function exists_in_path(prog::String) # Code taken from Alex Arslan with permision
    X_OK = 1 << 0 # Taken from unistd.h
    # If prog has a slash, we know the user wants to determine whether the given
    # file exists and is executable
    if '/' in prog
        isfile(prog) || return false
        return ccall(:access, Cint, (Ptr{UInt8}, Cint), prog, X_OK) == 0
    end
    path = get(ENV, "PATH", "")
    # Something is definitely wrong if the user's path is empty...
    @assert !isempty(path)
    sep = ':'
    for dir in split(path, sep), file in readdir(dir)
        if file == prog
            p = joinpath(dir, file)
            @assert isfile(p)
            return ccall(:access, Cint, (Ptr{UInt8}, Cint), p, X_OK) == 0
        end
    end
    false
end
isexecutable(prog::String) = isexecutable(Sys.KERNEL, prog)

end


text:
using Turing, Distributions, DifferentialEquations

# Import MCMCChain, Plots, and StatsPlots for visualizations and diagnostics.
using MCMCChains, Plots, StatsPlots

# Set a seed for reproducibility.
using Random
Random.seed!(14);
using MCMCTempering

function lotka_volterra(du,u,p,t)
    x, y = u
    α, β, γ, δ  = p
    du[1] = (α - β*y)x # dx =
    du[2] = (δ*x - γ)y # dy =
end

p = [1.5, 1.0, 3.0, 1.0]
u0 = [1.0,1.0]
prob1 = ODEProblem(lotka_volterra,u0,(0.0,10.0),p)
sol = solve(prob1,Tsit5())
plot(sol)

sol1 = solve(prob1,Tsit5(),saveat=0.1)
odedata = Array(sol1) + 0.8 * randn(size(Array(sol1)))
plot(sol1, alpha = 0.3, legend = false); scatter!(sol1.t, odedata')

Turing.setadbackend(:forwarddiff)

@model function fitlv(data, prob1)
    σ ~ InverseGamma(2, 3) # ~ is the tilde character
    α ~ truncated(Normal(1.5,0.5),0.5,2.5)
    β ~ truncated(Normal(1.0,0.5),0,2)
    γ ~ truncated(Normal(3.0,0.5),1,4)
    δ ~ truncated(Normal(1.0,0.5),0,2)

    p = [α,β,γ,δ]
    prob = remake(prob1, p=p)
    predicted = solve(prob,Tsit5(),saveat=0.1)

    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end

model = fitlv(odedata, prob1)

# This next command runs 3 independent chains without using multithreading.
chain1_mh = sample(model, MH(), MCMCThreads(), 10000, 4)
chain1_nuts = sample(model, NUTS(.65), 1000)
# chain1_mh_test = mapreduce(c -> sample(model, NUTS(), 1000), chainscat, 1:3)

chain2_mh = sample(model, Tempered(MH(), 2), MCMCThreads(), 10000, 4)
chain2_nuts = sample(model, Tempered(NUTS(.65), 2), 1000)

chain3_mh = sample(model, Tempered(MH(), 3), MCMCThreads(), 10000, 30)
chain3_nuts = sample(model, Tempered(NUTS(.65), 3), MCMCThreads(), 1000, 30)

chain4_mh = sample(model, Tempered(MH(), 4), MCMCThreads(), 10000, 30)
chain4_nuts = sample(model, Tempered(NUTS(.65), 4), MCMCThreads(), 1000, 30)

plot_swaps(chain2_mh)

plot(chain1_mh)
plot(chain2_mh)

interchain_stats(chain1_mh)
interchain_stats(chain2_mh)

# Pumas example

function theop_model_Depots1Central1(du, u, p, t)
    Depot, Central = u
    Ka, CL, Vc = p
    du[1] = -Ka * Depot # d Depot =
    du[2] =  Ka * Depot - (CL / Vc) * Central # d Central =
end

u0 = [1.0, 1.0]
p = [2.0, 0.2, 0.8, 2.0]
prob = ODEProblem(theop_model_Depots1Central1,u0,(0.0, 10.0),p)
sol = solve(prob, Tsit5())
plot(sol)

@model function theopmodel_bayes(dv, SEX, WT)

    N = length(dv)

    θ ~ arraydist(truncated.(Normal.([2.0, 0.2, 0.8, 2.0], 1.0), 0.0, 10.0))

    ωKa ~ Gamma(1.0, 0.2)
    ωCL ~ Gamma(1.0, 0.2)
    ωVc ~ Gamma(1.0, 0.2)

    σ ~ Gamma(1.0, 0.5)

    ηKa ~ filldist(Normal(0.0, ωKa), N)
    ηCL ~ filldist(Normal(0.0, ωCL), N)
    ηVc ~ filldist(Normal(0.0, ωVc), N)

    for i in 1:N
        Ka = (SEX[i] == 1 ? θ[1] : θ[4]) * exp(ηKa[i])
        CL = θ[2]*(WT[i]/70) * exp(ηCL[i])
        Vc = θ[3] * exp(ηVc[i])

        p = [Ka, CL, Vc]
        prob = remake(prob1, p=p)
        predicted = solve(prob, Tsit5(), saveat=0.1)

        μ[i] = predicted[i,2] / Vc
        dv[i] .~ Normal.(μ, σ)
    end
    dv

end

using Pumas

# BayesNeuralODE example

using BayesNeuralODE

N = 1
prior_std = likelihood_std = 1.0
model = BNO.generate_turing_model(:spiral, N, prior_std, likelihood_std)

bno_chain_1 = sample(model, NUTS(.6), 100)

bno_chain_2 = sample(model, Tempered(NUTS(.6), 4), MCMCThreads(), 1000, 4)


text:
include(joinpath(pwd(), "src/models/box.jl"))
include(joinpath(pwd(), "src/constraints/contact.jl"))

# Horizon
T = 26

# Time step
tf = 2.5
h = tf / (T-1)

# Bounds

_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 0.0
_ul = zeros(model.m)
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
mrp = MRP(UnitQuaternion(RotY(pi / 10.0)*RotX(pi / 15.0)))

q1 = [0.0, 0.0, 2.5, mrp.x, mrp.y, mrp.z]
v1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
v2 = v1 - G_func(model,q1) * h
q2 = q1 + 0.5 * h * (v1 + v2)

x1 = [q1; q2]

xl, xu = state_bounds(model, T, x1 = x1)

# Objective
obj = PenaltyObjective(1000.0, model.m)

# Constraints
con_contact = contact_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact
               )

# Trajectory initialization
X0 = [0.01 * rand(model.n) for t = 1:T] #linear_interp(x1, x1, T) # linear interpolation on state
U0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

@time Z̄ = solve(prob, copy(Z0), tol = 1.0e-3, c_tol = 1.0e-3)

check_slack(Z̄, prob)
X̄, Ū = unpack(Z̄, prob)

include(joinpath(pwd(), "src/models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model, state_to_configuration(X̄), Δt = h)

obj = Rect(Vec(-1.0 * model.r,
    -1.0 * model.r,
    -1.0 * model.r),
    Vec(2.0 * model.r, 2.0 * model.r, 2.0 * model.r))

setobject!(vis, obj)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
"""
"""
function multilabelprobabilitiestopredictions(
        probabilitiesassoc::AbstractDict;
        float_type::Type{<:AbstractFloat} = Float64,
        )
    result = Dict()
    label_names = sort(unique(collect(keys(probabilitiesassoc))))
    for j = 1:length(label_names)
        result[label_names[j]] = single_labelprobabilitiestopredictions(
            probabilitiesassoc[label_names[j]];
            float_type=float_type,
            )
    end
    result = Dict()
    return result
end

"""
"""
function multilabelprobabilitiestopredictions(
        probabilitiesassoc::AbstractDict,
        thresholds::AbstractDict;
        float_type::Type{<:AbstractFloat} = Float64,
        )
    result = Dict()
    label_names = sort(unique(collect(keys(probabilitiesassoc))))
    for j = 1:length(label_names)
        if haskey(thresholds, label_names[j])
            (positive_class, threshold) = thresholds[label_names[j]]
            result[label_names[j]] = single_labelprobabilitiestopredictions(
                probabilitiesassoc[label_names[j]],
                positive_class,
                threshold;
                float_type=float_type,
                )
        else
            result[label_names[j]] = single_labelprobabilitiestopredictions(
                probabilitiesassoc[label_names[j]];
                float_type=float_type,
                )
        end
    end
    result = Dict()
    return result
end

const probabilitiestopredictions = multilabelprobabilitiestopredictions

"""
"""
function single_labelprobabilitiestopredictions(
        probabilitiesassoc::AbstractDict;
        float_type::Type{<:AbstractFloat} = Float64,
        )
    classes = sort(unique(collect(keys(probabilitiesassoc))))
    numclasses = length(classes)
    numrows = size(probabilitiesassoc[classes[1]], 1)
    probabilitiesmatrix = Matrix{float_type}(
        undef,
        numrows,
        numclasses,
        )
    for j = 1:numclasses
        probabilitiesmatrix[:, j] = float_type.(probabilitiesassoc[classes[j]])
    end
    predictionsvector = Vector{String}(undef, numrows,)
    for i = 1:numrows
        predictionsvector[i] =
            string(classes[argmax(probabilitiesmatrix[i, :])])
    end
    return predictionsvector
end

function single_labelprobabilitiestopredictions(
        probabilitiesassoc::AbstractDict,
        positive_class,
        threshold::AbstractFloat;
        float_type::Type{<:AbstractFloat} = Float64,
        )
    classes = sort(unique(collect(keys(probabilitiesassoc))))
    numclasses = length(classes)
    numrows = size(probabilitiesassoc[classes[1]], 1)
    probabilitiesmatrix = Matrix{float_type}(
        undef,
        numrows,
        numclasses,
        )
    for j = 1:numclasses
        probabilitiesmatrix[:, j] = float_type.(probabilitiesassoc[classes[j]])
    end
    predictionsvector = Vector{String}(undef, numrows,)
    if numclasses == 2
        if positive_class in classes
            use_threshold = true
        else
            @error("", positive_class, classes)
            error("positive_class is not in the list of classes, so ignoring threshold")
            use_threshold = false
        end
    else
        error("numclasses is not 2, so ignoring threshold")
        use_threshold = false
    end
    if use_threshold
        positive_class_column::Int = findfirst(classes .== positive_class)
        negative_class = first(setdiff(classes, [positive_class]))
        for i = 1:numrows
            if probabilitiesmatrix[i, positive_class_column] > threshold
                predictionsvector[i] = string(positive_class)
            else
                predictionsvector[i] = string(negative_class)
            end
        end
    else
        for i = 1:numrows
            predictionsvector[i] =
                string(classes[argmax(probabilitiesmatrix[i, :])])
        end
    end
    return predictionsvector
end


text:
# Alignment
# =========
#
# Sequence alignment type.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Defines how to align a given sequence onto a reference sequence.
The alignment is represented as a sequence of elementary operations (match, insertion, deletion etc)
anchored to specific positions of the input and reference sequence.
"""
struct Alignment
    anchors::Vector{AlignmentAnchor}
    firstref::Int
    lastref::Int

    # https://github.com/JuliaLang/julia/issues/16730
    @doc """
        Alignment(anchors::Vector{AlignmentAnchor}, check=true)

    Create an alignment object from a sequence of alignment anchors.
    """ ->
    function Alignment(anchors::Vector{AlignmentAnchor}, check::Bool=true)
        # optionally check coherence of the anchors
        if check
            check_alignment_anchors(anchors)
        end

        # compute first and last aligned reference positions
        firstref = 0
        for i in 1:length(anchors)
            if ismatchop(anchors[i].op)
                firstref = anchors[i-1].refpos + 1
                break
            end
        end

        lastref = 0
        for i in length(anchors):-1:1
            if ismatchop(anchors[i].op)
                lastref = anchors[i].refpos
                break
            end
        end

        return new(anchors, firstref, lastref)
    end
end

"""
    Alignment(cigar::AbstractString, seqpos=1, refpos=1)

Make an alignment object from a CIGAR string.

`seqpos` and `refpos` specify the starting positions of two sequences.
"""
function Alignment(cigar::AbstractString, seqpos::Int=1, refpos::Int=1)
    # path starts prior to the first aligned position pair
    seqpos -= 1
    refpos -= 1
    alnpos = 0

    n = 0
    anchors = AlignmentAnchor[AlignmentAnchor(seqpos, refpos, alnpos, OP_START)]
    for c in cigar
        if isdigit(c)
            n = n * 10 + convert(Int, c - '0')
        else
            if n == 0
                error("CIGAR operations must be prefixed by a positive integer.")
            end
            op = Operation(c)
            if ismatchop(op)
                seqpos += n
                refpos += n
            elseif isinsertop(op)
                seqpos += n
            elseif isdeleteop(op)
                refpos += n
            else
                error("The $(op) CIGAR operation is not yet supported.")
            end
            alnpos += n

            push!(anchors, AlignmentAnchor(seqpos, refpos, alnpos, op))
            n = 0
        end
    end

    return Alignment(anchors)
end

function Base.:(==)(a::Alignment, b::Alignment)
    return a.anchors == b.anchors && a.firstref == b.firstref && a.lastref == b.lastref
end

function Base.show(io::IO, aln::Alignment)
    println(io, summary(aln), ':')
    println(io, "  aligned range:")
    println(io, "    seq: ", first(aln.anchors).seqpos, '-', last(aln.anchors).seqpos)
    println(io, "    ref: ", first(aln.anchors).refpos, '-', last(aln.anchors).refpos)
      print(io, "  CIGAR string: ", cigar(aln))
end

# generic function for mapping between sequence, reference and alignment positions
# getsrc specifies anchor source position getter
# getdest specifies anchor destination position getter
function pos2pos(aln::Alignment, i::Integer,
                 srcpos::Function, destpos::Function)::Tuple{Int,Operation}
    idx = findanchor(aln, i, srcpos)
    if idx == 0
        if srcpos === seqpos
            throw(ArgumentError("invalid sequence position: $i"))
        elseif srcpos === refpos
            throw(ArgumentError("invalid reference position: $i"))
        elseif srcpos === alnpos
            throw(ArgumentError("invalid alignment position: $i"))
        else
            throw(ArgumentError("Unknown position getter: $srcpos"))
        end
    end
    anchor = aln.anchors[idx]
    pos = destpos(anchor)
    if ismatchop(anchor.op) ||
        ((srcpos === alnpos) && ((destpos === seqpos) && isinsertop(anchor.op) || (destpos === refpos) && isdeleteop(anchor.op))) ||
        ((destpos === alnpos) && ((srcpos === seqpos) && isinsertop(anchor.op) || (srcpos === refpos) && isdeleteop(anchor.op)))
        pos += i - srcpos(anchor)
    end
    return pos, anchor.op
end

"""
    seq2ref(aln::Union{Alignment, AlignedSequence, PairwiseAlignment}, i::Integer)::Tuple{Int,Operation}

Map a position `i` from sequence to reference.
"""
seq2ref(aln::Alignment, i::Integer) = pos2pos(aln, i, seqpos, refpos)

"""
    ref2seq(aln::Union{Alignment, AlignedSequence, PairwiseAlignment}, i::Integer)::Tuple{Int,Operation}

Map a position `i` from reference to sequence.
"""
ref2seq(aln::Alignment, i::Integer) = pos2pos(aln, i, refpos, seqpos)

"""
    seq2aln(aln::Union{Alignment, AlignedSequence, PairwiseAlignment}, i::Integer)::Tuple{Int,Operation}

Map a position `i` from the input sequence to the alignment sequence.
"""
seq2aln(aln::Alignment, i::Integer) = pos2pos(aln, i, seqpos, alnpos)

"""
    ref2aln(aln::Union{Alignment, AlignedSequence, PairwiseAlignment}, i::Integer)::Tuple{Int,Operation}

Map a position `i` from the reference sequence to the alignment sequence.
"""
ref2aln(aln::Alignment, i::Integer) = pos2pos(aln, i, refpos, alnpos)

"""
    aln2seq(aln::Union{Alignment, AlignedSequence, PairwiseAlignment}, i::Integer)::Tuple{Int,Operation}

Map a position `i` from the alignment sequence to the input sequence.
"""
aln2seq(aln::Alignment, i::Integer) = pos2pos(aln, i, alnpos, seqpos)

"""
    aln2ref(aln::Union{Alignment, AlignedSequence, PairwiseAlignment}, i::Integer)::Tuple{Int,Operation}

Map a position `i` from the alignment sequence to the reference sequence.
"""
aln2ref(aln::Alignment, i::Integer) = pos2pos(aln, i, alnpos, refpos)

"""
    cigar(aln::Alignment)

Make a CIGAR string encoding of `aln`.

This is not entirely lossless as it discards the alignments start positions.
"""
function cigar(aln::Alignment)
    anchors = aln.anchors
    if isempty(anchors)
        return ""
    end
    seqpos = anchors[1].seqpos
    refpos = anchors[1].refpos
    @assert anchors[1].op == OP_START
    out = IOBuffer()
    for i in 2:length(anchors)
        n = max(anchors[i].seqpos - anchors[i-1].seqpos,
                anchors[i].refpos - anchors[i-1].refpos)
        print(out, n, convert(Char, anchors[i].op))
    end
    return String(take!(out))
end

# Check validity of a sequence of anchors.
function check_alignment_anchors(anchors)
    if isempty(anchors)
        # empty alignments are valid, representing an unaligned sequence
        return
    end

    if anchors[1].op != OP_START
        error("Alignments must begin with on OP_START anchor.")
    end

    for i in 2:lastindex(anchors)
        @inbounds acur, aprev = anchors[i], anchors[i-1]
        if acur.refpos < aprev.refpos || acur.seqpos < aprev.seqpos || acur.alnpos < aprev.alnpos
            error("Alignment anchors must be sorted.")
        end

        op = acur.op
        if !isvalid(op)
            error("Anchor at index $(i) has an invalid operation.")
        end

        # reference skip/delete operations
        if isdeleteop(op)
            if acur.seqpos != aprev.seqpos
                error("Invalid anchor sequence positions for reference deletion.")
            end
            if acur.alnpos - aprev.alnpos != acur.refpos - aprev.refpos
                error("Invalid anchor reference positions for reference deletion.")
            end
        # reference insertion operations
        elseif isinsertop(op)
            if acur.refpos != aprev.refpos
                error("Invalid anchor reference positions for reference insertion.")
            end
            if acur.alnpos - aprev.alnpos != acur.seqpos - aprev.seqpos
                error("Invalid anchor sequence positions for reference deletion.")
            end
        # match operations
        elseif ismatchop(op)
            if (acur.refpos - aprev.refpos != acur.seqpos - aprev.seqpos) ||
               (acur.alnpos - aprev.alnpos != acur.seqpos - aprev.seqpos)
               error("Invalid anchor positions for match operation.")
            end
        end
    end
end

# find the index of the first anchor that satisfies `i ≤ pos(anchor)`
function findanchor(aln::Alignment, i::Integer, pos::Function)
    anchors = aln.anchors
    lo = 1
    hi = lastindex(anchors)
    @inbounds if !(pos(anchors[lo]) < i ≤ pos(anchors[hi]))
        return 0
    end
    # binary search
    @inbounds while hi - lo > 2
        m = (lo + hi) >> 1
        if pos(anchors[m]) < i
            lo = m
        else  # i ≤ pos(anchors[m])
            hi = m
        end
        # invariant (activate this for debugging)
        #@assert pos(anchors[lo]) < i ≤ pos(anchors[hi])
    end
    # linear search
    @inbounds for j in lo+1:hi
        if i ≤ pos(aln.anchors[j])
            return j
        end
    end
    # do not reach here
    @assert false
    return 0
end


text:
# Expected reconstruction error for Bernoulli output with diagonal covariance

@defstruct BernoulliReconLossLayer Layer (
    name :: AbstractString = "bernoulli-recon-loss",
    (weight :: AbstractFloat = 1.0, weight >= 0),
    (bottoms :: Vector{Symbol} = Symbol[:p, :x], length(bottoms) == 2),
)

@characterize_layer(BernoulliReconLossLayer,
    has_loss  => true,
    can_do_bp => true,
    is_sink   => true,
    has_stats => true,
)

type BernoulliReconLossLayerState{T, B<:Blob} <: LayerState
    layer      :: BernoulliReconLossLayer
    loss       :: T
    loss_accum :: T
    n_accum    :: Int

    tmp_blobs :: Dict{Symbol, B}
end

function setup(backend::CPUBackend, layer::BernoulliReconLossLayer,
            inputs::Vector{Blob}, diffs::Vector{Blob})
    data_type = eltype(inputs[1])
    return BernoulliReconLossLayerState(layer, zero(data_type), zero(data_type),
                                                0, Dict{Symbol, CPUBlob}())
end

function shutdown(backend::Backend, state::BernoulliReconLossLayerState)
    for blob in values(state.tmp_blobs)
        destroy(blob)
    end
end

function reset_statistics(state::BernoulliReconLossLayerState)
    state.n_accum = 0
    state.loss_accum = zero(typeof(state.loss_accum))
end

function dump_statistics(storage, state::BernoulliReconLossLayerState, show::Bool)
    update_statistics(storage, "$(state.layer.name)-recon-loss", state.loss_accum)

    if show
      loss = @sprintf("%.4f", state.loss_accum)
      @info("  bernoulli-recon-loss (avg over $(state.n_accum)) = $loss")
    end
end

function forward(backend::CPUBackend, state::BernoulliReconLossLayerState,
            inputs::Vector{Blob})
    data_type = eltype(inputs[1])
    N = get_num(inputs[1])
    ND = length(inputs[1])
    p = inputs[1].data
    x = inputs[2].data

    state.loss = 0.
    for i in 1:ND
        state.loss -= x[i]log(p[i])
        state.loss -= (1 - x[i])log(1 - p[i])
    end
    state.loss *= state.layer.weight / N

    # accumulate statistics
    state.loss_accum *= state.n_accum
    state.loss_accum += state.loss * N
    state.loss_accum /= state.n_accum + N

    state.n_accum += N
end

function backward(backend::CPUBackend, state::BernoulliReconLossLayerState,
            inputs::Vector{Blob}, diffs::Vector{Blob})
    data_type = eltype(inputs[1])
    N = get_num(inputs[1])
    ND = length(inputs[1])
    p = inputs[1].data
    x = inputs[2].data

    if isa(diffs[1], CPUBlob)
        for i in 1:ND
            diffs[1].data[i] = x[i] / p[i] - (1 - x[i]) / (1 - p[i])
        end
        diffs[1].data[:] *= -state.layer.weight / N
    end

    if isa(diffs[2], CPUBlob)
        raise("last blob should be a data blob")
    end
end


text:
# Transmission & distribution coupling functions

## Functions replacing PowerModels or InfrastructureModels functions

""
function run_model(
        t_data::Dict{String,<:Any},
        d_data::Dict{String,<:Any},
        t_model_type::Type,
        d_model_type::Type{BF},
        optimizer::Union{MathOptInterface.AbstractOptimizer, MathOptInterface.OptimizerWithAttributes},
        build_method::Function;
        t_ref_extensions::Vector{<:Function} = Function[],
        d_ref_extensions::Vector{<:Function} = Function[],
        t_solution_processors::Vector{<:Function} = Function[],
        d_solution_processors::Vector{<:Function} = Function[],
        kwargs...
    ) where BF <: _PM.AbstractBFModel

    # Check that transmission and distribution network ids are different.
    if !isempty(intersect(Set(id for (id,nw) in t_data["nw"]), Set(id for (id,nw) in d_data["nw"])))
        Memento.error(_LOGGER, "Transmission and distribution data contain networks having the same IDs.")
    end

    # Instantiate models.
    start_time = time()
    t_pm, d_pm = instantiate_model(t_data, d_data, t_model_type, d_model_type, build_method; t_ref_extensions, d_ref_extensions, kwargs...)
    Memento.debug(_LOGGER, "combined T&D model build time: $(time() - start_time)")

    start_time = time()

    # Solve the optimization model and store the transmission result.
    t_result = _IM.optimize_model!(t_pm; optimizer, solution_processors=t_solution_processors)

    # Build the distribution result using the same model as above.
    d_result = _IM.build_result(d_pm, t_result["solve_time"]; solution_processors=d_solution_processors)

    # The asymmetric code above for building results produces inaccurate debugging messages;
    # this behavior can be fixed by writing a custom optimize_model!() that takes 2 models.

    Memento.debug(_LOGGER, "combined T&D model solution time: $(time() - start_time)")

    # Combine the result objects.
    result = t_result # All fields are pairwise equal, except for "solution".
    solution = result["solution"]
    for (k,v) in d_result["solution"]
        if k == "nw"
            # Merge transmission and distribution "nw" fields.
            # Network ids do not clash (checked before model instantiation).
            solution["nw"] = merge(solution["nw"], v)
        else
            if solution[k] != v
                Memento.warning(_LOGGER, "Transmission and distribution solutions differ on key \"$k\"; only transmission value is kept.")
            end
        end
    end

    return result
end

""
function instantiate_model(
        t_data::Dict{String,<:Any},
        d_data::Dict{String,<:Any},
        t_model_type::Type,
        d_model_type::Type{BF},
        build_method::Function;
        t_ref_extensions::Vector{<:Function} = Function[],
        d_ref_extensions::Vector{<:Function} = Function[],
        kwargs...
    ) where BF <: _PM.AbstractBFModel

    # Instantiate the transmission PowerModels struct, without building the model.
    t_pm = _PM.instantiate_model(t_data, t_model_type, method->nothing; ref_extensions=t_ref_extensions, kwargs...)

    # Instantiate the distribution PowerModels struct, without building the model.
    # Distribution and transmission structs share the same JuMP model. The `jump_model` parameter is used by _IM.InitializeInfrastructureModel().
    d_pm = _PM.instantiate_model(d_data, d_model_type, method->nothing; ref_extensions=d_ref_extensions, jump_model=t_pm.model, kwargs...)

    # Build the combined model.
    build_method(t_pm, d_pm)

    return t_pm, d_pm
end

## Functions that manipulate data structures

"""
    add_td_coupling_data!(d_data; sub_nw, qs_ratio_bound=0.48)

Add to distribution single-network data structures the data needed for T&D coupling.

- Add to distribution network data structure a coupling generator connected to the reference
  bus; if it already exists, overwrite all its parameters except cost.
- Define a bound on reactive power exchanged between transmission and distribution networks:
  it is computed as a fraction `qs_ratio_bound` of the rated power of the distribution
  network (based on the rated power of existing and candidate branches connected to the
  reference bus); default value is 0.48 as per “Network Code on Demand Connection” –
  Commission Regulation (EU) 2016/1388.
- Add to `dim_prop(d_data, :sub_nw, sub_nw)` the following entries:
  - `d_gen`: the id of the coupling generator of distribution network;
  - `qs_ratio_bound`: the aforementioned allowable fraction of the rated power.

Return `d_gen`.

This function is intended to be the last that edits distribution single-network data
structures: it should be called just before `make_multinetwork`.
"""
function add_td_coupling_data!(d_data::Dict{String,Any}; sub_nw::Int, qs_ratio_bound::Float64=0.48)

    # Get the reference bus id
    d_ref_buses = [b for (b,bus) in d_data["bus"] if bus["bus_type"] == 3]
    if length(d_ref_buses) != 1
        Memento.error(_LOGGER, "Distribution network data must have 1 ref bus, but $(length(d_ref_buses)) are present.")
    end
    d_ref_bus = parse(Int, first(d_ref_buses))

    # Get a list of the ids of generators connected to the reference bus
    d_ref_gens = [g for (g,gen) in d_data["gen"] if gen["gen_bus"] == d_ref_bus]
    if length(d_ref_gens) > 1
        Memento.error(_LOGGER, "Distribution network data must have 0 or 1 generators connected to ref bus, but $(length(d_ref_gens)) are present.")
    end

    # Define an upper bound on the rated apparent power based on the rated power of existing and candidate branches connected to the reference bus
    d_s_rate = (
        # `Float64[...]` is to prevent `ArgumentError: reducing over an empty collection is not allowed` in case at least one of the following collections is empty.
        # TODO: when dropping support for Julia < 1.6, remove `Float64[...]` and use `init=0.0` keyword argument of `sum`.
          sum(Float64[branch["rate_a"] for (b,branch) in d_data["branch"]    if (branch["f_bus"]==d_ref_bus || branch["t_bus"]==d_ref_bus) && branch["br_status"]==1]) # In t_bus here, the t stands for "to" (not for "transmission" as in the rest of the function)
        + sum(Float64[branch["rate_a"] for (b,branch) in d_data["ne_branch"] if (branch["f_bus"]==d_ref_bus || branch["t_bus"]==d_ref_bus) && branch["br_status"]==1]) # In t_bus here, the t stands for "to" (not for "transmission" as in the rest of the function)
    )

    # Add a coupling generator connected to d_ref_bus (or use existing one) to model the transmission network
    if isempty(d_ref_gens)
        d_gen_idx = length(d_data["gen"]) + 1 # Assumes that gens have contiguous indices starting from 1, as should be
        d_data["gen"]["$d_gen_idx"] = Dict{String,Any}(
            "gen_bus"      => d_ref_bus,
            "index"        => d_gen_idx,
            "dispatchable" => true,
            "model"        => 2, # Cost model (2 => polynomial cost)
            "ncost"        => 0, # Number of cost coefficients
            "cost"         => Any[]
        )
    else
        d_gen_idx = parse(Int, first(d_ref_gens))
    end

    # Set coupling generator parameters
    d_gen = d_data["gen"]["$d_gen_idx"]
    d_gen["mbase"]      = d_data["baseMVA"]
    d_gen["pmin"]       = -d_s_rate
    d_gen["pmax"]       =  d_s_rate
    d_gen["qmin"]       = -d_s_rate
    d_gen["qmax"]       =  d_s_rate
    d_gen["gen_status"] = 1

    # Store the T&D coupling parameters
    dim_prop(d_data, :sub_nw, sub_nw)["qs_ratio_bound"] = qs_ratio_bound
    dim_prop(d_data, :sub_nw, sub_nw)["d_gen"] = d_gen_idx
end

"""
    add_td_coupling_data!(t_data, d_data; t_bus, sub_nw, qs_ratio_bound=0.48)

Add to transmission and distribution single-network data structures the data needed for T&D coupling.

In addition to `add_td_coupling_data!(d_data; sub_nw, qs_ratio_bound)`, do the following:
- Add to transmission network data structure a coupling generator connected to `t_bus`, that
  is the bus to which the distribution network is to be connected.
- Add to `dim_prop(d_data, :sub_nw, sub_nw)` an entry `t_gen` for the id of the coupling
  generator of transmission network.

Return `t_gen`.

This function is intended to be the last that edits transmission and distribution
single-network data structures: it should be called just before `make_multinetwork`.
"""
function add_td_coupling_data!(t_data::Dict{String,Any}, d_data::Dict{String,Any}; t_bus::Int, sub_nw::Int, qs_ratio_bound::Float64=0.48)

    d_gen_idx = add_td_coupling_data!(d_data; sub_nw, qs_ratio_bound)

    # Check that t_bus exists
    if !haskey(t_data["bus"], "$t_bus")
        Memento.error(_LOGGER, "Bus $t_bus does not exist in transmission network data.")
    end

    # Add a coupling generator connected to t_bus, to model the distribution network
    t_s_rate = (d_data["baseMVA"]/t_data["baseMVA"]) * d_data["gen"]["$d_gen_idx"]["pmax"]
    t_gen_idx = length(t_data["gen"]) + 1 # Assumes that gens have contiguous indices starting from 1, as should be
    t_data["gen"]["$t_gen_idx"] = Dict{String,Any}(
        "gen_bus"      => t_bus,
        "index"        => t_gen_idx,
        "dispatchable" => true,
        "mbase"        => t_data["baseMVA"],
        "pmin"         => -t_s_rate,
        "pmax"         =>  t_s_rate,
        "qmin"         => -t_s_rate,
        "qmax"         =>  t_s_rate,
        "gen_status"   => 1,
        "model"        => 2, # Cost model (2 => polynomial cost)
        "ncost"        => 0, # Number of cost coefficients
        "cost"         => Any[]
    )

    # Add generator id to T&D coupling parameters
    dim_prop(d_data, :sub_nw, sub_nw)["t_gen"] = t_gen_idx
end

## Solution processors

"""
    sol_td_coupling!(pm, solution)

Add T&D coupling data to solution.

Report in `solution` the active and reactive power that distribution network `pm` exchanges with the
transmission network (positive if from transmission to distribution) in units of `baseMVA` of `pm`.
"""
function sol_td_coupling!(pm::_PM.AbstractBFModel, solution::Dict{String,Any})
    solution = _PM.get_pm_data(solution)
    if haskey(solution, "nw")
        nws_sol = solution["nw"]
    else
        nws_sol = Dict("0" => solution)
    end

    for (nw, nw_sol) in nws_sol
        n = parse(Int, nw)
        if !(haskey(dim_prop(pm), :sub_nw) && haskey(dim_prop(pm, n, :sub_nw), "d_gen"))
                Memento.error(_LOGGER, "T&D coupling data is missing from the model of nw $nw.")
        end
        d_gen = dim_prop(pm, n, :sub_nw, "d_gen")
        nw_sol["td_coupling"] = Dict{String,Any}()
        nw_sol["td_coupling"]["p"] = nw_sol["gen"]["$d_gen"]["pg"]
        nw_sol["td_coupling"]["q"] = nw_sol["gen"]["$d_gen"]["qg"]
    end
end

## Functions that group constraint templates, provided for convenience

"""
Connect each distribution nw to the corresponding transmission nw and apply coupling constraints.

The coupling constraint is applied to the two generators that each distribution nw indicates.
"""
function constraint_td_coupling(t_pm::_PM.AbstractPowerModel, d_pm::_PM.AbstractBFModel)
    t_nws = nw_ids(t_pm)

    for s in keys(dim_prop(d_pm, :sub_nw))
        d_nws = nw_ids(d_pm; sub_nw = s)
        for i in 1:length(t_nws)
            t_nw = t_nws[i]
            d_nw = d_nws[i]

            constraint_td_coupling_power_balance(t_pm, d_pm, t_nw, d_nw)
        end
    end
end

## Constraint templates

"""
State the power conservation between a distribution nw and the corresponding transmission nw.
"""
function constraint_td_coupling_power_balance(t_pm::_PM.AbstractPowerModel, d_pm::_PM.AbstractBFModel, t_nw::Int, d_nw::Int)
    sub_nw = dim_prop(d_pm, d_nw, :sub_nw)
    t_gen = sub_nw["t_gen"]
    d_gen = sub_nw["d_gen"]
    t_mbase = _PM.ref(t_pm, t_nw, :gen, t_gen, "mbase")
    d_mbase = _PM.ref(d_pm, d_nw, :gen, d_gen, "mbase")

    constraint_td_coupling_power_balance_active(t_pm, d_pm, t_nw, d_nw, t_gen, d_gen, t_mbase, d_mbase)
    constraint_td_coupling_power_balance_reactive(t_pm, d_pm, t_nw, d_nw, t_gen, d_gen, t_mbase, d_mbase)
end

"""
Apply bounds on reactive power exchange at the point of common coupling (PCC) of a distribution nw.
"""
function constraint_td_coupling_power_reactive_bounds(d_pm::_PM.AbstractBFModel; nw::Int=_PM.nw_id_default)
    sub_nw = dim_prop(d_pm, nw, :sub_nw)
    d_gen = sub_nw["d_gen"]
    qs_ratio_bound = sub_nw["qs_ratio_bound"] # Allowable fraction of rated apparent power

    constraint_td_coupling_power_reactive_bounds(d_pm, nw, d_gen, qs_ratio_bound)
end

## Constraint implementations

""
function constraint_td_coupling_power_balance_active(t_pm::_PM.AbstractPowerModel, d_pm::_PM.AbstractBFModel, t_nw::Int, d_nw::Int, t_gen::Int, d_gen::Int, t_mbase::Float64, d_mbase::Float64)
    t_p_in = _PM.var(t_pm, t_nw, :pg, t_gen)
    d_p_in = _PM.var(d_pm, d_nw, :pg, d_gen)
    JuMP.@constraint(t_pm.model, t_mbase*t_p_in + d_mbase*d_p_in == 0.0) # t_pm.model == d_pm.model
end

""
function constraint_td_coupling_power_balance_reactive(t_pm::_PM.AbstractPowerModel, d_pm::_PM.AbstractBFModel, t_nw::Int, d_nw::Int, t_gen::Int, d_gen::Int, t_mbase::Float64, d_mbase::Float64)
    t_q_in = _PM.var(t_pm, t_nw, :qg, t_gen)
    d_q_in = _PM.var(d_pm, d_nw, :qg, d_gen)
    JuMP.@constraint(t_pm.model, t_mbase*t_q_in + d_mbase*d_q_in == 0.0) # t_pm.model == d_pm.model
end

"Nothing to do because the transmission network model does not support reactive power."
function constraint_td_coupling_power_balance_reactive(t_pm::_PM.AbstractActivePowerModel, d_pm::_PM.AbstractBFModel, t_nw::Int, d_nw::Int, t_gen::Int, d_gen::Int, t_mbase::Float64, d_mbase::Float64)
end

""
function constraint_td_coupling_power_reactive_bounds(d_pm::_PM.AbstractBFModel, d_nw::Int, d_gen::Int, qs_ratio_bound)

    # Compute the rated apparent power of the distribution network, based on the rated power of
    # existing and candidate branches connected to its reference bus. This value depends on the
    # indicator variables of both existing (if applicable) and candidate branches (i.e. whether they
    # are built or not).
    if haskey(_PM.var(d_pm, d_nw), :z_branch) # Some `branch`es can be replaced by `ne_branch`es
        z_branch = _PM.var(d_pm, d_nw, :z_branch)
        s_rate = (
            sum(branch["rate_a"] * get(z_branch, b, 1.0) for (b,branch) in _PM.ref(d_pm, d_nw, :frb_branch))
            + sum(branch["rate_a"] * _PM.var(d_pm, d_nw, :branch_ne, b) for (b,branch) in _PM.ref(d_pm, d_nw, :frb_ne_branch))
        )
    else # No `ne_branch`es at all
        s_rate = sum(branch["rate_a"] for (b,branch) in _PM.ref(d_pm, d_nw, :frb_branch))
    end

    q = _PM.var(d_pm, d_nw, :qg, d_gen) # Exchanged reactive power (positive if from T to D)

    JuMP.@constraint(d_pm.model, q <=  qs_ratio_bound*s_rate)
    JuMP.@constraint(d_pm.model, q >= -qs_ratio_bound*s_rate)
end


text:
using BAT
using Test

import Cuba

using Base.Threads
using Distributions, ValueShapes

@testset "cuba_integration" begin
    function test_integration(target, algorithm::BAT.CubaIntegration)
        r = bat_integrate(target, algorithm).result
        val_expected = 1
        @test isapprox(r.val, val_expected, rtol = 100 * algorithm.rtol)
        # @test r.err < 100 * abs(r.val - val_expected)
    end

    uvprior = truncated(Normal(), -1, 3)

    mvprior = NamedTupleDist(
        a = Beta(0.5, 0.5),
        b = truncated(Normal(), -1, 3),
        c = Beta(5, 1),
        e = Uniform(-4, 7)
    )

    mvprior_simple = NamedTupleDist(
        a = truncated(Normal(), -1, 3),
        b = truncated(Normal(), 1, 3),
    )

    test_integration(uvprior, VEGASIntegration(trafo = NoDensityTransform()))
    test_integration(uvprior, VEGASIntegration(trafo = NoDensityTransform(), log_density_shift = 10))
    test_integration(mvprior, VEGASIntegration(trafo = NoDensityTransform(), nthreads = 1))
    test_integration(mvprior, VEGASIntegration(trafo = NoDensityTransform(), nthreads = nthreads()))
    test_integration(mvprior, VEGASIntegration())

    test_integration(mvprior, SuaveIntegration(trafo = NoDensityTransform()))

    test_integration(mvprior_simple, DivonneIntegration(trafo = NoDensityTransform(), rtol = 1e-4))

    test_integration(mvprior, CuhreIntegration(trafo = NoDensityTransform()))
end


text:
using Statistics, Econometrics
function Analyze(chain)
    posmean = vec(mean(chain,dims=1))
    posmedian = vec(median(chain,dims=1))
    inci = zeros(3)
    lower = zeros(3)
    upper = zeros(3)
    for i = 1:3
        lower[i] = quantile(chain[:,i],0.05)
        upper[i] = quantile(chain[:,i],0.95)
        inci[i] = θtrue[i] >= lower[i] && θtrue[i] <= upper[i]
    end
    return vcat(posmean[:], posmedian[:], lower[:], upper[:], inci[:])
end


text:
"""
    GibbsSampling(;ϵ::T=1e-5, nBurnin::Int=100, thinning::Int=1)

Draw samples from the true posterior via Gibbs Sampling.

## Keywords arguments
- `ϵ::T` : convergence criteria
- `nBurnin::Int` : Number of samples discarded before starting to save samples
- `thinning::Int` : Frequency at which samples are saved
"""
mutable struct GibbsSampling{T<:Real,N,Tx,Ty} <: SamplingInference{T}
    nBurnin::Int # Number of burnin samples
    thinning::Int # Frequency at which samples are saved
    ϵ::T # Convergence criteria
    nIter::Integer # Number of samples computed
    nSamples::Int # Number of data samples
    HyperParametersUpdated::Bool # Flag for updating kernel matrices
    opt::NTuple{N,SOptimizer}
    xview::Tx
    yview::Ty
    sample_store::Vector{Vector{Vector{T}}}
    function GibbsSampling{T}(nBurnin::Int, thinning::Int, ϵ::Real) where {T}
        nBurnin >= 0 || error("nBurnin should be positive")
        thinning >= 0 || error("thinning should be positive")
        return new{T,1,Vector{T},Vector{T}}(nBurnin, thinning, ϵ)
    end
    function GibbsSampling{T}(
        nBurnin::Int,
        thinning::Int,
        ϵ::Real,
        ::Vector{<:Int},
        nSamples::Int,
        nLatent::Int,
        xview::Tx,
        yview::Ty,
    ) where {T,Tx,Ty}
        opts = ntuple(_ -> SOptimizer{T}(nothing), nLatent)
        return new{T,nLatent,Tx,Ty}(
            nBurnin, thinning, ϵ, 0, nSamples, true, opts, xview, yview, []
        )
    end
end

isStochastic(::GibbsSampling) = false
getρ(::GibbsSampling{T}) where {T} = one(T)
nMinibatch(i::GibbsSampling) = i.nSamples

function GibbsSampling(; ϵ::T=1e-5, nBurnin::Int=100, thinning::Int=1) where {T<:Real}
    return GibbsSampling{T}(nBurnin, thinning, ϵ)
end

function Base.show(io::IO, ::GibbsSampling{T}) where {T<:Real}
    return print(io, "Gibbs Sampler")
end

function tuple_inference(
    inf::GibbsSampling{T},
    nLatent::Int,
    nFeatures::Vector{<:Int},
    nSamples::Int,
    ::Int,
    xview,
    yview,
) where {T}
    return GibbsSampling{T}(
        inf.nBurnin, inf.thinning, inf.ϵ, nFeatures, nSamples, nLatent, xview, yview
    )
end

function sample_local!(l::AbstractLikelihood, y, f::Tuple{<:AbstractVector{T}}) where {T}
    return sample_local!(l, y, first(f))
end
set_ω!(l::AbstractLikelihood, ω) = l.θ .= ω
get_ω(l::AbstractLikelihood) = l.θ

function sample_global!(
    ∇E_μ::AbstractVector, ∇E_Σ::AbstractVector, X::AbstractVector, gp::SampledLatent{T}
) where {T}
    gp.post.Σ .= inv(Symmetric(2.0 * Diagonal(∇E_Σ) + inv(pr_cov(gp))))
    rand!(MvNormal(cov(gp) * (∇E_μ + pr_cov(gp) \ pr_mean(gp, X)), cov(gp)), gp.post.f)
    return deepcopy(posterior(gp).f)
end


text:
# License for this file: MIT (expat)
# Copyright 2020, DLR Institute of System Dynamics and Control

"""
    module ModiaMathTest

Evaluate the basic structure of the ModiaMath simulator supporting DifferentialEquations.jl
"""
module ModiaMathTest

import DifferentialEquations
export SimulationState, solve!

#---------------------------------------- ModiaMath definitions

#=
Used Types

- StateType: Type of the state vector, of its derivative and of time:
                x::Vector{StateType}
                der_x::Vector{StateType}
				time::StateType
             and all variables that are used to compute der_x from x and time.
             Parameters can be either of StateType (e.g. to support parameters with uncertainty)
             or RealType (if not to be included in sensitivity analysis)

- RealType : Type of floating point variables that are not of type StateType, such as:
             tolerance, startTime, StopTime, interval, nominal, ...
=#

mutable struct SimulationState{StateType,RealType}
    modelName::AbstractString

    x_start::Vector{StateType}
    der_x_start::Vector{StateType}

    # Cache
    der_x::Vector{StateType}

    # Type of model
    dae::Bool    # = true: DAE; =false: ODE

    # Functions for ODE models
    setStates!::Function
    getDerivatives!::Function

    # Functions for DAE models
    setStatesAndDerivatives!::Function

    # Functions for ODE and DAE models
    getResiduals!::Function
    evaluate!::Function
    outputs!::Function
    Result                        # constructor for result data structure
    addToResult!::Function

    # Problem variables
    tolerance::RealType           # Relative integration tolerance
    startTime::RealType           # Start time of the simulation
    stopTime::RealType            # Stop time of the simulation
    interval::RealType            # Interval of the simulation

    # Time varying variables
    time::StateType
    storeResult::Bool

    # Simulation result
    result::Any                   # Instance of result

    SimulationState{StateType,RealType}(modelName, x_start,
                                        setStates!, getDerivatives!, setStatesAndDerivatives!, getResiduals!,
                                        evaluate!,  outputs!, Result, addToResult!;
                                        tolerance=1e-4, startTime=0.0, stopTime=1.0, interval=NaN, dae=true) where {StateType, RealType} =
         new(modelName, x_start, similar(x_start), similar(x_start), dae, setStates!, getDerivatives!, setStatesAndDerivatives!, getResiduals!,
             evaluate!, outputs!, Result, addToResult!, tolerance, startTime, stopTime, interval, RealType(0), false)
end

#
# x, der_x, residuals need not to be Vector{..}, e.g. Sundials.IDA provides N_Vector
#

"derivatives!: Called by ODE integrator for ODE model"
function derivatives!(der_x, x, simulationModel, t)::Nothing
    sim::SimulationState = simulationModel.simulationState
    sim.time = t
    sim.setStates!(simulationModel, x)
    sim.evaluate!(simulationModel)
    sim.getDerivatives!(simulationModel, der_x)
    return nothing
end

"ODEresiduals!: Called by DAE integrator for ODE model"
function ODEresiduals!(residuals, der_x, x, simulationModel, t)::Nothing
    sim::SimulationState = simulationModel.simulationState
    sim.time = t
    derivatives!(sim.der_x, x, simulationModel, t)
    for i = 1:length(der_x)
        residuals[i] = sim.der_x[i] - der_x[i]
    end
    return nothing
end

"ODEoutputs!: Called by ODE or DAE integrator for ODE model at every communication point"
function ODEoutputs!(x, t, integrator)::Nothing
    simulationModel = integrator.p
    sim::SimulationState = simulationModel.simulationState
    sim.time = t
    sim.setStates!(simulationModel, x)
    sim.storeResult = true
    sim.evaluate!(simulationModel)
    sim.storeResult = false
    sim.addToResult!(sim.result, t, simulationModel)
    return nothing
end

"initialize!: Called for DAE model to initialize"
function initialize!(simulationModel, sim::SimulationState{StateType,RealType})::Nothing where {StateType,RealType}
    if sim.dae
		error("... DAE initialization not yet implemented")
	else
		println("    ODE initialization function called")
		derivatives!(sim.der_x_start, sim.x_start, simulationModel, StateType(sim.startTime))
	end
    return nothing
end

function setStatesAndDerivativesDummy!(simulationModel, x, der_x)::Nothing
    error("setStatesAndDerivativesDummy! should never be called")
    return nothing
end

"SimulationState: Constructor for ODE Simulation State"
SimulationState{StateType,RealType}(modelName, x_start, setStates!, getDerivatives!, evaluate!, Result, addToResult!;
                                    tolerance=1e-4, startTime=0.0, stopTime=1.0, interval=NaN) where {StateType, RealType} =
    SimulationState{StateType,RealType}(modelName, x_start, setStates!, getDerivatives!, setStatesAndDerivativesDummy!, ODEresiduals!,
                                        evaluate!,  ODEoutputs!, Result, addToResult!; dae=false,
                                        tolerance=tolerance, startTime=startTime, stopTime=stopTime, interval=interval)

"solve!: Return DifferentialEquations.solve(...)"
function solve!(simulationModel, integrator; tolerance=NaN, startTime=NaN, stopTime=NaN, interval=NaN)
    sim::SimulationState = simulationModel.simulationState
    RealType  = typeof(sim.tolerance)
    StateType = typeof(sim.time)
    println("... Simulation started for model ", sim.modelName, " using integrator ", typeof(integrator), " with StateType = ", string(StateType) )

    # Construct a new result instance
    sim.result = sim.Result(simulationModel)

    # Change tolerance, startTime, stopTime, interval if required
    if !isnan(tolerance); @assert(tolerance > 0); sim.tolerance = RealType(tolerance); end
    if !isnan(startTime); sim.startTime = RealType(startTime); end
    if !isnan(stopTime) ; sim.stopTime  = RealType(stopTime) ; end
    if isnan(interval)
		if isnan(sim.interval)
			sim.interval = abs(sim.stopTime - sim.startTime)/500
		end
	else
		@assert(interval > 0)
		sim.interval = RealType(interval)
    end

    # Initialize model
    sim.time = StateType(sim.startTime)
    initialize!(simulationModel, sim)

    # Define problem based on integrator and model type
    if typeof(integrator) <: DifferentialEquations.IDA
        # DAE integrator
        if !sim.dae
            # ODE model
            problem = DifferentialEquations.DAEProblem(ODEresiduals!, sim.der_x_start, sim.x_start, (StateType(sim.startTime), sim.stopTime), simulationModel)
        else
            # DAE model
            error("DAE model not yet implemented")
        end
    else
        # Hope that it is an ODE integrator
        problem = DifferentialEquations.ODEProblem(derivatives!, sim.x_start, (StateType(sim.startTime), sim.stopTime), simulationModel)
    end

    # Callbacks
    tspan = sim.startTime:sim.interval:sim.stopTime
    cb1 = DifferentialEquations.FunctionCallingCallback(sim.outputs!, funcat=tspan)

    # Simulation
    solution = DifferentialEquations.solve(problem, integrator, reltol=sim.tolerance, saveat=sim.stopTime-sim.startTime, callback=cb1)
end

end


text:
@testset "Distribution Tests" begin
  normal = Normal(0, 1)
  @test JSON.json(normal) == "{\"mu\":0.0,\"sigma\":1.0,\"_type\":\"IncrementalInference.PackedNormal\"}"

  rayleigh = Rayleigh(1)
  @test JSON.json(rayleigh) == "{\"sigma\":1.0,\"_type\":\"IncrementalInference.PackedRayleigh\"}"

  fullnormal = FullNormal([1.0,2.0,3.0], [1.0 2.0 3.0; 4.0 5.0 6.0])
  @test JSON.json(fullnormal) == "{\"mu\":[1.0,2.0,3.0],\"cov\":[1.0,4.0,2.0,5.0,3.0,6.0],\"_type\":\"IncrementalInference.PackedFullNormal\"}"

  uniform = Uniform(0, 1)
  @test JSON.json(uniform) == "{\"a\":0.0,\"b\":1.0,\"_type\":\"IncrementalInference.PackedUniform\"}"

  categorical = Categorical([0.4, 0.6])
  @test JSON.json(categorical) == "{\"p\":[0.4,0.6],\"_type\":\"IncrementalInference.PackedCategorical\"}"
end


text:
module GPLikelihoods

using Distributions
using Functors
using InverseFunctions: InverseFunctions
using LinearAlgebra
using Random
using StatsFuns

export BernoulliLikelihood,
    CategoricalLikelihood,
    GaussianLikelihood,
    HeteroscedasticGaussianLikelihood,
    PoissonLikelihood,
    ExponentialLikelihood,
    GammaLikelihood,
    NegativeBinomialLikelihood
export NBParamI, NBParamII, NBParamPower, NBParamSuccess, NBParamFailure
export Link,
    ChainLink,
    BijectiveSimplexLink,
    ExpLink,
    LogLink,
    InvLink,
    SqrtLink,
    SquareLink,
    LogitLink,
    LogisticLink,
    ProbitLink,
    NormalCDFLink,
    SoftMaxLink
export expected_loglikelihood

# Links
include("links.jl")

# Likelihoods
abstract type AbstractLikelihood end

include("expectations.jl")
include("likelihoods/bernoulli.jl")
include("likelihoods/categorical.jl")
include("likelihoods/gaussian.jl")
include("likelihoods/poisson.jl")
include("likelihoods/gamma.jl")
include("likelihoods/exponential.jl")
include("likelihoods/negativebinomial.jl")

# TestInterface module
include("TestInterface.jl")

end # module


text:
struct GenericTimeGF{T, scalar, U <: AbstractTimeGrid} <: AbstractTimeGF{T, scalar}
  grid::U
  data::GenericStorage{T, scalar}
end

norbitals(G::GenericTimeGF) = G.data.norb

function GenericTimeGF(::Type{T}, grid::AbstractTimeGrid, norb=1, scalar=false) where T <: Number
  N = length(grid)
  data = GenericStorage(T, N, N, norb, scalar)
  GenericTimeGF(grid, data)
end
GenericTimeGF(grid::AbstractTimeGrid, norb=1, scalar=false) = GenericTimeGF(ComplexF64, grid, norb, scalar)

function GenericTimeGF(f::Function, ::Type{T}, grid::AbstractTimeGrid, norb=1, scalar=false) where T <: Number
  N = length(grid)
  G = GenericTimeGF(T, grid, norb, scalar)

  for t1 in grid
    for t2 in grid
      G[t1, t2] = f(t1, t2)
    end
  end
  return G
end
GenericTimeGF(f::Function, grid::AbstractTimeGrid, norb=1, scalar=false) = GenericTimeGF(f, ComplexF64, grid, norb, scalar)

@inline function Base.getindex(G::GenericTimeGF, k, l, t1::TimeGridPoint, t2::TimeGridPoint, greater=true)
  val = G.data[k, l, t1.cidx, t2.cidx]
  (!greater && t1.cidx == t2.cidx) && (val += jump(G))
  return val
end

function Base.setindex!(G::GenericTimeGF, v, k, l, t1::TimeGridPoint, t2::TimeGridPoint)
  G.data[k, l, t1.cidx, t2.cidx] = v
end

function jump(G::GenericTimeGF)
  t0_plus = branch_bounds(G.grid, forward_branch)[1]
  t0_minus = branch_bounds(G.grid, backward_branch)[2]
  return G[t0_plus, t0_minus] - G[t0_plus, t0_plus]
end

function Base.similar(G::T) where T <: GenericTimeGF
  T(G.grid, similar(G.data))
end

function Base.zero(G::T) where T <: GenericTimeGF
  T(G.grid, zero(G.data))
end

function GenericTimeGF(dos::AbstractDOS, grid::FullTimeGrid)
  G = TimeInvariantFullTimeGF(dos, grid)
  GenericTimeGF(grid,1,true) do t1, t2
    G[t1,t2]
  end
end

function GenericTimeGF(dos::AbstractDOS, grid::ImaginaryTimeGrid)
  G = ImaginaryTimeGF(dos, grid)
  GenericTimeGF(grid,1,true) do t1, t2
    G[t1,t2]
  end
end

function GenericTimeGF(dos::AbstractDOS, β, grid::KeldyshTimeGrid)
  G = TimeInvariantKeldyshTimeGF(dos, β, grid)
  GenericTimeGF(grid,1,true) do t1, t2
    G[t1,t2]
  end
end


text:
using Cameras
using PylonCameras

const images_to_grab = 30

function acquire_images()
    camera = PylonCamera(
        max_num_buffer = 8,
        grab_result_wait_timeout_ms = 100,
        grab_result_retrieve_timeout_ms = 1)
    try
        println("Using $(PylonCameras.info(camera))")
        open!(camera)
        start!(camera, images_to_grab)

        while isrunning(camera)
            img = take!(camera)
            print("ID $(id(img)), No. $(image_number(img)) @ $(timestamp(img)), size $(size(img)) : ")
            @show img[1,1,1]
            @time sum(img)
            release!(img)
        end
    catch e
        println(e)
    finally
        stop!(camera)
        close!(camera)
    end
end

acquire_images()


text:
using Distributions
using CSVFiles
using DataFrames

function getsim()
        mcs = @defsim begin

                ############################################################################
                # Define random variables (RVs)
                ############################################################################

                #The folllowing RVs are in more than one component.  For clarity they are
                #set here as opposed to below within the blocks of RVs separated by component
                #so that they are not set more than once.

                save_savingsrate = TriangularDist(10, 20, 15) # components: MarketDamages, NonMarketDamages, GDP, SLRDamages
                tcal_CalibrationTemp = TriangularDist(2.5, 3.5, 3.) # components: MarketDamages, NonMarketDamages

                wincf_weightsfactor["USA"] = TriangularDist(.6, 1, .8) # components: MarketDamages, NonMarketDamages, SLRDamages, Discountinuity
                wincf_weightsfactor["OECD"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor["USSR"] = TriangularDist(.2, .6, .4)
                wincf_weightsfactor["China"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor["SEAsia"] = TriangularDist(.4, 1.2, .8)
                wincf_weightsfactor["Africa"] = TriangularDist(.4, .8, .6)
                wincf_weightsfactor["LatAmerica"] = TriangularDist(.4, .8, .6)

                automult_autonomouschange = TriangularDist(0.5, 0.8, 0.65)  #components: AdaptationCosts, AbatementCosts

                #The following RVs are divided into blocks by component

                # CO2cycle
                air_CO2fractioninatm = TriangularDist(57, 67, 62)
                res_CO2atmlifetime = TriangularDist(50, 100, 70)
                ccf_CO2feedback = TriangularDist(4, 15, 10)
                ccfmax_maxCO2feedback = TriangularDist(30, 80, 50)
                stay_fractionCO2emissionsinatm = TriangularDist(0.25,0.35,0.3)

                # SulphateForcing
                d_sulphateforcingbase = TriangularDist(-0.8, -0.2, -0.4)
                ind_slopeSEforcing_indirect = TriangularDist(-0.8, 0, -0.4)

                # ClimateTemperature
                rlo_ratiolandocean = TriangularDist(1.2, 1.6, 1.4)
                pole_polardifference = TriangularDist(1, 2, 1.5)
                frt_warminghalflife = TriangularDist(10, 65, 30)
                tcr_transientresponse = TriangularDist(1, 2.8, 1.3)

                # SeaLevelRise
                s0_initialSL = TriangularDist(0.1, 0.2, 0.15)
                sltemp_SLtemprise = TriangularDist(0.7, 3., 1.5)
                sla_SLbaselinerise = TriangularDist(0.5, 1.5, 1.)
                sltau_SLresponsetime = TriangularDist(500, 1500, 1000)

                # GDP
                isat0_initialimpactfxnsaturation = TriangularDist(20, 50, 30)

                # MarketDamages
                iben_MarketInitialBenefit = TriangularDist(0, .3, .1)
                W_MarketImpactsatCalibrationTemp = TriangularDist(.2, .8, .5)
                pow_MarketImpactExponent = TriangularDist(1.5, 3, 2)
                ipow_MarketIncomeFxnExponent = TriangularDist(-.3, 0, -.1)

                # NonMarketDamages
                iben_NonMarketInitialBenefit = TriangularDist(0, .2, .05)
                w_NonImpactsatCalibrationTemp = TriangularDist(.1, 1, .5)
                pow_NonMarketExponent = TriangularDist(1.5, 3, 2)
                ipow_NonMarketIncomeFxnExponent = TriangularDist(-.2, .2, 0)

                # SLRDamages
                scal_calibrationSLR = TriangularDist(0.45, 0.55, .5)
                #iben_SLRInitialBenefit = TriangularDist(0, 0, 0) # only usable if lb <> ub
                W_SatCalibrationSLR = TriangularDist(.5, 1.5, 1)
                pow_SLRImpactFxnExponent = TriangularDist(.5, 1, .7)
                ipow_SLRIncomeFxnExponent = TriangularDist(-.4, -.2, -.3)

                # Discountinuity
                rand_discontinuity = Uniform(0, 1)
                tdis_tolerabilitydisc = TriangularDist(2, 4, 3)
                pdis_probability = TriangularDist(10, 30, 20)
                wdis_gdplostdisc = TriangularDist(5, 25, 15)
                ipow_incomeexponent = TriangularDist(-.3, 0, -.1)
                distau_discontinuityexponent = TriangularDist(20, 200, 50)

                # EquityWeighting
                civvalue_civilizationvalue = TriangularDist(1e10, 1e11, 5e10)
                ptp_timepreference = TriangularDist(0.1,2,1)
                emuc_utilityconvexity = TriangularDist(0.5,2,1)

                # AbatementCosts
                AbatementCostParametersCO2_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,75,0)
                AbatementCostParametersCH4_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-25,100,0)
                AbatementCostParametersN2O_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,50,0)
                AbatementCostParametersLin_emit_UncertaintyinBAUEmissFactorinFocusRegioninFinalYear = TriangularDist(-50,50,0)

                AbatementCostParametersCO2_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,40,20)
                AbatementCostParametersCH4_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
                AbatementCostParametersN2O_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)
                AbatementCostParametersLin_q0propinit_CutbacksinNegativeCostinFocusRegioninBaseYear = TriangularDist(0,20,10)

                AbatementCostParametersCO2_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)
                AbatementCostParametersCH4_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-8000,-1000,-4000)
                AbatementCostParametersN2O_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-15000,0,-7000)
                AbatementCostParametersLin_c0init_MostNegativeCostCutbackinBaseYear = TriangularDist(-400,-100,-200)

                AbatementCostParametersCO2_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)
                AbatementCostParametersCH4_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
                AbatementCostParametersN2O_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(35,70,50)
                AbatementCostParametersLin_qmaxminusq0propinit_MaxCutbackCostatPositiveCostinBaseYear = TriangularDist(60,80,70)

                AbatementCostParametersCO2_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,700,400)
                AbatementCostParametersCH4_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(3000,10000,6000)
                AbatementCostParametersN2O_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(2000,60000,20000)
                AbatementCostParametersLin_cmaxinit_MaximumCutbackCostinFocusRegioninBaseYear = TriangularDist(100,600,300)

                AbatementCostParametersCO2_ies_InitialExperienceStockofCutbacks = TriangularDist(100000,200000,150000)
                AbatementCostParametersCH4_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)
                AbatementCostParametersN2O_ies_InitialExperienceStockofCutbacks = TriangularDist(30,80,50)
                AbatementCostParametersLin_ies_InitialExperienceStockofCutbacks = TriangularDist(1500,2500,2000)

                #the following variables need to be set, but set the same in all 4 abatement cost components
                #note that for these regional variables, the first region is the focus region (EU), which is set in the preceding code, and so is always one for these variables

                emitf_uncertaintyinBAUemissfactor["USA"] = TriangularDist(0.8,1.2,1.0)
                emitf_uncertaintyinBAUemissfactor["OECD"] = TriangularDist(0.8,1.2,1.0)
                emitf_uncertaintyinBAUemissfactor["USSR"] = TriangularDist(0.65,1.35,1.0)
                emitf_uncertaintyinBAUemissfactor["China"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["SEAsia"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["Africa"] = TriangularDist(0.5,1.5,1.0)
                emitf_uncertaintyinBAUemissfactor["LatAmerica"] = TriangularDist(0.5,1.5,1.0)

                q0f_negativecostpercentagefactor["USA"] = TriangularDist(0.75,1.5,1.0)
                q0f_negativecostpercentagefactor["OECD"] = TriangularDist(0.75,1.25,1.0)
                q0f_negativecostpercentagefactor["USSR"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["China"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["SEAsia"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["Africa"] = TriangularDist(0.4,1.0,0.7)
                q0f_negativecostpercentagefactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)

                cmaxf_maxcostfactor["USA"] = TriangularDist(0.8,1.2,1.0)
                cmaxf_maxcostfactor["OECD"] = TriangularDist(1.0,1.5,1.2)
                cmaxf_maxcostfactor["USSR"] = TriangularDist(0.4,1.0,0.7)
                cmaxf_maxcostfactor["China"] = TriangularDist(0.8,1.2,1.0)
                cmaxf_maxcostfactor["SEAsia"] = TriangularDist(1,1.5,1.2)
                cmaxf_maxcostfactor["Africa"] = TriangularDist(1,1.5,1.2)
                cmaxf_maxcostfactor["LatAmerica"] = TriangularDist(0.4,1.0,0.7)

                q0propmult_cutbacksatnegativecostinfinalyear = TriangularDist(0.3,1.2,0.7)
                qmax_minus_q0propmult_maxcutbacksatpositivecostinfinalyear = TriangularDist(1,1.5,1.3)
                c0mult_mostnegativecostinfinalyear = TriangularDist(0.5,1.2,0.8)
                curve_below_curvatureofMACcurvebelowzerocost = TriangularDist(0.25,0.8,0.45)
                curve_above_curvatureofMACcurveabovezerocost = TriangularDist(0.1,0.7,0.4)
                cross_experiencecrossoverratio = TriangularDist(0.1,0.3,0.2)
                learn_learningrate = TriangularDist(0.05,0.35,0.2)

                # AdaptationCosts
                AdaptiveCostsSeaLevel_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
                AdaptiveCostsSeaLevel_ci_costimpact_eu = TriangularDist(0.0005, 0.002, 0.001)
                AdaptiveCostsEconomic_cp_costplateau_eu = TriangularDist(0.005, 0.02, 0.01)
                AdaptiveCostsEconomic_ci_costimpact_eu = TriangularDist(0.001, 0.008, 0.003)
                AdaptiveCostsNonEconomic_cp_costplateau_eu = TriangularDist(0.01, 0.04, 0.02)
                AdaptiveCostsNonEconomic_ci_costimpact_eu = TriangularDist(0.002, 0.01, 0.005)

                cf_costregional["USA"] = TriangularDist(0.6, 1, 0.8)
                cf_costregional["OECD"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["USSR"] = TriangularDist(0.2, 0.6, 0.4)
                cf_costregional["China"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["SEAsia"] = TriangularDist(0.4, 1.2, 0.8)
                cf_costregional["Africa"] = TriangularDist(0.4, 0.8, 0.6)
                cf_costregional["LatAmerica"] = TriangularDist(0.4, 0.8, 0.6)

                ############################################################################
                # Indicate which parameters to save for each model run
                ############################################################################

                save(EquityWeighting.td_totaldiscountedimpacts,
                        EquityWeighting.tpc_totalaggregatedcosts,
                        EquityWeighting.tac_totaladaptationcosts,
                        EquityWeighting.te_totaleffect,
                        co2cycle.c_CO2concentration,
                        TotalForcing.ft_totalforcing,
                        ClimateTemperature.rt_g_globaltemperature,
                        SeaLevelRise.s_sealevel,
                        SLRDamages.rgdp_per_cap_SLRRemainGDP,
                        MarketDamages.rgdp_per_cap_MarketRemainGDP,
                        NonMarketDamages.rgdp_per_cap_NonMarketRemainGDP,
                        Discontinuity.rgdp_per_cap_NonMarketRemainGDP)

                end #def
        return mcs
end

#Reformat the RV results into the format used for testing
function reformat_RV_outputs(samplesize::Int; output_path::String = joinpath(@__DIR__, "../output"))

    #create vectors to hold results of Monte Carlo runs
    td=zeros(samplesize);
    tpc=zeros(samplesize);
    tac=zeros(samplesize);
    te=zeros(samplesize);
    ft=zeros(samplesize);
    rt_g=zeros(samplesize);
    s=zeros(samplesize);
    c_co2concentration=zeros(samplesize);
    rgdppercap_slr=zeros(samplesize);
    rgdppercap_market=zeros(samplesize);
    rgdppercap_nonmarket=zeros(samplesize);
    rgdppercap_disc=zeros(samplesize);

    #load raw data
    #no filter
    td      = load_RV("EquityWeighting_td_totaldiscountedimpacts", "td_totaldiscountedimpacts"; output_path = output_path)
    tpc     = load_RV("EquityWeighting_tpc_totalaggregatedcosts", "tpc_totalaggregatedcosts"; output_path = output_path)
    tac     = load_RV("EquityWeighting_tac_totaladaptationcosts", "tac_totaladaptationcosts"; output_path = output_path)
    te      = load_RV("EquityWeighting_te_totaleffect", "te_totaleffect"; output_path = output_path)

    #time index
    c_co2concentration = load_RV("co2cycle_c_CO2concentration", "c_CO2concentration"; output_path = output_path)
    ft      = load_RV("TotalForcing_ft_totalforcing", "ft_totalforcing"; output_path = output_path)
    rt_g    = load_RV("ClimateTemperature_rt_g_globaltemperature", "rt_g_globaltemperature"; output_path = output_path)
    s       = load_RV("SeaLevelRise_s_sealevel", "s_sealevel"; output_path = output_path)

    #region index
    rgdppercap_slr          = load_RV("SLRDamages_rgdp_per_cap_SLRRemainGDP", "rgdp_per_cap_SLRRemainGDP"; output_path = output_path)
    rgdppercap_market       = load_RV("MarketDamages_rgdp_per_cap_MarketRemainGDP", "rgdp_per_cap_MarketRemainGDP"; output_path = output_path)
    rgdppercap_nonmarket    =load_RV("NonMarketDamages_rgdp_per_cap_NonMarketRemainGDP", "rgdp_per_cap_NonMarketRemainGDP"; output_path = output_path)
    rgdppercap_disc         = load_RV("NonMarketDamages_rgdp_per_cap_NonMarketRemainGDP", "rgdp_per_cap_NonMarketRemainGDP"; output_path = output_path)

    #resave data
    df=DataFrame(td=td,tpc=tpc,tac=tac,te=te,c_co2concentration=c_co2concentration,ft=ft,rt_g=rt_g,sealevel=s,rgdppercap_slr=rgdppercap_slr,rgdppercap_market=rgdppercap_market,rgdppercap_nonmarket=rgdppercap_nonmarket,rgdppercap_di=rgdppercap_disc)
    save(joinpath(output_path, "mimipagemontecarlooutput.csv"),df)
end

function do_monte_carlo_runs(samplesize::Int; output_dir::String = joinpath(@__DIR__, "../output"))

        # get simulation
        mcs = getsim()

        # get a model
        m = get_model()
        run(m)

        # Run
        res = run(mcs, m, samplesize; trials_output_filename = joinpath(output_dir, "trialdata.csv"), results_output_dir = output_dir)

        # reformat outputs for testing and analysis
        MimiPAGE2009.reformat_RV_outputs(samplesize, output_path = output_dir)
end


text:
struct GeneralConstraint{T}
    evaluate::Any
    jacobian::Any
    hessian::Any
    num_variables::Int
    num_parameter::Int
    num_constraint::Int
    num_jacobian::Int
    num_hessian::Int
    jacobian_sparsity::Vector{Vector{Int}}
    hessian_sparsity::Vector{Vector{Int}}
    evaluate_cache::Vector{T}
    jacobian_cache::Vector{T}
    hessian_cache::Vector{T}
    indices_inequality::Vector{Int}
end

function GeneralConstraint(f::Function, num_variables::Int, num_parameter::Int;
    indices_inequality=collect(1:0),
    evaluate_hessian=false)

    #TODO: option to load/save methods
    @variables z[1:num_variables], w[1:num_parameter]
    evaluate = f(z, w)
    jacobian = Symbolics.sparsejacobian(evaluate, z)
    evaluate_func = eval(Symbolics.build_function(evaluate, z, w)[2])
    jacobian_func = eval(Symbolics.build_function(jacobian.nzval, z, w)[2])
    num_constraint = length(evaluate)
    num_jacobian = length(jacobian.nzval)
    jacobian_sparsity = [findnz(jacobian)[1:2]...]
    if evaluate_hessian
        @variables λ[1:num_constraint]
        lag_con = dot(λ, evaluate)
        hessian = Symbolics.sparsehessian(lag_con, z)
        hessian_func = eval(Symbolics.build_function(hessian.nzval, z, w, λ)[2])
        hessian_sparsity = [findnz(hessian)[1:2]...]
        num_hessian = length(hessian.nzval)
    else
        hessian_func = Expr(:null)
        hessian_sparsity = [Int[]]
        num_hessian = 0
    end

    return GeneralConstraint(
        evaluate_func,
        jacobian_func,
        hessian_func,
        num_variables,
        num_parameter,
        num_constraint,
        num_jacobian,
        num_hessian,
        jacobian_sparsity,
        hessian_sparsity,
        zeros(num_constraint),
        zeros(num_jacobian),
        zeros(num_hessian),
        indices_inequality)
end

function GeneralConstraint()
    return GeneralConstraint(
        (constraint, variables, parameters) -> nothing,
        (jacobian,   variables, parameters) -> nothing,
        (hessian,    variables, parameters) -> nothing,
        0, 0, 0, 0, 0,
        [Int[], Int[]],
        [Int[], Int[]],
        Float64[], Float64[], Float64[],
        collect(1:0))
end

function constraints!(violations, indices, general::GeneralConstraint{T}, variables, parameters) where T
    general.evaluate(general.evaluate_cache, variables, parameters)
    @views violations[indices] .= general.evaluate_cache
    fill!(general.evaluate_cache, 0.0) # TODO: confirm this is necessary
end

function jacobian!(jacobians, indices, general::GeneralConstraint{T}, variables, parameters) where T
    general.jacobian(general.jacobian_cache, variables, parameters)
    @views jacobians[indices] .= general.jacobian_cache
    fill!(general.jacobian_cache, 0.0) # TODO: confirm this is necessary
end

function hessian_lagrangian!(hessian, indices, general::GeneralConstraint{T}, variables, parameters, duals) where T
    if !isempty(general.hessian_cache)
        general.hessian(general.hessian_cache, variables, duals)
        @views hessian[indices] .+= general.hessian_cache
        fill!(general.hessian_cache, 0.0) # TODO: confirm this is necessary
    end
end

function sparsity_jacobian(general::GeneralConstraint{T}, num_variables::Int;
    row_shift=0) where T

    row = Int[]
    col = Int[]

    push!(row, (general.jacobian_sparsity[1] .+ row_shift)...)
    push!(col,  general.jacobian_sparsity[2]...)

    return collect(zip(row, col))
end

function sparsity_hessian(general::GeneralConstraint{T}, num_variables::Int) where T
    row = Int[]
    col = Int[]

    if !isempty(general.hessian_sparsity[1])
        shift = 0
        push!(row, (general.hessian_sparsity[1] .+ shift)...)
        push!(col, (general.hessian_sparsity[2] .+ shift)...)
    end

    return collect(zip(row, col))
end

num_constraint(general::GeneralConstraint{T}) where T = general.num_constraint
num_jacobian(general::GeneralConstraint{T}) where T = general.num_jacobian
# num_hessian(general::GeneralConstraint{T}) where T = general.num_hessian

constraint_indices(general::GeneralConstraint{T}; shift=0) where T = collect(shift .+ (1:general.num_constraint))

jacobian_indices(general::GeneralConstraint{T}; shift=0) where T = collect(shift .+ (1:general.num_jacobian))

function hessian_indices(general::GeneralConstraint{T}, key::Vector{Tuple{Int,Int}}, num_variables::Int) where T
    if !isempty(general.hessian_sparsity[1])
        row = Int[]
        col = Int[]
        shift = 0
        push!(row, (general.hessian_sparsity[1] .+ shift)...)
        push!(col, (general.hessian_sparsity[2] .+ shift)...)
        rc = collect(zip(row, col))
        indices = [findfirst(x -> x == i, key) for i in rc]
    else
        indices = Vector{Int}[]
    end
    return indices
end