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

using DelimitedFiles

@testset "trust_region_utils" begin
  @testset "solve_bound_constrained_trust_region unbounded" begin
    for solve_approximately in [true false]
      # min -x
      # || x || <= 5.0
      variable_lower_bounds = [-Inf]
      variable_upper_bounds = [Inf]
      current_point = [0.0]
      current_gradient = [-1.0]
      target_radius = 5.0
      norm_weights = [1.0]
      result = FirstOrderLp.solve_bound_constrained_trust_region(
        current_point,
        current_gradient,
        variable_lower_bounds,
        variable_upper_bounds,
        norm_weights,
        target_radius,
        solve_approximately,
      )
      @test result.value == -5.0
      @test result.solution == [5.0]

      # min x + y
      # 2 * x^2 + y^2 <= 6.0
      # [x*, y*] = [-1.0, -2.0]
      variable_lower_bounds = [-Inf, -Inf]
      variable_upper_bounds = [Inf, Inf]
      current_point = [0.0, 0.0]
      current_gradient = [1.0, 1.0]
      norm_weights = [2.0, 1.0]
      target_radius = sqrt(6.0)
      result = FirstOrderLp.solve_bound_constrained_trust_region(
        current_point,
        current_gradient,
        variable_lower_bounds,
        variable_upper_bounds,
        norm_weights,
        target_radius,
        solve_approximately,
      )
      @test result.solution ≈ [-1.0, -2.0] atol = 1.0e-8
      @test result.value ≈ -3.0 atol = 1.0e-8
    end
  end

  @testset "solve_bound_constrained_trust_region" begin
    # min -x
    # || x || <= 5.0
    variable_lower_bounds = [-Inf]
    variable_upper_bounds = [Inf]
    current_point = [0.0]
    current_gradient = [-1.0]
    target_radius = 5.0
    norm_weights = [1.0]
    result = FirstOrderLp.solve_bound_constrained_trust_region(
      current_point,
      current_gradient,
      variable_lower_bounds,
      variable_upper_bounds,
      norm_weights,
      target_radius,
      false,
    )
    @test result.solution == [5.0]

    # min -x
    # || x || <= 5.0
    # x <= 0.0
    variable_lower_bounds = [-Inf]
    variable_upper_bounds = [0.0]
    current_point = [0.0]
    current_gradient = [-1.0]
    result = FirstOrderLp.solve_bound_constrained_trust_region(
      current_point,
      current_gradient,
      variable_lower_bounds,
      variable_upper_bounds,
      norm_weights,
      target_radius,
      false,
    )
    @test result.solution == [0.0]

    # min -x
    # || x || <= 5.0
    # x <= 2.0
    variable_lower_bounds = [-Inf]
    variable_upper_bounds = [2.0]
    result = FirstOrderLp.solve_bound_constrained_trust_region(
      current_point,
      current_gradient,
      variable_lower_bounds,
      variable_upper_bounds,
      norm_weights,
      target_radius,
      false,
    )
    @test result.solution == [2.0]

    # min -2 * x - y
    # || [x; y] || <= 5.0
    # x <= 3.0
    # [x*, y*] = [3.0, 4.0]
    variable_lower_bounds = [-Inf, -Inf]
    variable_upper_bounds = [3.0, Inf]
    current_point = [0.0, 0.0]
    current_gradient = [-2.0, -1.0]
    norm_weights = [1.0, 1.0]
    result = FirstOrderLp.solve_bound_constrained_trust_region(
      current_point,
      current_gradient,
      variable_lower_bounds,
      variable_upper_bounds,
      norm_weights,
      target_radius,
      false,
    )
    @test result.solution ≈ [3.0, 4.0] atol = 1.0e-8

    # min -x
    # || [x; y] || <= 5.0
    # x <= 2.0
    # [x*, y*] = [2.0, 0.0]
    variable_lower_bounds = [-Inf, -Inf]
    variable_upper_bounds = [2.0, Inf]
    current_point = [0.0, 0.0]
    current_gradient = [-1.0, 0.0]
    result = FirstOrderLp.solve_bound_constrained_trust_region(
      current_point,
      current_gradient,
      variable_lower_bounds,
      variable_upper_bounds,
      norm_weights,
      target_radius,
      false,
    )
    @test result.solution == [2.0, 0.0]

    # Test norm_weights:
    # min -4.0 * x - 3.0 * y
    # 4.0^2 * x^2 + 3.0^2 * y^2 <= 2.0
    # [x*, y*] = [1.0 / 4.0, 1.0 / 3.0]
    # to see this substitute a = 4.0 x, b = 3.0 y, yielding:
    # min -a - b
    # a^2 + b^2 <= 2.0
    # [a*, b*] = [1.0, 1.0]
    variable_lower_bounds = [-Inf, -Inf]
    variable_upper_bounds = [Inf, Inf]
    current_point = [0.0, 0.0]
    current_gradient = [-4.0, -3.0]
    norm_weights = [4.0^2, 3.0^2]
    target_radius = sqrt(2.0)
    result = FirstOrderLp.solve_bound_constrained_trust_region(
      current_point,
      current_gradient,
      variable_lower_bounds,
      variable_upper_bounds,
      norm_weights,
      target_radius,
      false,
    )
    @test FirstOrderLp.weighted_norm(result.solution, norm_weights) ≈ sqrt(2.0) atol =
      1.0e-8
    @test result.solution ≈ [1.0 / 4.0, 1.0 / 3.0] atol = 1.0e-8

    # Hundred dimensional test
    # min -sum(x[i])
    # || x || <= r
    # x[i] <= i
    #
    # x*[i] = min(i, m)
    # Set r = || x* ||.
    n = 100
    variable_lower_bounds = zeros(n)
    variable_upper_bounds = 1.0 * collect(1:n)
    current_point = zeros(n)
    current_gradient = -ones(n)
    acceptable_relative_error = 1e-6
    norm_weights = ones(length(variable_lower_bounds))
    for m in [10.0, 50.0]
      target_radius = sqrt(sum([min(i, m)^2 for i in 1:n]))
      result = FirstOrderLp.solve_bound_constrained_trust_region(
        current_point,
        current_gradient,
        variable_lower_bounds,
        variable_upper_bounds,
        norm_weights,
        target_radius,
        false,
      )
      @test result.solution ≈ [min(i, m) for i in 1:n] atol = 1.0e-8
    end
  end

  @testset "bound_primal_and_dual_objective" begin
    lp = example_lp()
    primal_norm_weights = ones(4)
    dual_norm_weights = ones(3)
    # Test cases with zero or close to zero duality gap
    for norm in [FirstOrderLp.MAX_NORM, FirstOrderLp.EUCLIDEAN_NORM]
      primal_solution = [1.0, 0.0, 6.0, 2.0]
      dual_solution = [0.5, 4.0, 0.0]
      distance_to_optimality = 10.0
      result = FirstOrderLp.bound_optimal_objective(
        lp,
        primal_solution,
        dual_solution,
        primal_norm_weights,
        dual_norm_weights,
        distance_to_optimality,
        FirstOrderLp.MAX_NORM;
        solve_approximately = false,
      )

      @test result.lower_bound_value == -1.0
      @test result.upper_bound_value == -1.0

      primal_solution = [1.0, 0.0, 5.99999, 2.0]
      dual_solution = [0.50001, 4.0, 0.0]
      result = FirstOrderLp.bound_optimal_objective(
        lp,
        primal_solution,
        dual_solution,
        primal_norm_weights,
        dual_norm_weights,
        distance_to_optimality,
        FirstOrderLp.MAX_NORM;
        solve_approximately = false,
      )
      @test -1.01 < result.lower_bound_value < -1.0
      @test -1.0 < result.upper_bound_value < -0.99
    end

    primal_solution = [1.0, 0.0, 6.0, 1.0]
    dual_solution = [0.0, 4.0, 0.0]
    distance_to_optimality = 2.0
    result = FirstOrderLp.bound_optimal_objective(
      lp,
      primal_solution,
      dual_solution,
      primal_norm_weights,
      dual_norm_weights,
      distance_to_optimality,
      FirstOrderLp.MAX_NORM;
      solve_approximately = false,
    )
    @test result.lower_bound_value == -4.0
    @test result.upper_bound_value == 2.0
    @test result.lower_bound_value ==
          FirstOrderLp.corrected_dual_obj(lp, primal_solution, dual_solution)

    primal_solution = [3.0, 0.0, 6.0, 0.0]
    dual_solution = [0.0, 4.0, 0.0]
    distance_to_optimality = 5.0
    result = FirstOrderLp.bound_optimal_objective(
      lp,
      primal_solution,
      dual_solution,
      primal_norm_weights,
      dual_norm_weights,
      distance_to_optimality,
      FirstOrderLp.EUCLIDEAN_NORM;
      solve_approximately = false,
    )
    @test result.lower_bound_value == -4.0
    @test result.lagrangian_value == -1.0
    @test distance_to_optimality^2 ==
          norm(result.primal_solution - primal_solution, 2)^2 +
          norm(result.dual_solution - dual_solution, 2)^2

    @test result.upper_bound_value == 7.0

    primal_solution = [1.0, 1.0, 4.0, 1.0]
    dual_solution = [0.0, 0.0, 0.0]
    distance_to_optimality = 10.0
    result = FirstOrderLp.bound_optimal_objective(
      example_lp(),
      primal_solution,
      dual_solution,
      primal_norm_weights,
      dual_norm_weights,
      distance_to_optimality,
      FirstOrderLp.MAX_NORM;
      solve_approximately = false,
    )
    corrected_dual_obj =
      FirstOrderLp.corrected_dual_obj(lp, primal_solution, dual_solution)
    @test result.lower_bound_value ==
          FirstOrderLp.corrected_dual_obj(lp, primal_solution, dual_solution)

    primal_solution = [0.5, 0.5, 0.5, 1.0, 1.0, 1.0] # Solution in the interior.
    dual_solution = zeros(3)
    primal_norm_weights = ones(6)
    dual_norm_weights = ones(3)

    result = FirstOrderLp.bound_optimal_objective(
      example_cc_star_lp(),
      primal_solution,
      dual_solution,
      primal_norm_weights,
      dual_norm_weights,
      distance_to_optimality,
      FirstOrderLp.MAX_NORM;
      solve_approximately = false,
    )
    @test result.lagrangian_value == result.upper_bound_value
    @test result.lower_bound_value < result.lagrangian_value
  end
end


text:
"""
htc_fs:
TODO: we will use a consant heat transfer coefficient here, which is
not ideal for the future, because the heat transfer coefficient depends
on Re and Pr and therefore is dependet of flow speed (turbulence) and
the thermal properties near the "wall".
Therefore, we will have to expand this function to be dependend of
(Tf,Tg,material_solid, material_fluid) in the future ..
"""
function htc_fs(ht_type::ConstantHeatTransfer)
    return ht_type.val
end

ext: jl


text:
include("network_model.jl")
using Statistics

# %%

n = Network(1000, 0.1, 800)

# %%

function post(neuron,connection)
    return findall(x->x==1,connection[neuron,:]) # change for new_net_col
end

const n1 = 2
const n2 = post(n1,n.connection)[1]
const interval = 20 #msec

n.s[n1,n2] = 0.0
n1f = [-100]
n2f = Int64[]
rew = Int64[]

function reward!(n1f::Array{Int64},n2f::Array{Int64},rew::Array{Int64},spiked::Array{Int64},time::Int64) # Module DA_STDP
    if n1 in spiked
        append!(n1f,time)
    end
    if n2 in spiked
        append!(n2f,time)
        if (time-last(n1f)<interval) && (last(n2f)>last(n1f))
            append!(rew,time+1000+rand(1:2000))
        end
    end
end

shist = zeros(10*3601,3)

idx = findall(x -> x != 0, n.connection)

# %%

@inbounds for sec in 0:3600
    @show sec
    @time @inbounds for msec in 1:1000
        time = 1000*sec+msec
        train = false
        if msec%10 == 0
            train = true
        end
        spiked = step!(n,train)
        reward!(n1f,n2f,rew,spiked,time)
        if time in rew
            n.da += 0.5
        end
        if msec%100==0
            shist[div(time,100),1:2] .= n.s[n1,n2],n.sd[n1,n2] # change for new_net_col
            shist[div(time,100),3] = mean(n.s[idx])
        end
    end
end

# %% Plot learning of the targeted synapse
using Plots
gr()
x1 = 0.1.*collect(1:length(shist[:,1]))
y1 = shist[:,1]
x2 = x1
y2 = shist[:,2]
y3 = shist[:,3]
fig = plot()
plot!(x1,y1,color="blue",label="synapse weight", legend = true)
plot!(x2,y2,color="green",label="eligibilty trace", legend = true)
plot!(x2, y3, color = "grey", label = "mean weight", legend = true)
xlabel!("Time (sec)")

max_stars_repo_name: t-cormier/da-stdp

max_issues_repo_name: t-cormier/da-stdp

max_forks_repo_name: t-cormier/da-stdp


text:
# A symbol is a name and the concolicly computed type
struct Sym
    name::Symbol
    _type::DataType
    Sym(base, _type::DataType) = new(Base.gensym(Symbol(base, '#', _type)), _type)
    Sym(_type::DataType) = new(Base.gensym(Symbol(_type)), _type)
end
Cassette.metadatatype(::Type{<:TraceCtx}, ::DataType) = Sym

record(args, ctx) = map(args) do arg
    if istagged(arg, ctx)
        metadata(arg, ctx)
    else
        arg
    end
end

# A callsite is a function with a concolic argument set
# children are all the callsites within the function
mutable struct Callsite
    f::Any
    args::Any
    retval::Union{Some{Any}, Nothing}
    children::Vector{Callsite}
    Callsite(f, args) = new(f, args, nothing, Callsite[])
end
Base.push!(trace::Callsite, call::Callsite) = push!(trace.children, call)

struct Metadata
    trace::Callsite
    # For current execution
    record::Vector{Sym}
    substitutes::Dict{DataType, Any}
end
Metadata() = Metadata(
    Callsite(:toplevel, ()),
    Vector{Sym}(),
    Dict{DataType, Vector}())

Metadata(m::Metadata, call::Callsite) = Metadata(call, m.record, m.substitutes)
Metadata(subs) = Metadata(Callsite(:toplevel, ()), Vector{Sym}(), subs)

function record!(ctx, type)
    m = ctx.metadata
    sym = Sym(:fval, type)
    push!(m.record, sym)

    if haskey(m.substitutes, type)
        subs = m.substitutes[type]
        if !isempty(subs)
            return pop!(subs), sym
        end
    end
    return nothing, sym
end
Base.push!(m::Metadata, call) = push!(m.trace, call)

function augment(record, subs)
    substitutions = Dict{DataType, Vector}()
    function insert!(type, val)
        if !haskey(substitutions, type)
            substitutions[type] = Vector{type}()
        end
        push!(substitutions[type], val)
    end

    for sym in record
        insert!(sym._type, subs[sym.name])
    end
    return substitutions
end

function Cassette.overdub(ctx::TraceCtx, ::typeof(rand), ::Type{T}) where T<:INTEGERS
    val, sym = record!(ctx, T)
    if val === nothing
        val = rand(T)
    end
    return tag(val::T, ctx, sym)
end

# oftype will drop tag on the x value since we are only looking at it'd type.
function Cassette.overdub(ctx::TraceCtx, ::typeof(oftype), x::Tagged, y)
    x = untag(x, ctx)
    return oftype(x, y)
end

##
# We need to manually override for IntrinsicFunctions (which are the leaf-nodes we are interested in)
# Since in tagged contexts there is an automatic fallback available.
##
function Cassette.overdub(ctx::TraceCtx{Metadata, <:Cassette.Tag}, f::Core.IntrinsicFunction, args...)
    call = Callsite(f, record(args, ctx))
    push!(ctx.metadata, call)
    retval = fallback(ctx, f, args...)
    if any(a -> istagged(a, ctx), args)
        retval = tag(retval, ctx, Sym(typeof(retval)))
    end
    call.retval = Some(istagged(retval, ctx) ? metadata(retval, ctx) : retval)
    return retval
end

##
# Recursivly trace though our program.
#
# Note: This won't trace functions that are defined primitive and there are several
#       fallbacks that Cassette provides for tagged contexts.
#
# Question: Shouldn't canoverdub for primitives be false?
##
function Cassette.overdub(ctx::TraceCtx, f, args...)
    # We need to push the callsite first into metadata
    # otherwise we run into issues with functions that throw errors
    call = Callsite(f, record(args, ctx))
    push!(ctx.metadata, call)

    # Recurse into the next function
    retval = if canrecurse(ctx, f, args...)
        newctx = similarcontext(ctx, metadata = Metadata(ctx.metadata, call))
        Cassette.recurse(newctx, f, args...)
    else
        retval = fallback(ctx, f, args...)
        # If any of our inputs was tagged we want to tag the return value so that
        # the symbolic nature of this operation is captured.
        if any(a -> istagged(a, ctx), args) && !istagged(retval, ctx)
            retval = tag(retval, ctx, Sym(typeof(retval)))
        end
        retval
    end
    call.retval = Some(istagged(retval, ctx) ? metadata(retval, ctx) : retval)
    return retval
end


text:
using UniqueVectors: UniqueVector
using LinearAlgebra: issymmetric, diag, I
using SparseArrays: spzeros

export Circuit
export partial_copy, matrices
export get_inv_inductance, get_inductance
export get_conductance, get_resistance
export get_elastance, get_capacitance
export set_inv_inductance!, set_inductance!
export set_conductance!, set_resistance!
export set_capacitance!, set_elastance!

export SpanningTree, coordinate_matrix, unite_vertices

SparseMat{T} = SparseMatrixCSC{T, Int}

"""
    struct Circuit{T, K<:Number, G<:Number, C<:Number}
Model of a circuit, containing matrix representations of the inverse inductance,
conductance, and capacitance between any pair of named vertices. All vertices are
named.

## Fields
- `k`: inverse inductance matrix
- `g`: conductance matrix
- `c`: capacitance matrix
- `vertices`: vector of uniquely named vertices

## Note
The matrices are not in the conventional form, i.e. `c` is not a Maxwell capacitance
matrix. Rather, `c[i,j]` is only the direct capacitance between vertices `i` and `j`.
`iszero(c[i,i])` is demanded for all `i`.
"""
struct Circuit{T, K<:Number, G<:Number, C<:Number}
    k::SparseMat{K} # inverse inductances
    g::SparseMat{G} # conductances
    c::SparseMat{C} # capacitances
    vertices::UniqueVector{T}
    function Circuit{T}(k::SparseMat, g::SparseMat, c::SparseMat,
            vertices::UniqueVector{T}) where {T}
        for m in (k, g, c)
            @assert isreal(m) && issymmetric(m)
            @assert size(m) == (length(vertices), length(vertices))
            @assert all(iszero.(diag(m)))
            @assert all(m .>= zero(eltype(m)))
        end
        kk, gg, cc = float(k), float(g), float(c)
        K, G, C = eltype(kk), eltype(gg), eltype(cc)
        return new{T,K,G,C}(kk, gg, cc, vertices)
    end
end

"""
    Circuit(k::AbstractMatrix, g::AbstractMatrix, c::AbstractMatrix,
        vertices::AbstractVector{T}) where {T}
    Circuit(vertices::AbstractVector)

Construct a circuit with named vertices, defaulting to zero inverse inductance,
conductance, and capacitance between the vertices if the matrices are not given.
"""
function Circuit(k::AbstractMatrix, g::AbstractMatrix, c::AbstractMatrix,
        vertices::AbstractVector{T}) where {T}
    return Circuit{T}(sparse(float(k)), sparse(float(g)), sparse(float(c)),
        UniqueVector(vertices))
end

function Circuit(vertices::AbstractVector)
    z() = spzeros(length(vertices), length(vertices))
    return Circuit(z(), z(), z(), vertices)
end

"""
    partial_copy(circ::Circuit;
       k::Union{AbstractMatrix, Nothing} = nothing,
       g::Union{AbstractMatrix, Nothing} = nothing,
       c::Union{AbstractMatrix, Nothing} = nothing,
       vertices::Union{AbstractVector, Nothing} = nothing)

Create a new `Circuit` with the same fields except those given as keyword arguments.
"""
function partial_copy(circ::Circuit;
       k::Union{AbstractMatrix, Nothing} = nothing,
       g::Union{AbstractMatrix, Nothing} = nothing,
       c::Union{AbstractMatrix, Nothing} = nothing,
       vertices::Union{AbstractVector, Nothing} = nothing)
   k = isnothing(k) ? circ.k : k
   g = isnothing(g) ? circ.g : g
   c = isnothing(c) ? circ.c : c
   vertices = isnothing(vertices) ? circ.vertices : vertices
   return Circuit(k, g, c, vertices)
end

"""
    ==(circ1::Circuit, circ2::Circuit)
Test whether two circuits have the same circuit components and named vertices
(in order).
"""
# don't remove space between Base. and ==
function Base. ==(circ1::Circuit, circ2::Circuit)
    return all(getfield(circ1, name) == getfield(circ2, name)
        for name in fieldnames(Circuit))
end

"""
    isapprox(circ1::Circuit, circ2::Circuit)
Tests for approximate equality of the circuits' `k`, `g`, `c` matrices.
"""
function Base.isapprox(circ1::Circuit, circ2::Circuit)
    return all(getfield(circ1, name) ≈ getfield(circ2, name)
        for name in fieldnames(Circuit)[1:end-1])
end

matrices(c::Circuit) = (c.k, c.g, c.c)

"""
    get_matrix_element(c::Circuit{T}, matrix_name::Symbol, v0::T, v1::T) where {T}
Returns a matrix element from one of the circuit matrices corresponding to the
edge between vertices with names `v0` and `v1`. `matrix_name` must be in
`(:k, :g, :c)`.
"""
function get_matrix_element(c::Circuit{T}, matrix_name::Symbol, v0::T, v1::T) where {T}
    @assert matrix_name in (:k, :g, :c)
    i0 = findfirst(isequal(v0), c.vertices)
    i1 = findfirst(isequal(v1), c.vertices)
    return getfield(c, matrix_name)[i0, i1]
end

"""
    get_inv_inductance(c::Circuit{T,K}, v0, v1) where {T,K}
Returns the inverse of the inductance on the edge between between vertices `v0` and `v1`.
"""
get_inv_inductance(c::Circuit{T, K}, v0, v1) where {T, K} =
    get_matrix_element(c, :k, v0, v1)::K

"""
    get_inductance(c::Circuit, v0, v1)
Returns the inductance on the edge between between vertices `v0` and `v1`.
"""
get_inductance(c::Circuit, v0, v1) = inv(get_inv_inductance(c, v0, v1))

"""
    get_conductance(c::Circuit{T,K,G}, v0, v1) where {T,K,G}
Returns the conductance on the edge between between vertices `v0` and `v1`.
"""
get_conductance(c::Circuit{T, K, G}, v0, v1) where {T, K, G} =
    get_matrix_element(c, :g, v0, v1)::G

"""
    get_resistance(c::Circuit, v0, v1)
Returns the resistance on the edge between between vertices `v0` and `v1`.
"""
get_resistance(c::Circuit, v0, v1) = inv(get_conductance(c, v0, v1))

"""
    get_capacitance(c::Circuit{T,K,G,C}, v0, v1) where {T,K,G,C}
Returns the capacitance on the edge between between vertices `v0` and `v1`.
"""
get_capacitance(c::Circuit{T,K,G,C}, v0, v1) where {T,K,G,C} =
    get_matrix_element(c, :c, v0, v1)::C

"""
    get_elastance(c::Circuit, v0, v1)
Returns the elastance on the edge between between vertices `v0` and `v1`.
"""
get_elastance(c::Circuit, v0, v1) = inv(get_capacitance(c, v0, v1))

"""
    set_matrix_element!(c::Circuit, matrix_name::Symbol, v0, v1, value)
Sets a matrix element in one of the circuit matrices corresponding to the
edge between vertices with names `v0` and `v1`. `matrix_name` must be in
`(:k, :g, :c)`.
"""
function set_matrix_element!(c::Circuit, matrix_name::Symbol, v0, v1, value)
    @assert matrix_name in (:k, :g, :c)
    @assert value >= zero(value)
    @assert !isinf(value)
    i0 = findfirst(isequal(v0), c.vertices)
    i1 = findfirst(isequal(v1), c.vertices)
    @assert i0 != i1
    matrix = getfield(c, matrix_name)
    matrix[i0, i1] = value
    matrix[i1, i0] = value
    return c
end

"""
    set_inv_inductance!(c::Circuit, v0, v1)
Sets the inverse of the inductance on the edge between between vertices `v0` and `v1`.
"""
set_inv_inductance!(c::Circuit, v0, v1, value) =
    set_matrix_element!(c, :k, v0, v1, value)

"""
    set_inductance!(c::Circuit, v0, v1)
Sets the inductance on the edge between between vertices `v0` and `v1`.
"""
set_inductance!(c::Circuit, v0, v1, value) =
    set_inv_inductance!(c, v0, v1, inv(value))

"""
    set_conductance!(c::Circuit, v0, v1)
Sets the conductance on the edge between between vertices `v0` and `v1`.
"""
set_conductance!(c::Circuit, v0, v1, value) =
    set_matrix_element!(c, :g, v0, v1, value)

"""
    set_resistance!(c::Circuit, v0, v1)
Sets the resistance on the edge between between vertices `v0` and `v1`.
"""
set_resistance!(c::Circuit, v0, v1, value) =
    set_conductance!(c, v0, v1, inv(value))

"""
    set_capacitance!(c::Circuit, v0, v1)
Sets the capacitance on the edge between between vertices `v0` and `v1`.
"""
set_capacitance!(c::Circuit, v0, v1, value) =
    set_matrix_element!(c, :c, v0, v1, value)

"""
    set_elastance!(c::Circuit, v0, v1)
Sets the elastance on the edge between between vertices `v0` and `v1`.
"""
set_elastance!(c::Circuit, v0, v1, value) =
    set_capacitance!(c, v0, v1, inv(value))

#######################################################
# spanning trees
#######################################################

"""
    struct SpanningTree{T}
A spanning tree of a graph with vertices of type `T`.

## Fields
- `root`: the root of the spanning tree
- `edges`: a vector of tuples `(vertex1, vertex2)` defining tree edges.
"""
struct SpanningTree{T}
    root::T
    edges::Vector{Tuple{T, T}}
end

"""
    SpanningTree(vertices::AbstractVector{T}) where T
Construct a depth-1 spanning tree from a vector of vertices of a complete graph,
with the first vertex taken as the root of the spanning tree.
"""
function SpanningTree(vertices::AbstractVector{T}) where T
    @assert !isempty(vertices)
    root = first(vertices)
    edges = [(v, root) for v in vertices[eachindex(vertices)[2:end]]]
    return SpanningTree{T}(root, edges)
end

"""
    coordinate_matrix(c::Circuit{T}, tree::SpanningTree{T}) where T
Returns the `T` matrix described in section IIC of arXiv:1810.11510 for a
circuit `c` with a given spanning `tree`. The direction of the edges passed in
is ignored and enforced as being outwards from the root.
"""
function coordinate_matrix(c::Circuit{T}, tree::SpanningTree{T}) where T
    V = length(c.vertices)
    @assert tree.root in c.vertices
    @assert length(tree.edges) == V-1 # the correct number of edges for a spanning tree
    mapping = Dict(tree.root => spzeros(V-1)) # the function l in equation 2
    edge_set = Set(enumerate(tree.edges))
    while !isempty(edge_set)
        for e in edge_set
            (index, (v0, v1)) = e
            @assert !(v0 in keys(mapping) && v1 in keys(mapping)) # edges contains a loop
            if v0 in keys(mapping) || v1 in keys(mapping)
                delete!(edge_set, e)
                if v0 in keys(mapping) # enforces all edges are oriented away from the root
                    mapping[v1] = 1.0 * mapping[v0]
                    mapping[v1][index] += 1
                else
                    mapping[v0] = 1.0 * mapping[v1]
                    mapping[v0][index] += 1
                end
            end
        end
    end
    @assert length(mapping) == V # all vertices have been assigned a coordinate
    return hcat([mapping[v] for v in c.vertices]...)
end

"""
    coordinate_matrix(num_vertices::Int)
Returns the `T` matrix described in section IIC of arXiv:1810.11510 for a
circuit with `num_vertices`, where the spanning tree is depth-1 and rooted at
the first vertex of the circuit.
"""
coordinate_matrix(num_vertices::Int) = [spzeros(num_vertices - 1) I]

"""
    coordinate_matrix(num_vertices::Int)
Returns the `T` matrix described in section IIC of arXiv:1810.11510 for a
circuit `c`, where the spanning tree is depth-1 and rooted at the first vertex
of the circuit.
"""
coordinate_matrix(c::Circuit) = coordinate_matrix(length(c.vertices))

#######################################################
# cascade and unite
#######################################################

"""
    cascade(circs::AbstractVector{<:Circuit})
    cascade(circs::Circuit...)
Cascade the circuits into one larger block diagonal circuit.
"""
function cascade(circs::AbstractVector{<:Circuit})
    circuit_matrices = [cat(m..., dims=(1,2))
        for m in zip([matrices(circ) for circ in circs]...)]
    vertices = vcat([circ.vertices for circ in circs]...)
    return Circuit(circuit_matrices..., vertices)
end
cascade(circs::Circuit...) = cascade(collect(circs))

"""
    unite_vertices(circ::Circuit{T}, vertices::AbstractVector{T}) where T
    unite_vertices(circ::Circuit{T}, vertices::T...) where T
Unite the given vertices of a circuit.
"""
function unite_vertices(circ::Circuit{T}, vertices::AbstractVector{T}) where T
    length(vertices) <= 1 && return circ
    n = length(circ.vertices)
    vertex_inds = [findfirst(isequal(v), circ.vertices) for v in vertices]
    keep_inds = filter(!in(vertex_inds[2:end]), 1:n)
    function unite_indices(mat::AbstractMatrix)
        m = copy(mat)
        ind1 = vertex_inds[1]
        for ind2 in vertex_inds[2:end]
            m[ind1, :] += mat[ind2, :]
            m[:, ind1] += mat[:, ind2]
        end
        m = m[keep_inds, keep_inds]
        for i in 1:size(m, 1)
            m[i, i] = 0
        end
        return m
    end
    circuit_matrices = [unite_indices(mat) for mat in matrices(circ)]
    return Circuit(circuit_matrices..., circ.vertices[keep_inds])
end
unite_vertices(circ::Circuit{T}, vertices::T...) where T =
    unite_vertices(circ, collect(vertices))

"""
    connect(models::AbstractVector{<:Circuit})
    connect(models::Circuit...)
Cascade all models and unite ports with the same name. This results in a
combined circuit with common ports connected.
"""
function connect(circs::AbstractVector{<:Circuit})
    @assert !isempty(circs)
    if length(circs) == 1
        return first(circs)
    end
    # number the vertices so that the names are all distinct and then cascade
    vertex_number = 1
    function rename(circ)
        vertices = [(vertex_number + i - 1, vertex)
            for (i, vertex) in enumerate(circ.vertices)]
        vertex_number += length(vertices)
        return Circuit(matrices(circ)..., vertices)
    end
    circ = cascade(map(rename, circs))
    # merge all vertices with the same name
    original_vertices = vcat([c.vertices for c in circs]...)
    for vertex in unique(original_vertices)
        inds = findall([v[2] == vertex for v in circ.vertices])
        circ = unite_vertices(circ, circ.vertices[inds])
    end
    # remove numbering
    return Circuit(matrices(circ)..., [v[2] for v in circ.vertices])
end
connect(circs::Circuit...) = connect(collect(circs))


text:
# model initializations (normal or xaiver)

# xavier initialization
function initx(d...; ftype=Float32)
    if gpu() >=0
        KnetArray{ftype}(xavier(d...))
    else
        Array{ftype}(xavier(d...))
    end
end

# random normal initialization
function initr(d...; ftype=Float32, GPUFEATS=false)
    if GPUFEATS && gpu() >=0
        KnetArray{ftype}(0.1*randn(d...))
    else
        Array{ftype}(0.1*randn(d...))
    end
end

# zero initialization
function initzeros(d...; ftype=Float32, GPUFEATS=false)
    if GPUFEATS && gpu() >=0
        KnetArray{ftype}(zeros(d...))
    else
        Array{ftype}(zeros(d...))
    end
end

# initialize buffer lstm
# Possible bufemebed: fvec o bvec o wvec o upos (o xpos o feats)
function initbuf(bufembed, bufhidden; buftype=:lstm, numLayers=1)
    r, wr = rnninit(bufembed, bufhidden, rnnType=buftype, numLayers=numLayers)
end

# commonly used decision module
function mlp(w,x; pdrop=(0,0))
    x = dropout(x,pdrop[1])
    for i=1:2:length(w)-2
        x = relu.(w[i]*x .+ w[i+1])
        x = dropout(x,pdrop[2])
    end
    return w[end-1]*x .+ w[end]
end


text:
using Test

@testset "./src/strategies.jl" begin

using BellScenario

@testset "Strategy()" begin
    @testset "with scenario" begin
    probs = [1 0 0;0 1 0;0 0 1]

        scenario = BlackBox(3,3)

        S = Strategy(probs, scenario)

        @test S == probs
        @test S isa Strategy
        @test S.scenario isa BlackBox
        @test S.scenario == scenario
    end

    @testset "without scenario, initializes as BlackBox" begin
        S = Strategy([0.5 0.5 0.5;0.5 0.5 0.5])

        @test S == 1/2*[1 1 1;1 1 1]

        @test S.scenario == BlackBox(2,3)
    end

    @testset "Strategy DomainError" begin
        scenario_error = BlackBox(3,4)

        @test_throws DomainError Strategy([1 0;0 1], scenario_error)
        @test_throws DomainError Strategy([0.5 0.5;1 0.5])
    end

    @testset "absolute tolerance" begin
        @test_throws DomainError Strategy([1/2 1/2; 1/2 1/3])
        @test Strategy([1/2 1/2; 1/2 1/3], atol=0.2) == [1/2 1/2;1/2 1/3]

        scenario = BlackBox(2,3)
        @test_throws DomainError Strategy([1 0 0;0 0 0], scenario)
        @test Strategy([1 0 0;0 0 0], scenario, atol=1.1) == [1 0 0;0 0 0]
    end
end

@testset "Strategy multiplication" begin
    @testset "Strategy * Strategy" begin
        S1 = Strategy([1 0;0 1;0 0])
        S2 = Strategy([1 0 1;0 1 0])

        S_prod1 = S1*S2
        S_prod2 = S2*S1

        @test S_prod1 == [1 0 1;0 1 0;0 0 0]
        @test S_prod1 isa Strategy
        @test S_prod1.scenario == BlackBox(3,3)

        @test S_prod2 == [1 0;0 1]
        @test S_prod2 isa Strategy
        @test S_prod2.scenario == BlackBox(2,2)

        scenario = BipartiteSignaling((3,1),(1,3))
        S_prod = *(S1, S2, scenario)

        @test S_prod == [1 0 1;0 1 0;0 0 0]
        @test S_prod isa Strategy
        @test S_prod.scenario == BipartiteSignaling((3,1),(1,3))

        @test_throws DomainError (*(S2, S1, scenario))
    end

    @testset "DeterministicStrategy * DeterministicStrategy" begin
        S1 = DeterministicStrategy([1 0;0 1;0 0])
        S2 = DeterministicStrategy([1 0 1;0 1 0])

        S_prod1 = S1*S2
        S_prod2 = S2*S1

        @test S_prod1 == [1 0 1;0 1 0;0 0 0]
        @test S_prod1 isa DeterministicStrategy
        @test S_prod1.scenario == BlackBox(3,3)

        @test S_prod2 == [1 0;0 1]
        @test S_prod2 isa DeterministicStrategy
        @test S_prod2.scenario == BlackBox(2,2)

        scenario = BipartiteSignaling((3,1),(1,3))
        S_prod = *(S1, S2, scenario)

        @test S_prod == [1 0 1;0 1 0;0 0 0]
        @test S_prod isa DeterministicStrategy
        @test S_prod.scenario == BipartiteSignaling((3,1),(1,3))

        @test_throws DomainError (*(S2, S1, scenario))
    end

    @testset "DeterministicStrategy * Strategy is stochastic" begin
        S1 = Strategy([1 0;0 1;0 0])
        S2 = DeterministicStrategy([1 0 1;0 1 0])

        S_prod1 = S1*S2
        S_prod2 = S2*S1

        @test S_prod1 == [1 0 1;0 1 0;0 0 0]
        @test S_prod1 isa Strategy
        @test S_prod1.scenario == BlackBox(3,3)

        @test S_prod2 == [1 0;0 1]
        @test S_prod2 isa Strategy
        @test S_prod2.scenario == BlackBox(2,2)

        scenario = BipartiteSignaling((3,1),(1,3))
        S_prod = *(S1, S2, scenario)

        @test S_prod == [1 0 1;0 1 0;0 0 0]
        @test S_prod isa Strategy
        @test S_prod.scenario == BipartiteSignaling((3,1),(1,3))

        @test_throws DomainError (*(S2, S1, scenario))
    end
end

@testset "random_strategy()" begin
    @testset "simple cases"  begin
        strat = random_strategy(3,3)

        @test strat isa Strategy
        @test size(strat) == (3,3)

        strat = random_strategy(9,7)

        @test strat isa Strategy
        @test size(strat) == (7,9)
    end
end

@testset "strategy_dims()" begin
    @testset "black box scenario" begin
        @test strategy_dims(BlackBox(5,2)) == (5,2)
    end

    @testset "prepare and measure" begin
        @test strategy_dims(LocalSignaling(3,5,2)) == (5,3)
    end

    @testset "bipartite non-signaling scenario" begin
        @test strategy_dims(BipartiteNonSignaling(2,2,2,2)) == (4,4)
        @test strategy_dims(BipartiteNonSignaling(2,3,4,5)) == (6,20)
    end

    @testset "bipartite scenario" begin
        @test strategy_dims(BipartiteSignaling((3,1),(1,4))) == (3,4)
        @test strategy_dims(BipartiteSignaling((2,2),(2,2))) == (4,4)
    end
end

end


text:
function module_dir()
    joinpath(@__DIR__, "..")
end

using Pkg
Pkg.activate(module_dir())
using RoboDojo


text:
@testset "Ensure transformations to the real line/model space are valid" begin
    for T in subtypes(Transform)
        global u = parameter(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.), T(), fixed=false)
        @test ( transform_to_real_line(u) |> x -> transform_to_model_space(u,x) ) == u.value

        if !isa(T,Type{ModelConstructors.Untransformed})
            # check transform_to_real_line and transform_to_model_space to different things if T is not ModelConstructors.Untransformed
            @test transform_to_real_line(u,u.value) != transform_to_model_space(u,u.value)
        end
    end
end

tomodel_answers = zeros(3)
toreal_answers  = zeros(3)
a = 1e-8; b = 5.; c = 1.
x = 2.5230
cx = 2 * (x - (a+b)/2)/(b-a)
toreal_answers[1] = 1. / (c * (x - a))
toreal_answers[2] = (1/c) * (1. / (1. - cx^2)^(-3/2)) * (2/(b-a))
toreal_answers[3] = 1.
x = transform_to_real_line(parameter(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.),
                                     ModelConstructors.Exponential(), fixed=false))
tomodel_answers[1] = c * exp(c * (x - b))
x = transform_to_real_line(parameter(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.),
                                     ModelConstructors.SquareRoot(), fixed=false))
tomodel_answers[2] = (b - a) / 2 * c / (1 + c^2 * x^2)^(3/2)
tomodel_answers[3] = 1.
@testset "Ensure derivatives of transformations to the real line/model space are valid" begin
    for (i,T) in enumerate(subtypes(Transform))
        global u = parameter(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.), T(), fixed=false)
        @test differentiate_transform_to_real_line(u, u.value) == toreal_answers[i]
        x = transform_to_real_line(u)
        @test differentiate_transform_to_model_space(u,x) == tomodel_answers[i]

        if !isa(T,Type{ModelConstructors.Untransformed})
            # check transform_to_real_line and transform_to_model_space to different things if T is not ModelConstructors.Untransformed
            @test differentiate_transform_to_real_line(u,u.value) != differentiate_transform_to_model_space(u,u.value)
        end
    end

end

# probability
N = 10^2
u = parameter(:bloop, 2.5230, (1e-8, 5.), (1e-8, 5.), ModelConstructors.SquareRoot(); fixed = true)
v = parameter(:cat, 2.5230, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), Gamma(2.00, 0.1))

pvec =  ParameterVector{Float64}(undef, N)
for i in 1:length(pvec)
	pvec[i] = (i%2 == 0) ? u : v
end
@testset "Check logpdf/pdf function approximations" begin
    @test logpdf(pvec) ≈ 50*logpdf(v)
    @test pdf(pvec) ≈ exp(50*logpdf(v))
end

updated = update(pvec, ones(length(pvec)))
ModelConstructors.update!(pvec, ones(length(pvec)))

@testset "Check if update! preserves dimensions and values" begin
    @test all(updated .== pvec)
    @test logpdf(pvec) == logpdf(updated)
end

# test we only update non-fixed parameters
@testset "Ensure only non-fixed parameters are updated" begin
    for p in pvec
        if p.fixed
            @test p.value == 2.5230
        elseif isa(p, Parameter)
            @test p.value == one(Float64)
        end
    end
end

# vector of new values must be the same length
@testset "Ensure update! enforces the same length of the parameter vector being updated" begin
    @test_throws AssertionError ModelConstructors.update!(pvec, ones(length(pvec)-1))
end

@testset "Ensure parameters being updated are of the same type." begin
    for w in [parameter(:moop, 3.0, fixed=false), parameter(:moop, 3.0; scaling = log, fixed=false)]
        # new values must be of the same type
        @test_throws ErrorException parameter(w, one(Int))

        # new value is out of bounds
        @test_throws ParamBoundsError parameter(w, -1.)
    end
end

@testset "Ensure parameter value types can be changed when forced by keyword." begin
    for w in [parameter(:moop, 1.0, fixed=false), parameter(:moop, 1.0; scaling = log, fixed=false)]
        # new values must be of the same type
        @test typeof(parameter(w, one(Int); change_value_type = true).value) == Int
    end
end

# subspecs
function sstest(m::AnSchorfheide)

    # Change all the fields of an unfixed parameter
    m <= parameter(:ι_w, 0.000, (0.0, .9999), (0.0,0.9999), ModelConstructors.Untransformed(), Normal(0.0,1.0), fixed=false,
                   description="ι_w: A new parameter.",
                   tex_label="\\iota_w")

    # Change an unfixed parameter to be fixed
    m <= parameter(:ι_p, 0.000, fixed=true,
                   description= "ι_p: Another new parameter",
                   tex_label="\\iota_p")

    # Change a fixed parameter
    m <= parameter(:δ, 0.02,  fixed=true,
                   description="δ: The capital depreciation rate.", tex_label="\\delta" )

    # Overwrite a fixed parameter with an unfixed parameter
    m <= parameter(:ϵ_p, 0.750, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.75, 0.4),        fixed=false,  scaling = x -> 1 + x/100,
                   description="ϵ_p: No description available.",
                   tex_label="\\varepsilon_{p}")

    steadystate!(m)
end

m = AnSchorfheide()
sstest(m)

@testset "Test steady-state parameters" begin
    @test m[:ι_w].value == 0.0
    @test m[:ι_w].valuebounds == (0.0, .9999)
    @test m[:ι_w].transform == ModelConstructors.Untransformed()
    @test m[:ι_w].transform_parameterization == (0.0,0.9999)
    @test isa(m[:ι_w].prior.value, Normal)

    @test m[:ι_p].value == 0.0
    @test m[:ι_p].valuebounds == (0.0, 0.0)
    @test isnull(m[:ι_p].prior)
    @test m[:ι_p].fixed == true

    @test m[:δ].value == 0.02
    @test m[:δ].valuebounds == (0.02, 0.02)
    @test isnull(m[:δ].prior)
    @test m[:δ].fixed == true

    @test m[:ϵ_p].value == 0.750
    @test m[:ϵ_p].transform == ModelConstructors.Exponential()
    @test isa(m[:ϵ_p].prior.value, Gamma)
    @test m[:ϵ_p].fixed==false
end

# Tests moved from DSGE.jl core.jl
# UnscaledParameter, fixed=false
α =  parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.05), fixed=false)
@testset "Test non-fixed UnscaledParameter" begin
    @test isa(α, UnscaledParameter)
    @test α.key == :α
    @test isa(α.prior.value, Normal)
    @test α.prior.value.μ == 0.3
    @test α.description == "No description available."
    @test α.tex_label == ""
    @test isa(α.transform, ModelConstructors.SquareRoot)
end

# UnscaledParameter, fixed = true
α_fixed =  parameter(:α_fixed, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.Untransformed(), Normal(0.30, 0.05), fixed=true)
@testset "Test fixed UnscaledParameter" begin
    @test α_fixed.transform_parameterization == (0.1596,0.1596)
    @test isa(α_fixed.transform, ModelConstructors.Untransformed)
end

# UnscaledParameter, fixed = true, transform should be overwritten given fixed
α_fixed =  parameter(:α_fixed, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.05), fixed=true)
@testset "Test fixed UnscaledParameter, ensuring transform is overwritten" begin
    @test isa(α_fixed.transform, ModelConstructors.Untransformed)
end

# Fixed UnscaledParameter, minimal constructor
δ = parameter(:δ, 0.025)
@testset "Test fixed UnscaledParameter minimal constructor" begin
    @test δ.fixed
    @test δ.transform_parameterization == (0.025, 0.025)
    @test δ.valuebounds == (0.025, 0.025)
end

# Scaled parameter
β = parameter(:β, 0.1402, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.25, 0.1), fixed=false,  scaling = x -> (1 + x/100)\1, description="β: Discount rate.", tex_label="\\beta ")
@testset "Test ScaledParameter constructor" begin
    @test isa(β, ScaledParameter)
    @test isa(β.prior.value, Gamma)
    @test isa(β.transform, ModelConstructors.Exponential)
end

# Invalid transform
@testset "Ensure error thrown on invalid transform" begin
    @test_throws UndefVarError α_bad = parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999),
                                                 InvalidTransform(), Normal(0.30, 0.05), fixed=false)
end

# Arithmetic with parameters
@testset "Check arithmetic with parameters" begin
    @test promote_type(AbstractParameter{Float64}, Float16) == Float64
    @test promote_type(AbstractParameter{Float64}, Int8) == Float64
    ## @test promote_rule(AbstractParameter{Float64}, Float16) == Float64
    ## @test promote_rule(AbstractParameter{Float64}, Int8) == Float64
    @test δ + δ == 0.05
    @test δ^2 == 0.025^2
    @test -δ == -0.025
    @test log(δ) == log(0.025)
end

# transform_to_real_line and transform_to_model_space
cx = 2 * (α - 1/2)
@testset "Check parameter transformations for optimization" begin
    @test abs(transform_to_real_line(α) - cx / sqrt(1 - cx^2)) < .001
    @test transform_to_real_line(δ) == 0.025
end

m = AnSchorfheide()
let lastparam = parameter(:p, 0.0)
    for θ in m.parameters
        isa(θ, Parameter) && (lastparam = θ)
    end
    @testset "Check AnSchorfheide last parameter" begin
        @test isa(lastparam, Parameter)
        @test lastparam.value == 0.20*2.237937
    end
end
# transform_to_real_line and transform_to_model_space, acting on the entire parameter vector. they should be inverses!
pvec = m.parameters
vals = transform_to_real_line(pvec)
transform_to_model_space!(m, vals)
@testset "Check parameter transformations for optimization part 2" begin
    @test pvec == m.parameters
end

# all fixed parameters should be unchanged by both transform_to_real_line and transform_to_model_space
@testset "Check fixed parameters are unchanged by optimization transformations" begin
    for θ in m.parameters
        if θ.fixed
            @test θ.value == transform_to_real_line(θ)
            @test θ.value == transform_to_model_space(θ, θ.value)
        end
    end
end

# Check parameters can be declared as vectors
@testset "Test vector-valued parameters" begin
    p = parameter(:test, ones(5), (0.0, 5.), (0., 5.0), Untransformed(), MvNormal(ones(5), ones(5)), fixed = false, scaling = x -> x/2)
    @test all(ones(5) ./ 2 .== parameter(p, ones(5)).scaledvalue)
end

# Check we can update only specific parameter values in a ParameterVector
@testset "Test updating specific values of a ParameterVector" begin
    # Check you can overwrite values which are unfixed
    pvec = ParameterVector{Float64}(undef, 3)
    pvec[1] = parameter(:a, 1., (0., 3.), (0., 3.), fixed = false)
    pvec[2] = parameter(:b, 1.)
    pvec[3] = parameter(:c, 1., (0., 3.), (0., 3.), fixed = false)
    vals = [2., 2.]
    update!(pvec, vals, BitArray([true, false, true]))
    @test map(x -> x.value, pvec) == [2., 1., 2.]

    # Check that overwriting fixed parameter values does not work, even if you pass
    # the BitArray telling update! to change the values.
    pvec[1] = parameter(:a, 1.)
    pvec[2] = parameter(:b, 1.)
    pvec[3] = parameter(:c, 1.)
    vals = [2., 2.]
    update!(pvec, vals, BitArray([true, false, true]))
    @test map(x -> x.value, pvec) == [1., 1., 1.]
end

nothing


text:
using InteractiveFixedEffectModels

tests = ["factormodel.jl",
		"interactivefixedeffectsmodel.jl"
		]

println("Running tests:")

for test in tests
	try
		include(test)
		println("\t\033[1m\033[32mPASSED\033[0m: $(test)")
	 catch e
	 	println("\t\033[1m\033[31mFAILED\033[0m: $(test)")
	 	showerror(stdout, e, backtrace())
	 	rethrow(e)
	 end
end


text:
function psiNext(p::Vector,psi::Vector)
    psi = cat(1,zero(Complex128),psi)
    N = length(psi)
    psik = zero(psi)
    if length(p)==1
        p = cat(1,zero(p[1]),p)
        c = sqrt(1-abs2(p[end]))
        #in this case we assume psi = u(t), rather than the previous psi
        for j=0:(N-1)
            psik[j+1] = c*sum([p[end]^(i)*psi[j-i+1] for i=0:(j-1)])
        end
    else
        c = sqrt((1-abs2(p[end]))/(1-abs2(p[end-1])))
        for j=0:(N-1)
            psik[j+1] = c*sum([p[end]^(i)*(psi[j-i]-conj(p[end-1])*psi[j-i+1]) for i=0:(j-1)])
        end
    end
    return psik[2:end]
end

function psiNextRecursive(p::Vector,psi::Vector)
    psi = cat(1,zero(Complex128),psi)
    N = length(psi)
    psik = zero(psi)
    if length(p)==1
        p = cat(1,zero(p[1]),p)
        c = sqrt(1-abs2(p[end]))
        #in this case we assume psi = u(t), rather than the previous psi
        for j=1:(N-1)
            psik[j+1] = p[end]*psik[j]+c*psi[j+1]
        end
    else
        c = sqrt((1-abs2(p[end]))/(1-abs2(p[end-1])))
        for j=1:(N-1)
            psik[j+1] = p[end]*psik[j] + c*(psi[j] - conj(p[end-1])*psi[j+1])
        end
    end
    return psik[2:end]
end

function genPsiBasis(p::Vector,u::Vector)
    N = length(p)
    B = zeros(Complex128,(length(u),N+1))
    B[:,1] = u
    for k=1:N
        B[:,k+1] = psiNext([p[1:k]],B[:,k])
    end
    return B[:,2:end]
end

function genPsiBasisRecursive(p::Vector,u::Vector)
    N = length(p)
    B = zeros(Complex128,(length(u),N+1))
    B[:,1] = u
    for k=1:N
        B[:,k+1] = psiNextRecursive([p[1:k]],B[:,k])
    end
    return B[:,2:end]
end

function pulseFun(k::Int)
    a = one(Complex128);
    return cat(1,a,zeros(Complex128,k-1))
end

function gaussianPulse(k::Int,fwhm::Float64)
    a = one(Complex128);
    x = [a*exp(-4*log(2)*(i^2/fwhm^2)) for i=1:k]
    return x/norm(x)
end

function unitStep(k::Int)
    return k>=0 ? 1 : 0
end

function dictElem(p::Complex128,K::Int)
    x = [(sqrt(1-abs2(p))*(p^(k-1))*unitStep(k-1))::Complex128 for k=0:(K)]
    return x[2:end]
end

function dictElemRecursive(p::Complex128,u::Vector)
    c = sqrt(1-abs2(p))
    u = cat(1,zero(Complex128),u)
    N = length(u)
    r = zero(u)
    for k=1:(N-1)
        r[k+1] = p*r[k]+c*u[k]
    end
    return r[2:end]
end

function filterR(p::Complex128,R::Vector)
    ai = 1/conj(p)
    R = cat(1,zero(Complex128),R)
    N = length(R)
    Rkp1 = zero(R)
    for k=0:(N-1)
        Rkp1[k+1]  = sum([ai^(i+1)*(p*R[k-i]-R[k-i+1]) for i=0:(k-1)])
    end
    return Rkp1[2:end]
end

function filterRrecursive(p::Complex128,R)
    ai = 1/conj(p)
    Rp = cat(1,zeros(Complex128,(1,size(R,2))),R)
    N = size(Rp,1)
    Rkp1 = zero(R)
    for j=1:size(R,2)
        Rkp1j = zero(Rp[:,j])
        for k=1:(N-1)
            Rkp1j[k+1]  = ai*Rkp1[k]-ai*(Rp[k+1]-p*Rp[k])
        end
        Rkp1[:,j] = Rkp1j[2:end]
    end
    return Rkp1
end

function dictObj(x,g,target)
    p = x[1]*exp(im*2*pi*x[2])
    trial = dictElem(p,size(target,1))
    r = sum(abs2(trial'*target),2)
    return r[1]
end

function dictObjNbest(x,g,target)
    p = [(x[i]*exp(im*2*pi*x[i+1]))::Complex128 for i=1:2:length(x)]
    R = copy(target)
    for k=1:length(p)
        trial = dictElem(p[k],size(target,1))
        #R = filterR(p[k],R - trial*((trial'*R)[1]))
        R = filterRrecursive(p[k],R - trial*((trial'*R)[1]))
    end
    r = norm(R)
    return r
end

function nthRemainder(p,target)
    R = copy(target)
    L = size(target,2);
    for k=1:length(p)
        trial = dictElem(p[k],size(target,1))
        for j=1:L
            R[:,j] = filterRrecursive(p[k],R[:,j] - dot(trial,R[:,j])*trial)
        end
    end
    return R
end

function unitDiskGrid(n::Int,delta::Float64)
    N = nextpow2(n)
    M = N>>1
    v = linspace(delta,one(Float64)-delta,M)
    Zjk = zeros(Complex128,N,N)
    for j=1:(M)
             #bottom, running counter-clockwise
            Zjk[(M+j),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(5*pi/4,7*pi/4,2*j))
            #top, running clockwise
            Zjk[(M-j+1),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(3*pi/4,1*pi/4,2*j))
    end
    for j=1:(M)
            #right, running clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M+j] = v[j]*exp(im*linspace(pi/4,-pi/4,2*j)) #right
            #left, running counter-clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M-j+1] = v[j]*exp(im*linspace(3*pi/4,5*pi/4,2*j)) #left
    end
    return reshape(Zjk,N^2)
end

function afdT(target::Array,Npoles::Int,Nstarts::Int)
    Z = unitDiskGrid(int(sqrt(Nstarts)),0.2)
    opt = Opt(:LN_SBPLX,2)
    ftol_rel!(opt,1e-6) #local stopping criteria
    lb = [0.01; 0]
    ub = [0.999; 1-eps(Float64)]
    lower_bounds!(opt,lb)
    upper_bounds!(opt,ub)
    avec = zeros(Complex128,(length(Z),2))
    ovec = zeros(Float64,length(Z))
    p = zeros(Complex128,Npoles)
    a = zeros(Complex128,(Npoles,size(target,2)))
    gk = copy(target)
    for j=1:Npoles
        max_objective!(opt,(x,g)->dictObj(x,g,gk)[1])
        for k=1:length(Z)
            beta0 = [abs(Z[k]); ((angle(Z[k])+pi)/(2*pi))]
            (ovec[k],avec[k,:],e) = optimize(opt,beta0)
        end
        aopt = avec[indmax(ovec),:]
        p[j] = aopt[1]*exp(im*2*pi*aopt[2])
        ep = dictElem(p[j],size(target,1))
        a[j,:] = ep'*gk
        #for i=1:size(target,2)
        #    gk[:,i] = filterRrecursive(p[j],(gk-a[j,i]*ep))
        #end
        gk = filterRrecursive(p[j],(gk-ep*a[j,:]))
    end
    return (p,a)
end

function afdTCyclic(target::Array,polesIn::Vector,Nstarts::Int,CycleTol::Float64,CycleMax::Int)
    Z = unitDiskGrid(int(sqrt(Nstarts)),0.01)
    Npoles = length(polesIn)
    opt = Opt(:LN_SBPLX,2)
    ftol_rel!(opt,1e-12) #local stopping criteria
    lb = [0.01; 0]
    ub = [0.9999999999; 1-eps(Float64)]
    lower_bounds!(opt,lb)
    upper_bounds!(opt,ub)
    avec = zeros(Complex128,(length(Z),2))
    ovec = zeros(Float64,length(Z))
    p = zero(polesIn)
    prel = CycleTol+1;
    cycle_count = 0;
    while ((prel>CycleTol) && (cycle_count<CycleMax))
        cycle_count += 1
        for j=1:Npoles
            pold = copy(p)
            pprime = pold[(1:Npoles).!=j]
            gk = nthRemainder(pprime,target)
            max_objective!(opt,(x,g)->dictObj(x,g,gk)[1])
            for k=1:length(Z)
                beta0 = [abs(Z[k]); ((angle(Z[k])+pi)/(2*pi))]
                (ovec[k],avec[k,:],errcode) = optimize(opt,beta0)
            end
            aopt = avec[indmax(ovec),:]
            p[j] = aopt[1]*exp(im*2*pi*aopt[2])
            prel = norm(p-pold)/norm(pold)
            #println(prel)
        end
    end
    B = genPsiBasisRecursive(p,pulseFun(size(target,1)))
    a = B\target
    #println(cycle_count)
    return (p,a)
end

function poles2vec(p::Vector{Complex128})
    r = zeros(Float64,length(p)*2)
    r[1:2:end] = abs(p)
    r[2:2:end] = (angle(p)+pi)/(2*pi)
    return r
end

function afdTNbest(target::Array,Npoles::Int,gWallClockTime::Number,beta0::Vector)
    opt = Opt(:GN_MLSL_LDS,Npoles*2)
    lopt = Opt(:LN_SBPLX,Npoles*2)
    ftol_rel!(lopt,1e-6) #local stopping criteria
    maxtime!(opt,gWallClockTime) #global stopping critera
    lb = [0.01; 0]
    ub = [0.99999; 1-eps(Float64)]
    lb = repmat(lb,Npoles)
    ub = repmat(ub,Npoles)
    lower_bounds!(opt,lb)
    upper_bounds!(opt,ub)
    local_optimizer!(opt,lopt)
    p = zeros(Complex128,Npoles)
    a = zeros(Complex128,(Npoles,size(target,2)))
    gk = copy(target)
    min_objective!(opt,(x,g)->dictObjNbest(x,g,gk)[1])
    (r,popt,e) = optimize(opt,beta0)
    p = [(popt[i]*exp(im*2*pi*popt[i+1]))::Complex128 for i=1:2:(Npoles*2)]
    B = genPsiBasis(p,pulseFun(size(target,1)))
    a = B\target
    return (p,a)
end

function akaikeInformationCriterion(RSS,numParams,dataLen)
    #This gives the sample size corrected Akaike Information Criteria
    #commonly known as AICc in the literature
    K = numParams+1
    return dataLen*log(RSS/dataLen)+2*K+2*K*((K+1)/(dataLen-K-1))
end

function processDataStack(Y::Array{Complex128,2},minNpoles::Int,maxNpoles::Int,Nstarts::Int,Textrap::Int)
    pmat = {}
    amat = {}
    Bmat = {}
    for j=1:size(Y,2)
        AICc = {}
        ptemp   = {}
        atemp   = {}
        for k=minNpoles:maxNpoles
            (p,a) = afdT(Y[:,j],k,Nstarts)
            push!(ptemp,p)
            push!(atemp,a)
            B = genPsiBasis(p,pulseFun(size(Y,1)))
            Yhat = squeeze(B*a,2)
            R = (Y[:,j]-Yhat)
            RSS = abs(dot(R,R))
            aicval = akaikeInformationCriterion(RSS,4*k,2*size(Y,1))
            push!(AICc,aicval)
        end
        optInd = indmin(AICc)
        push!(pmat,ptemp[optInd])
        push!(amat,atemp[optInd])
        B = genPsiBasis(pmat[j],pulseFun(Textrap))
        push!(Bmat,B)
        println(optInd+minNpoles-1)
        println(AICc[optInd])
    end
    That = zeros(Complex128,(Textrap,size(Y,2)))
    for k=1:size(Y,2)
        That[:,k] = Bmat[k]*amat[k]
    end
    return (That,pmat,amat)
end

function processDataStackCyclic(Y::Array{Complex128,2},minNpoles::Int,maxNpoles::Int,maxNcycles::Int,Nstarts::Int,CycleTol::Float64,Textrap::Int)
    pmat = {}
    amat = {}
    Bmat = {}
    for j=1:size(Y,2)
        AICc = {}
        ptemp   = {}
        atemp   = {}
        for k=minNpoles:maxNpoles
            (pinit,ainit) = afdT(Y[:,j],k,Nstarts)
            (p,a) = afdTCyclic(Y[:,j],pinit,Nstarts,CycleTol,maxNcycles)
            push!(ptemp,p)
            push!(atemp,a)
            B = genPsiBasis(p,pulseFun(size(Y,1)))
            Yhat = B*a
            R = (Y[:,j]-Yhat)
            RR = R'*R
            RSS = abs(RR)[1]
            aicval = akaikeInformationCriterion(RSS,2*length(p)+2*length(a),2*size(Y,1))
            push!(AICc,aicval)
        end
        optInd = indmin(AICc)
        push!(pmat,ptemp[optInd])
        push!(amat,atemp[optInd])
        B = genPsiBasis(pmat[j],pulseFun(Textrap))
        push!(Bmat,B)
        println(optInd+minNpoles-1)
        println(AICc)
    end
    That = zeros(Complex128,(Textrap,size(Y,2)))
    for k=1:size(Y,2)
        That[:,k] = Bmat[k]*amat[k]
    end
    return (That,pmat,amat)
end

function processDataStackCyclicGlobal(Y::Array{Complex128,2},minNpoles::Int,maxNpoles::Int,maxNcycles::Int,Nstarts::Int,CycleTol::Float64,Textrap::Int)
    Bmat = {}
    AICc = {}
    ptemp   = {}
    atemp   = {}
    for k=minNpoles:maxNpoles
        (pinit,ainit) = afdT(Y,k,Nstarts)
        (p,a) = afdTCyclic(Y,pinit,Nstarts,CycleTol,maxNcycles)
        push!(ptemp,p)
        push!(atemp,a)
        B = genPsiBasis(p,pulseFun(size(Y,1)))
        Yhat = B*a
        R = (Y-Yhat)
        RSS = sum((abs(R'*R)))
        aicval = akaikeInformationCriterion(RSS,2*length(p)+2*length(a),2*length(Y))
        println(aicval)
        push!(AICc,aicval)
    end
    optInd = indmin(AICc)
    popt = ptemp[optInd]
    aopt = atemp[optInd]
    B = genPsiBasis(popt,pulseFun(Textrap))
    #println(optInd+minNpoles-1)
    #println(AICc[optInd])
    That = B*aopt
    return (That,popt,aopt)
end

function processDataStackMLSL(Y::Array{Complex128,2},minNpoles::Int,maxNpoles::Int,Nstarts::Int,SearchTime::Number,Textrap::Int)
    pmat = {}
    amat = {}
    Bmat = {}
    for j=1:size(Y,2)
        AICc = {}
        ptemp   = {}
        atemp   = {}
        for k=minNpoles:maxNpoles
            (pinit,ainit) = afdT(Y[:,j],k,Nstarts)
            (p,a) = afdTNbest(Y[:,j],k,SearchTime,poles2vec(pinit))
            push!(ptemp,p)
            push!(atemp,a)
            B = genPsiBasis(p,pulseFun(size(Y,1)))
            Yhat = B*a
            R = (Y[:,j]-Yhat)
            RSS = abs(dot(R,R))
            aicval = akaikeInformationCriterion(RSS,4*k,2*size(Y,1))
            push!(AICc,aicval)
        end
        optInd = indmin(AICc)
        push!(pmat,ptemp[optInd])
        push!(amat,atemp[optInd])
        B = genPsiBasis(pmat[j],pulseFun(Textrap))
        push!(Bmat,B)
        println(optInd+minNpoles-1)
        println(AICc[optInd])
    end
    That = zeros(Complex128,(Textrap,size(Y,2)))
    for k=1:size(Y,2)
        That[:,k] = Bmat[k]*amat[k]
    end
    return (That,pmat,amat)
end

function processDataStackNpole(Y::Array{Complex128,2},Nsing::Int,Npoles::Int,Nstarts::Int,Textrap::Int)
    pmat = {}
    amat = {}
    Bmat = {}
    for j=1:Nsing
        (p,a) = afdT(Y[:,j],k,Nstarts)
        B = genPsiBasis(p,pulseFun(size(Y,1)))
        push!(pmat,p)
        push!(amat,a)
        B = genPsiBasis(pmat[j],pulseFun(Textrap))
        push!(Bmat,B)
    end
    That = zeros(Complex128,(Textrap,Nsing))
    for k=1:Nsing
        That[:,k] = Bmat[k]*amat[k]
    end
    return (That,pmat,amat)
end

function refineDataStack(Y::Array{Complex128,2},Nsing::Int,pmat,SearchTime::Number,Textrap::Int)
    Bmat = {}
    pref = {}
    aref = {}
    for j=1:Nsing
        pin = poles2vec(pmat[j])
        (p,a) = afdTNbest(Y[:,j],length(pmat[j]),SearchTime,pin)
        push!(pref,p)
        push!(aref,a)
        B = genPsiBasis(pref[j],pulseFun(Textrap))
        push!(Bmat,B)
    end
    That = zeros(Complex128,(Textrap,Nsing))
    for k=1:Nsing
        That[:,k] = Bmat[k]*aref[k]
    end
    return (That,pref,aref)
end

function refineDataStackCyclic(Y::Array{Complex128,2},Nsing::Int,pmat,Nstarts::Int,CycleTol::Float64,CycleMax::Int,Textrap::Int)
    Bmat = {}
    pref = {}
    aref = {}
    for j=1:Nsing
        (p,a) = afdTCyclic(Y[:,j],pmat[j],Nstarts,CycleTol,CycleMax)
        push!(pref,p)
        push!(aref,a)
        B = genPsiBasis(pref[j],pulseFun(Textrap))
        push!(Bmat,B)
    end
    That = zeros(Complex128,(Textrap,Nsing))
    for k=1:Nsing
        That[:,k] = Bmat[k]*aref[k]
    end
    return (That,pref,aref)
end

function unfoldMatrix(M::Array,N::Int)
    uM = M[:,1:N]
    for k=2:int(size(M,2)/N)
        uM = cat(1,uM,M[:,((k-1)*N+1):(k*N)]);
    end
    return uM
end

function cubeMatrix(M::Array,N::Int)
    cM = zeros(typeof(M[1]),(N,size(M,2),int(size(M,1)/N)))
    for k=1:int(size(M,1)/N)
        cM[:,:,k] = M[((k-1)*N+1):(k*N),:]
    end
    return cM
end

function fftaxis(N::Int,TSpac=1.0)
    if abs(mod(N,2))>0
        a = [0:1:int((N-1)/2)]
        b = [-(N-1)/2:1:-1]
        r = cat(1,a,b)/(N*TSpac)
    else
        a = [0:1:int(N/2-1)]
        b = [-int(N/2):1:-1]
        r = cat(1,a,b)/(N*TSpac)
    end
    return r
end

function unfoldCube(M::Array)
    R = zeros(typeof(M[1]),(size(M,1)*size(M,2),size(M,3)))
    N = size(M,1);
    for j=1:size(M,2)
        R[((j-1)*N+1):(j*N),:] = squeeze(M[:,j,:],2);
    end
    return R
end

function optimalSVDTrunc(u,s,v)
    N = size(v,1)
    M = size(v,1)
    if isreal(u[1])
        f = 1
    else
        f = 2;
    end
    AICc = {}
    Y = u*diagm(s)*v';
    for k=1:(N-1)
        Yhat = u[:,1:k]*diagm(s[1:k])*v[:,1:k]'
        R = Y-Yhat
        RSS = abs(R'*R)[1]
        push!(AICc,akaikeInformationCriterion(RSS,k,f*N*M))
    end
    optInd = indmin(AICc)
    return (AICc,optInd)
end

function w2wn(w::Vector)
    c = 299.792458; #nm/fs
    return 1E7*w/(2*pi*c)
end


text:
"""
    ArrayAllocators

Defines an array allocator interface and concrete array allocators using `malloc`, `calloc`, and memory alignment.

# Examples

```julia
using ArrayAllocators

Array{UInt8}(malloc, 100)
Array{UInt8}(calloc, 1024, 1024)
Array{UInt8}(MemAlign(2^16), (1024, 1024, 16))
```

See also `NumaAllocators`, `SafeByteCalculators`
"""
module ArrayAllocators

import Core: Array

export malloc, calloc

include("ByteCalculators.jl")

using .ByteCalculators

"""
    DefaultByteCalculator

Alias for [`ByteCalculators.CheckedMulByteCalculator`](@ref) representing the
byte calculator used with subtypes of `AbstractArrayAllocator` when one is not
specified.
"""
const DefaultByteCalculator = CheckedMulByteCalculator

export AbstractArrayAllocator, UndefAllocator, MallocAllocator, CallocAllocator
export MemAlign

"""
    lineage_finalizer(f, x)

Attach a finalizer function `f` to `x` or an ancestor of `x` using
`parent`.
"""
function lineage_finalizer(f, x)
    try
        finalizer(f, x)
    catch err
        if applicable(parent, x)
            lineage_finalizer(f, parent(x))
        else
            rethrow(err)
        end
    end
end

"""
    AbstractArrayAllocator{B}

Parent abstract type for array allocators. Parameter `B` is an AbstractByteCalculator
Defines `Array{T}(allocator, dims...) where T = Array{T}(allocator, dims)`
"""
abstract type AbstractArrayAllocator{B} end

# Allocate arrays when the dims are given as arguments rather than as a tuple
function (::Type{ArrayType})(alloc::A, dims...) where {T, ArrayType <: AbstractArray{T}, A <: AbstractArrayAllocator}
    return ArrayType(alloc, dims)
end

function (::Type{ArrayType})(alloc::A, dims) where {T, ArrayType <: AbstractArray{T}, B, A <: AbstractArrayAllocator{B}}
    num_bytes = nbytes(B{T}(dims))
    ptr = allocate(alloc, T, num_bytes)
    return unsafe_wrap(alloc, ArrayType, ptr, dims)
end

function allocate(alloc::A, ::Type{T}, num_bytes) where {A <: AbstractArrayAllocator, T}
    iszeroinit(A) || isbitstype(T) || throw(ArgumentError("$T is not a bitstype"))
    return Ptr{T}(allocate(alloc, num_bytes))
end
(::Type{A})(args...) where A <: AbstractArrayAllocator = A{DefaultByteCalculator}(args...)
iszeroinit(::Type{A}) where A <: AbstractArrayAllocator = false

"""
    UndefAllocator{B}

Allocate arrays using the builtin `undef` method. The `B` parameter is a `ByteCalculator`
"""
struct UndefAllocator{B} <: AbstractArrayAllocator{B}
end
allocate(::UndefAllocator, num_bytes) = C_NULL
Base.unsafe_wrap(::UndefAllocator, ::Type{ArrayType}, ::Ptr, dims::Dims) where {T, ArrayType <: AbstractArray{T}} = ArrayType(Core.undef, dims)

abstract type LibcArrayAllocator{B} <: AbstractArrayAllocator{B} end

"""
    wrap_libc_pointer(::Type{A}, ptr::Ptr{T}, dims) where {T, A <: AbstractArray{T}}
    wrap_libc_pointer(ptr::Ptr{T}, dims) where {T, A <: AbstractArray{T}}

Checks to see if `ptr` is C_NULL for an OutOfMemoryError.
Owns the array such that `Libc.free` is used.
"""
function wrap_libc_pointer(::Type{A}, ptr::Ptr{T}, dims) where {T, A <: AbstractArray{T}}
    if ptr == C_NULL
        throw(OutOfMemoryError())
    end
    # We use own = true, so we do not need Libc.free explicitly
    arr = unsafe_wrap(A, ptr, dims; own = true)
    return arr
end
wrap_libc_pointer(ptr::Ptr{T}, dims) where T = wrap_libc_pointer(Array{T}, ptr, dims)
Base.unsafe_wrap(alloc::A, args...) where A <: LibcArrayAllocator = wrap_libc_pointer(args...)

"""
    MallocAllocator()

Allocate array using [`Libc.malloc`](https://docs.julialang.org/en/v1/base/libc/#Base.Libc.malloc). This is not meant to be useful
but rather just to prototype the concept for a custom array allocator
concept. This should be similar to using `undef`.

See also https://en.cppreference.com/w/c/memory/malloc .
"""
struct MallocAllocator{B} <: LibcArrayAllocator{B}
end
allocate(::MallocAllocator, num_bytes) = Libc.malloc(num_bytes)

"""
    malloc

[`MallocAllocator`](@ref) singleton instance. `malloc` will only allocate memory. It does not initialize memory is is similar
in use as `undef`. See the type and the [C standard library function](https://en.cppreference.com/w/c/memory/malloc) for details.

# Example

```jldoctest
julia> Array{UInt8}(malloc, 16, 16);

```

"""
const malloc = MallocAllocator()

"""
    CallocAllocator()

Use [`Libc.calloc`](https://docs.julialang.org/en/v1/base/libc/#Base.Libc.calloc) to allocate an array. This is similar to `zeros`, except
that the Libc implementation or the operating system may allocate and
zero the memory in a lazy fashion.

See also https://en.cppreference.com/w/c/memory/calloc .
"""
struct CallocAllocator{B} <: LibcArrayAllocator{B}
end
allocate(::CallocAllocator, num_bytes) = Libc.calloc(num_bytes, 1)
iszeroinit(::Type{A}) where A <: CallocAllocator = true

"""
    calloc

[`CallocAllocator`](@ref) singleton instance. `calloc` will allocate memory and guarantee initialization to `0`.
See the type for details and the [C standard library function](https://en.cppreference.com/w/c/memory/calloc) for further details.

# Example

```jldoctest
julia> A = Array{UInt8}(calloc, 16, 16);

julia> sum(A)
0x0000000000000000
```
"""
const calloc = CallocAllocator()

"""
    MemAlign([alignment::Integer])

Allocate aligned memory. Alias for platform specific implementations.

`alignment` must be a power of 2.

On POSIX systems, `alignment` must be a multiple of `sizeof(Ptr)`.
On Windows, `alignment` must be a multiple of 2^16.

If `alignment` is not specified, it will be set to `min_alignment(MemAlign)`.

`MemAlign` is a constant alias for one the following platform specific implementations.
* POSIX (Linux and macOS): [`POSIX.PosixMemAlign`](@ref)
* Windows: [`Windows.WinMemAlign`](@ref).

"""
MemAlign

"""
    AbstractMemAlign{B} <: AbstractArrayAllocator{B}

Abstract supertype for aligned memory allocators.
"""
abstract type AbstractMemAlign{B} <: AbstractArrayAllocator{B} end

"""
    alignment(alloc::AbstractMemAlign)

Get byte alignment of the AbstractMemAlign array allocator.
"""
alignment(alloc::AbstractMemAlign) = alloc.alignment

"""
    min_alignment(::AbstractMemAlign)

Get the minimum byte alignment of the AbstractMemAlign array allocator.
"""
function min_alignment() end

include("Windows.jl")
include("POSIX.jl")

import .Windows: WinMemAlign
import .POSIX: PosixMemAlign

@static if Sys.iswindows()
    const MemAlign = WinMemAlign
elseif Sys.isunix()
    const MemAlign = PosixMemAlign
end

include("zeros.jl")

end # module ArrayAllocators


text:
using Test
using OrdinaryDiffEq, Calculus, ForwardDiff

function f(du,u,p,t)
  du[1] = -p[1]
  du[2] = p[2]
end

for x in 0:0.001:5
  called = false
  if x in [1.0, 2.0, 3.0, 4.0, 5.0]
    print("AD Ping $x")
  end
  function test_f(p)
    cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(called=true;integrator.p[2]=zero(integrator.p[2])))
    prob = ODEProblem(f,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-14,reltol=1e-14,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end
  p = [2.0, x]
  called = false
  findiff = Calculus.finite_difference_jacobian(test_f,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_f,p)
  @test called
  @test findiff ≈ fordiff
end

function f2(du,u,p,t)
  du[1] = -u[2]
  du[2] = p[2]
end

for x in 2.1:0.001:5
  called = false
  if x in [3.0, 4.0, 5.0]
    print("AD Ping $x")
  end
  function test_f2(p)
    cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(called=true;integrator.p[2]=zero(integrator.p[2])))
    prob = ODEProblem(f2,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end
  p = [2.0, x]
  findiff = Calculus.finite_difference_jacobian(test_f2,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_f2,p)
  @test called
  @test findiff ≈ fordiff
end

#=
#x = 2.0 is an interesting case

x = 2.0

function test_f2(p)
  cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(@show(x,integrator.t);called=true;integrator.p[2]=zero(integrator.p[2])))
  prob = ODEProblem(f2,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
  integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
  step!(integrator)
  solve!(integrator).u[end]
end

p = [2.0, x]
findiff = Calculus.finite_difference_jacobian(test_f2,p)
@test called
called = false
fordiff = ForwardDiff.jacobian(test_f2,p)
@test called

# At that value, it shouldn't be called, but a small perturbation will make it called, so finite difference is wrong!
=#

for x in 1.0:0.001:2.5
  if x in [1.5, 2.0, 2.5]
    print("AD Ping $x")
  end
  function lotka_volterra(du,u,p,t)
    x, y = u
    α, β, δ, γ = p
    du[1] = dx = α*x - β*x*y
    du[2] = dy = -δ*y + γ*x*y
  end
  u0 = [1.0,1.0]
  tspan = (0.0,10.0)
  p = [x,1.0,3.0,1.0]
  prob = ODEProblem(lotka_volterra,u0,tspan,p)
  sol = solve(prob,Tsit5())

  called=false
  function test_lotka(p)
    cb = ContinuousCallback((u,t,i) -> u[1]-2.5, (integrator)->(called=true;integrator.p[4]=1.5))
    prob = ODEProblem(lotka_volterra,eltype(p).([1.0,1.0]),eltype(p).((0.0,10.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end

  findiff = Calculus.finite_difference_jacobian(test_lotka,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_lotka,p)
  @test called
  @test findiff ≈ fordiff
end

# Gradients and Hessians

function myobj(θ)
  f(u,p,t) = -θ[1]*u
  u0, _ = promote(10.0, θ[1])
  prob = ODEProblem(f, u0, (0.0, 1.0))
  sol = solve(prob, Tsit5())
  diff = sol.u - 10*exp.(-sol.t)
  return diff'diff
end

ForwardDiff.gradient(myobj, [1.0])
ForwardDiff.hessian(myobj, [1.0])

function myobj2(θ)
  f(du,u,p,t) = (du[1]=-θ[1]*u[1])
  u0, _ = promote(10.0, θ[1])
  prob = ODEProblem(f, [u0], (0.0, 1.0))
  sol = solve(prob, Tsit5())
  diff = sol[:,1] .- 10 .*exp.(-sol.t)
  return diff'diff
end

ForwardDiff.gradient(myobj2, [1.0])
ForwardDiff.hessian(myobj2, [1.0])

function myobj3(θ)
  f(u,p,t) = -θ[1]*u
  u0, _ = promote(10.0, θ[1])
  tspan_end, _ = promote(1.0, θ[1])
  prob = ODEProblem(f, u0, (0.0, tspan_end))
  sol = solve(prob, Tsit5())
  diff = sol.u - 10*exp.(-sol.t)
  return diff'diff
end

ForwardDiff.gradient(myobj3, [1.0])
ForwardDiff.hessian(myobj3, [1.0])

function myobj4(θ)
  f(du,u,p,t) = (du[1] = -θ[1]*u[1])
  u0, _ = promote(10.0, θ[1])
  tspan_end, _ = promote(1.0, θ[1])
  prob = ODEProblem(f, [u0], (0.0, tspan_end))
  sol = solve(prob, Tsit5())
  diff = sol[:,1] .- 10 .* exp.(-sol.t)
  return diff'diff
end

ForwardDiff.gradient(myobj4, [1.0])
ForwardDiff.hessian(myobj4, [1.0])

f1s = function (du,u,p,t)
    du[1] = p[2] * u[1]
    du[2] = p[3] * u[2]
    du[3] = p[4] * u[3]
    du[4] = p[5] * u[4]
    nothing
end
f2s = function (du,u,p,t)
    du[1] = p[1]*u[2]
    du[2] = p[1]*u[3]
    du[3] = p[1]*u[4]
    du[4] = p[1]*u[1]
    nothing
end
u0 = [3.4, 3.3, 3.2, 3.1]
params = [0.002, -0.005, -0.004, -0.003, -0.002]
tspan = (7.0, 84.0)
times = collect(minimum(tspan):0.5:maximum(tspan))
prob = SplitODEProblem(f1s,f2s, u0, tspan, params)
sol2 = solve(prob, KenCarp4(); dt=0.5, saveat=times)

function difffunc(p)
    tmp_prob = remake(prob,p=p)
    vec(solve(tmp_prob,KenCarp4(),saveat=times))
end
ForwardDiff.jacobian(difffunc,ones(5))


text:
using Test, FinancialMonteCarlo, DifferentialEquations;
@show "Differential Equation Junctor"
Nsim = 10000;
Nstep = 30;
toll = 0.8;
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
S0 = 100.0
K = 100.0;
D = 90.0;
r = 0.02
sigma = 0.2
T = 1.0;
d = 0.01;
u0 = S0;
#Drift
dr_(u, p, t) = (r - d) * u
#Diffusion
g_1(u, p, t) = sigma * u
#Time Window
tspan = (0.0, T)
#Definition of the SDE
prob = SDEProblem(dr_, g_1, u0, tspan, rng = mc.rng)
monte_prob = EnsembleProblem(prob)
rfCurve = ZeroRate(r);
model = FinancialMonteCarlo.MonteCarloDiffEqModel(monte_prob, Underlying(S0, d))

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)

@show FwdPrice = pricer(model, rfCurve, mc, FwdData);
@show EuPrice = pricer(model, rfCurve, mc, EUData);
@show AmPrice = pricer(model, rfCurve, mc, AMData);
@show BarrierPrice = pricer(model, rfCurve, mc, BarrierData);
@show AsianPrice1 = pricer(model, rfCurve, mc, AsianFloatingStrikeData);
@show AsianPrice2 = pricer(model, rfCurve, mc, AsianFixedStrikeData);

@test abs(FwdPrice - 98.8436678850961) < toll
@test abs(EuPrice - 8.17917706833563) < toll
@test abs(BarrierPrice - 7.08419994877601) < toll
@test abs(AsianPrice1 - 4.70419181812687) < toll

Tneg = -T;
tspanNeg = (0.0, Tneg)
probNeg = SDEProblem(dr_, g_1, u0, tspanNeg, rng = mc.rng)
monte_probNeg = EnsembleProblem(probNeg)


text:
# Environments of the CMPS
function leftenv(Ψ::FiniteCMPS{<:AbstractPiecewise}; Kmax = 50, tol = defaulttol(Ψ))
    Q, Rs, vL, vR = Ψ
    a, b = domain(Ψ)
    ρL, infoL = lefttransfer(vL*vL', nothing, Ψ; Kmax = Kmax, tol = tol)
    Z = vR'*ρL(b)*vR
    imag(Z) < defaulttol(Z) && real(Z) > 0 ||
        @warn "Non-positive left environment: Z = vR'*ρL(b)*vR = $Z"
    λ = log(abs(Z))/(b-a)/2
    return ρL, λ, infoL
end

function leftenv!(Ψ; Kmax = 50, tol = defaulttol(Ψ))
    ρL, λ, infoL = leftenv(Ψ)
    (a,b) = domain(Ψ)
    Ψ.Q -= λ * one(Ψ.Q)
    for i = 1:length(ρL)
        xᵢ = ρL.nodes[i]
        xⱼ = ρL.nodes[i+1]
        xc = offset(ρL[i])
        Δx = max(abs(xⱼ - xc), abs(xc - xᵢ))
        exptaylor = [exp(-2*λ*(xc-a))]
        for k = 1:Kmax
            push!(exptaylor, exptaylor[k]*(-2*λ)/k)
            abs(exptaylor[k+1]*Δx^k) < tol && break
        end
        ρL[i] = truncmul(TaylorSeries(exptaylor, xc), ρL[i];
                            Kmax = Kmax, tol = tol, dx = Δx)
    end
    return ρL, infoL
end

function rightenv(Ψ::FiniteCMPS{<:AbstractPiecewise}; Kmax = 50, tol = defaulttol(Ψ))
    Q, Rs, vL, vR = Ψ
    a, b = domain(Q)
    ρR, infoR = righttransfer(vR*vR', nothing, Ψ; Kmax = Kmax, tol = tol)
    Z = vL'*ρR(a)*vL
    imag(Z) < defaulttol(Z) && real(Z) > 0 ||
        @warn "Non-positive right environment: Z = vL'*ρR(a)*vL = $Z"
    λ = log(abs(Z))/(b-a)/2
    return ρR, λ, infoR
end

function rightenv!(Ψ; Kmax = 50, tol = defaulttol(Ψ))
    ρR, λ, infoR = rightenv(Ψ)
    Ψ.Q -= λ * one(Ψ.Q)
    for i = 1:length(ρR)
        xᵢ = ρR.nodes[i]
        xⱼ = ρR.nodes[i+1]
        xc = offset(ρR[i])
        Δx = max(abs(xⱼ - xc), abs(xc - xᵢ))
        exptaylor = [exp(2*λ*(xc-b))]
        for k = 1:Kmax
            push!(exptaylor, exptaylor[k]*(2*λ)/k)
            abs(exptaylor[k+1]*Δx^k) < tol && break
        end
        ρR[i] = truncmul(TaylorSeries(exptaylor, xc), ρR[i];
                            Kmax = Kmax, tol = tol, dx = Δx)
    end
    return ρR, infoR
end

function environments!(Ψ::FiniteCMPS; Kmax = 50, tol = defaulttol(Ψ))
    ρL, infoL = leftenv!(Ψ; Kmax = Kmax, tol = tol)
    ρR, λR, infoR = rightenv(Ψ; Kmax = Kmax, tol = tol)
    λR < 1000*tol ||
        @warn "Incompatible normalization between left and right environments: λR = $λR"
    localZ = localdot(ρL, ρR)
    (a, b) = domain(Ψ)
    Za = localZ(a)
    Zb = localZ(b)
    abs(Za-1) < sqrt(defaulttol(Za)) && abs(Zb-1) < sqrt(defaulttol(Zb)) ||
        @warn "Incompatible environment normalizations: Z = ⟨ρL|ρR⟩ = $Za = $Zb"
    return ρL, ρR, infoL, infoR
end

# Environments of the CMPS with the Hamiltonian
const FiniteCMPSData = Tuple{FiniteCMPS{<:PiecewiseLinear}, Piecewise, Piecewise}

function leftenv!(H::LocalHamiltonian, Ψ::FiniteCMPS; kwargs...)
    domain(H) == domain(Ψ) || throw(DomainMismatch())
    ρL, ρR, infoL, infoR = environments!(Ψ; kwargs...)
    return leftenv(H, (Ψ,ρL,ρR); kwargs...)
end

function rightenv!(H::LocalHamiltonian, Ψ::FiniteCMPS; kwargs...)
    domain(H) == domain(Ψ) || throw(DomainMismatch())
    ρL, ρR, infoL, infoR = environments!(Ψ; kwargs...)
    return rightenv(H, (Ψ,ρL,ρR); kwargs...)
end

# assumes Ψ is normalized and ⟨ρL|ρR⟩ = 1
function leftenv(H::LocalHamiltonian, Ψρs::FiniteCMPSData;
                    Kmax = 50, tol = defaulttol(Ψρs[1]))
    Ψ, ρL, ρR = Ψρs
    (a,b) = domain(Ψ)
    domain(H) == (a,b) || throw(DomainMismatch())
    hL = leftreducedoperator(H.h, Ψ, ρL)
    eL = localdot(hL, ρR)
    EL = integrate(eL, (a,b))
    HL, infoL = lefttransfer(zero(hL(a)), hL, Ψ; Kmax = Kmax, tol = tol)
    return HL, EL, eL, hL, infoL
end

function rightenv(H::LocalHamiltonian, Ψρs::FiniteCMPSData;
                    Kmax = 50, tol = defaulttol(Ψρs[1]))
    Ψ, ρL, ρR = Ψρs
    (a,b) = domain(Ψ)
    domain(H) == (a,b) || throw(DomainMismatch())
    hR = rightreducedoperator(H.h, Ψ, ρR)
    eR = localdot(ρL, hR)
    ER = integrate(eR, (a,b))
    HR, infoR = righttransfer(zero(hR(b)), hR, Ψ; Kmax = Kmax, tol = tol)
    return HR, ER, eR, hR, infoR
end

function environments!(H::LocalHamiltonian, Ψ::FiniteCMPS; kwargs...)
    ρL, ρR, infoρL, infoρR = environments!(Ψ; kwargs...)
    HL, HR, hL, hR, E, e, infoHL, infoHR = environments(H,(Ψ,ρL,ρR); kwargs...)
    return HL, HR, ρL, ρR, hL, hR, E, e, infoHL, infoHR, infoρL, infoρR
end

function environments(H::LocalHamiltonian, Ψρs::FiniteCMPSData; kwargs...)
    HL, EL, eL, hL, infoHL = leftenv(H, Ψρs; kwargs...)
    HR, ER, eR, hR, infoHR = rightenv(H, Ψρs; kwargs...)
    EL ≈ ER ||
        @warn "non-matching energy from left and right environments"
    E = (EL+ER)/2
    e = rmul!(eL + eR, 1//2)
    return HL, HR, hL, hR, E, e, infoHL, infoHR
end


text:
mutable struct ListContainer{T}
    store::Gtk.GtkListStore
    selection::Gtk.GtkTreeSelection
    tree::Gtk.GtkTreeView
    cols::Array{Gtk.GtkTreeViewColumn, 1}
    data::Array{T, 1}
    dataType::DataType
    visible::Array{Bool, 1}
end

const columnTypes = Dict{Type, Gtk.GObject}(
    Array{String, 1} => GtkCellRendererCombo(),
    Bool => GtkCellRendererToggle(),
    Gtk.GdkPixbufLeaf => GtkCellRendererPixbuf()
)

const columnNames = Dict{Type, String}(
    Array{String, 1} => "combo",
    Bool => "active",
    Gtk.GdkPixbufLeaf => "pixbuf"
)

const listViewCallbackUnion = Union{Bool, Function}
const listViewSelectUnion = Union{Integer, Function, Nothing}

# Broken in Gtk.jl, redefined here
index_from_iter(store::GtkListStore, iter::Gtk.TRI) = parse(Int, Gtk.get_string_from_iter(GtkTreeModel(store), iter)) + 1

sanitizeListData(t::Tuple) = t
sanitizeListData(a::Array{<: Tuple, 1}) = sanitizeListData.(a)
sanitizeListData(a::Array{T, 1}) where T = Tuple{T}[sanitizeListData(v) for v in a]
sanitizeListData(v) = (v,)

function getSelected(list::ListContainer, default::Union{Integer, Nothing}=nothing)
    if hasselection(list.selection)
        return index_from_iter(list.store, selected(list.selection))

    else
        return default
    end
end

function unselectRow!(list::ListContainer)
    if hasselection(list.selection)
        row = selected(list.selection)

        Gtk.unselect!(list.selection, row)
    end
end

function currentRow(list::ListContainer)
    if hasselection(list.selection)
        return index_from_iter(list.store, selected(list.selection))

    else
        return 0
    end
end

function selectRow!(list::ListContainer, i::Integer=1; force::Bool=false)
    if 1 <= i <= length(list.store) && (i != getSelected(list, -1) || force)
        Gtk.GLib.@sigatom Gtk.select!(list.selection, Gtk.iter_from_index(list.store, i))

        return true
    end

    return false
end

function selectRow!(list::ListContainer, f::Function, default::Number=1; force::Bool=false)
    for i in 1:length(list.store)
        row = list.store[i]

        if f(row)
            selectRow!(list, i, force=force)

            return true
        end
    end

    return selectRow!(list, default, force=true)
end

function getListData(container::ListContainer)
    return container.dataType[container.store[i] for i in 1:length(container.store)]
end

# Default to currently selected
function unsafeUpdateTreeView!(container::ListContainer{T}, data::Array{T, 1}, select::listViewSelectUnion=1; setData::Bool=true, updateByReplacement::Bool=false) where T <: Tuple
    rows, cols = size(container.store)

    if setData
        container.data = deepcopy(data)
        container.visible = fill(true, length(data))
    end

    # This makes sure nothing new is automaticly selected while emptying the store
    container.selection = GAccessor.mode(container.selection, Gtk.GConstants.GtkSelectionMode.NONE)

    if updateByReplacement
        if rows > length(data)
            for i in 1:rows - length(data)
                Gtk.GLib.@sigatom deleteat!(container.store, rows - i + 1)
            end
        end

        for (i, row) in enumerate(data)
            if i <= rows
                for j in 1:length(row)
                    container.store[i, j] = row[j]
                end

            else
                push!(container.store, row)
            end
        end

    else
        empty!(container.store)

        for row in data
            push!(container.store, row)
        end
    end

    container.selection = GAccessor.mode(container.selection, Gtk.GConstants.GtkSelectionMode.SINGLE)

    if !isempty(data) && select !== nothing
        Gtk.GLib.@sigatom selectRow!(container, select)
    end
end

function updateTreeView!(container::ListContainer, data::Array{T, 1}, select::listViewSelectUnion=1; setData::Bool=true, updateByReplacement::Bool=false) where T
    data = isempty(data) ? container.dataType[] : sanitizeListData.(data)
    Gtk.GLib.@sigatom unsafeUpdateTreeView!(container, data, select; setData=setData, updateByReplacement=updateByReplacement)
end

function filterContainer!(container::ListContainer, select::listViewSelectUnion=1)
    updateTreeView!(container, container.data[container.visible], select, setData=false)
end

function placeholderCallback(store, col, row, text)

end

function textRenderer(editable::Bool=false, callback::Tuple{GtkListStore, Integer, Function}=(GtkListStore(String), 1, placeholderCallback))
    renderer = GtkCellRendererText(editable=editable)
    @guarded signal_connect(renderer, :edited) do widget, row, text
        store, col, func = callback

        func(store, col, parse(Int, row) + 1, text)
    end

    return renderer
end

correctTupleType(::Type{Bool}) = Bool
correctTupleType(::Type{<: Integer}) = Int
correctTupleType(::Type{<: Real}) = Float64
correctTupleType(type::Type{T}) where T = type

function generateTreeView(header::Union{NTuple{N, String}, String}, data::Array{T, 1}; resize::Bool=true, sortable::Bool=true, editable::Array{Bool, 1}=fill(false, length(header)), callbacks::Array{listViewCallbackUnion, 1}=Array{listViewCallbackUnion, 1}(fill(false, length(header))), visible::Array{Bool, 1}=fill(true, length(header))) where {N, T <: Tuple}
    header = sanitizeListData(header)
    data = sanitizeListData(data)

    tupleTypesRaw = eltype(data).parameters
    tupleTypes = correctTupleType.(tupleTypesRaw)
    store = GtkListStore(tupleTypes...)

    for d in data
        push!(store, d)
    end

    treeView = GtkTreeView(GtkTreeModel(store))

    cols = GtkTreeViewColumn[]

    for c in 1:length(tupleTypes)
        col = GtkTreeViewColumn(
            header[c],
            get(columnTypes, tupleTypes[c], textRenderer(editable[c], (store, c, isa(callbacks[c], Function) ? callbacks[c] : placeholderCallback))),
            Dict(get(columnNames, tupleTypes[c], "text") => c - 1)
        )

        push!(cols, col)
    end

    if resize
        for col in cols
            GAccessor.resizable(col, true)
        end
    end

    if sortable
        for (i, col) in enumerate(cols)
            GAccessor.sort_column_id(col, i - 1)
            GAccessor.reorderable(col, i - 1)
        end
    end

    append!(treeView, cols)

    selection = GAccessor.selection(treeView)
    set_gtk_property!.(cols, :visible, visible)

    return ListContainer{T}(store, selection, treeView, cols, deepcopy(data), Tuple{tupleTypes...}, fill(true, length(data)))
end

function connectChanged(f::Function, container::ListContainer)
    @guarded signal_connect(container.selection, "changed") do widget
        if hasselection(container.selection)
            row = selected(container.selection)

            if applicable(f, container, container.store[row]...)
                f(container, container.store[row]...)

            else
                f(container, container.store[row])
            end
        end
    end
end

function connectDoubleClick(f::Function, container::ListContainer)
    @guarded signal_connect(container.tree, "row-activated") do widget, box, col
        if hasselection(container.selection)
            row = selected(container.selection)

            if applicable(f, container, container.store[row]...)
                f(container, container.store[row]...)

            else
                f(container, container.store[row])
            end
        end
    end
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn