text:
using JSom
using Base.Test

srand(1)

@testset "Grid SOM" begin
    gsom = GridSOM(5, 5, 1)
    gsom.weights = zeros(5, 5, 1)
    set_unit_weight(gsom, 3, 4, 5.0)
    set_unit_weight(gsom, 2, 2, 2.0)

    @testset "neighboring units" begin
        n1 = sort(JSom.neighbor_units(gsom, (2,2)))
        n2 = sort(JSom.neighbor_units(gsom, (1,5)))
        x, y = (2,2)
        @test n1 == sort([
            (x-1, y-1), (x-1, y), (x-1, y+1),
            (x, y-1),             (x, y+1),
            (x+1, y-1), (x+1, y), (x+1, y+1)
        ])
        x, y = (1,5)
        @test n2 == sort([
            (x, y-1),
            (x+1, y-1), (x+1, y)
        ])
    end

    @testset "U-matrix" begin
        u = umatrix(gsom)
        @test u[3,4] == get_unit_weight(gsom, (3,4))[1]
        @test u[2,2] == get_unit_weight(gsom, (2,2))[1]
        @test u[2,3] == mean([2, 5, 0, 0, 0, 0, 0, 0])
    end

    @testset "neighborhood functions" begin
        L = 5
        indices = [(x,y) for x in 1:L, y in 1:L]
        @testset "gaussian" begin
            gsom.ħ = Gaussian_Neighborhood()
            bell = reshape([JSom.ħ(gsom, u, (2,2), 1.0) for u in indices], (L,L))
            @test maximum(bell) == 1.0
            @test indmax(bell) == 7
            @test bell[2,1] == bell[2,3] == bell[1,2] == bell[3,2]
            @test bell[1,1] == bell[1,3] == bell[3,1] == bell[3,3]
            @test bell[2,2] ≥ bell[2,3] ≥ bell[2,4] ≥ bell[2,5]
        end

        @testset "mexican hat" begin
            gsom.ħ = Ricker_Neighborhood()
            bell = reshape([JSom.ħ(gsom, u, (2,2), 1.0) for u in indices], (L,L))
            @test maximum(bell) == 1.0
            @test indmax(bell) == 7
            @test bell[2,1] == bell[2,3] == bell[1,2] == bell[3,2] == 0
            @test bell[1,1] == bell[1,3] == bell[3,1] == bell[3,3] < 0
            @test bell[3,3] ≤ bell[4,4] ≤ bell[5,5]
        end

        @testset "triangular" begin
            gsom.ħ = Triangular_Neighborhood()
            bell = reshape([JSom.ħ(gsom, u, (2,2), 2.0) for u in indices], (L,L))
            @test maximum(bell) == 1.0
            @test indmax(bell) == 7
            @test bell[2,1] == bell[2,3] == bell[1,2] == bell[3,2] == 0.5
            @test bell[3,4] == bell[3,5] == 0.0
        end
    end
end

@testset "Hex SOM" begin
    hsom = HexSOM(5, 5, 1)
    hsom.weights = zeros(5, 5, 1)
    set_unit_weight(hsom, 3, 4, 5.0)
    set_unit_weight(hsom, 2, 2, 2.0)

    @testset "U-matrix" begin
        u = umatrix(hsom)
        @test u[3,4] == get_unit_weight(hsom, (3,4))[1]
        @test u[2,2] == get_unit_weight(hsom, (2,2))[1]
        @test u[2,3] == mean([2, 0, 0, 0, 0, 0])
    end

    @testset "neighboring units" begin
        n1 = sort(JSom.neighbor_units(hsom, (2,2)))
        n2 = sort(JSom.neighbor_units(hsom, (1,5)))
        x, y = (2,2)
        @test n1 == sort([
                        (x-1, y), (x-1, y+1),
            (x, y-1),             (x, y+1),
            (x+1, y-1), (x+1, y)
        ])
        x, y = (1,5)
        @test n2 == sort([
            (x, y-1),
            (x+1, y-1), (x+1, y)])
    end

    @testset "neighborhood functions" begin
        L = 5
        indices = [(x,y) for x in 1:L, y in 1:L]
        @testset "gaussian" begin
            hsom.ħ = Gaussian_Neighborhood()
            bell = reshape([JSom.ħ(hsom, u, (2,2), 1.0) for u in indices], (L,L))
            @test maximum(bell) == 1.0
            @test indmax(bell) == 7
            @test bell[2,1] == bell[2,3] == bell[1,2] == bell[3,2] == bell[1,3] == bell[3,1]
            @test bell[2,2] ≥ bell[2,3] ≥ bell[2,4] ≥ bell[2,5]
        end

        @testset "mexican hat" begin
            hsom.ħ = Ricker_Neighborhood()
            bell = reshape([JSom.ħ(hsom, u, (2,2), 1.0) for u in indices], (L,L))
            @test maximum(bell) == 1.0
            @test indmax(bell) == 7
            @test bell[2,1] == bell[2,3] == bell[1,2] == bell[3,2] == bell[1,3] == bell[3,1] == 0
            @test bell[1,1] == bell[3,3] < 0
            @test bell[3,3] ≤ bell[4,4] ≤ bell[5,5]
        end

        @testset "triangular" begin
            hsom.ħ = Triangular_Neighborhood()
            bell = reshape([JSom.ħ(hsom, u, (2,2), 2.0) for u in indices], (L,L))
            @test maximum(bell) == 1.0
            @test indmax(bell) == 7
            @test bell[2,1] == bell[2,3] == bell[1,2] == bell[3,2] == bell[1,3] == bell[3,1] == 0.5
            @test bell[3,4] == bell[3,5] == 0.0
        end
    end
end

@testset "Generic SOM" begin
    som = GridSOM(5, 5, 1)
    som.weights = zeros(5, 5, 1)
    set_unit_weight(som, 3, 4, 5.0)
    set_unit_weight(som, 2, 2, 2.0)
    data = rand(4, 2)

    @testset "bmu_map" begin
        bmus = bmu_map(som, [5,2])
        @test bmus[(3,4)] == Vector[[5]]
        @test bmus[(2,2)] == Vector[[2]]
        @test bmus[(1,1)] == []
        @test bmus[(2,4)] == []
    end

    @testset "activation_response" begin
        response = activation_response(som, [5,2])
        @test response[3,4] == 1
        @test response[2,2] == 1
        @test response[1,1] == 0
        @test response[2,4] == 0
    end

    @testset "activate" begin
        activate(som, [5.0,])
        @test indmin(som.activation_map) == 18
    end

    @testset "quantize" begin
        q = quantize(som, [4,2])
        @test q[1] == 5.0
        @test q[2] == 2.0
    end

    @testset "quantization_error" begin
        @test quantization_error(som, [5,2]) == 0.0
        @test quantization_error(som, [4,2]) == 0.5
    end

    @testset "random seed" begin
        som1 = GridSOM(5, 5, 2, seed=1)
        som2 = GridSOM(5, 5, 2, seed=1)
        @test som1.weights == som2.weights
        sequential_random_epoch(som1, data, 5)
        sequential_random_epoch(som2, data, 5)
        @test som1.weights == som2.weights
    end

    @testset "reset" begin
        som = GridSOM(5, 5, 2, seed=1)
        sequential_random_epoch(som, data, 5)
        q1 = quantization_error(som, data)
        reset_state(som)
        sequential_random_epoch(som, data, 5)
        q2 = quantization_error(som, data)
        @test q1 == q2
    end

    @testset "training" begin
        @testset "sequential_epoch" begin
            som = GridSOM(5, 5, 2)
            data = rand(4, 2)
            q1 = quantization_error(som, data)
            sequential_epoch(som, data)
            q2 = quantization_error(som, data)
            sequential_epoch(som, data)
            q3 = quantization_error(som, data)
            @test q1 ≥ q2
            @test q2 ≥ q3
        end

        @testset "sequential_random_epoch" begin
            som = GridSOM(5, 5, 2)
            data = rand(4, 2)
            q1 = quantization_error(som, data)
            sequential_random_epoch(som, data, 5)
            q2 = quantization_error(som, data)
            sequential_random_epoch(som, data, 5)
            q3 = quantization_error(som, data)
            @test q1 ≥ q2
            @test q2 ≥ q3
        end
    end
end

@testset "Decay functions" begin
    @testset "inverse" begin
        @test JSom.τ_inverse(1., 2, 3.) == 1 / (1 + 2/3)
    end

    @testset "exponential" begin
        @test JSom.τ_exponential(1., 2, 3.) == 1 * exp(-2/3)
    end
end

max_stars_repo_name: agrimaldi/JSom

max_issues_repo_name: grmwld/JSom

max_forks_repo_name: agrimaldi/JSom


text:
using Logging
using Random
using LinearAlgebra
using StatsBase
using Distributions
using CubicSplines
using LightGraphs
using Plots
using DelimitedFiles
using GraphPlot
using Compose
using Cairo
using HypothesisTests
using Clustering

# plot results of `D` against `C` for the 1D ring
d_data = readdlm("D-results-matrix.csv", ',')
d_data[1,1] = 4
d_data = Array{Integer}(d_data)
trendline = mean(d_data, dims=2)
c_range = 0.0:0.1:0.7
scatter(c_range, d_data, legend=false, c=:grey, xlabel="C", ylabel="D", markerstrokecolor = :grey,
	xtickfontsize=15, ytickfontsize=15, linewidth=2, xguidefontsize=18, yguidefontsize=18, fmt=:pdf)
plot!(c_range, trendline, c=:black, linewidth=2)
savefig(string("D-C-curve.pdf"))

# plot `E` FI curve
phi_x_train = [-0.015,      0, 0.025, 0.05, 0.075,  0.1, 0.15]
phi_y_train = [     0,  0.005, 0.033, 0.05, 0.06, 0.068, 0.08]
phi_spline = CubicSpline(phi_x_train, phi_y_train)
phi_curve = phi_spline[range(phi_x_train[1], stop=phi_x_train[end], length=1_000)]
phi_curve[phi_curve.<0.0] .= 0.0
xs = range(phi_x_train[1], stop=phi_x_train[end], length=1_000)
plot(xs, phi_curve, legend=false, c=:black, xlabel="I", ylabel="V",
	xtickfontsize=15, ytickfontsize=15, linewidth=2,
	xguidefontsize=18, yguidefontsize=18, fmt=:pdf)
savefig(string("E_F-I-curve.pdf"))

# plot `I` FI curve
function psi(x)
	return max(0, (x-0.05)*0.1)
end
xs_psi = range(0, stop=0.6, length=1_000)
psi_curve = psi.(xs_psi)

plot(xs_psi, psi_curve, legend=false, c=:black, xlabel="I", ylabel="V",
	xtickfontsize=15, ytickfontsize=15, linewidth=2,
	xguidefontsize=18, yguidefontsize=18, fmt=:pdf)
savefig(string("I_F-I-curve.pdf"))

# distance plots in arbitrary graphs

function cd_scatter(seed, g)
	correlations = readdlm(string("seed",seed,"correlations_array_selective.csv"), ',', Float64)
	diam_g = diameter(g)
	geodesics_array = zeros(Int, nv(g), nv(g))
	distances = zeros(nv(g)^2 - nv(g), 2)

	for v in vertices(g)
		geodesics_array[v,:] = gdistances(g,v)
	end

	k = 1
	for i in vertices(g)
		for j in vertices(g)
			if i!=j
				distance = geodesics_array[i,j]
				correlation = correlations[i,j]
				distances[k,1] = distance
				distances[k,2] = correlation
				k+=1
			end
		end
	end

	scatter(distances[:,1], distances[:,2], legend=false, c=:black, xlabel="Distance", ylabel="Correlation")
	savefig(string("seed",seed,"_distances-correlations-scatter.pdf"))
end

# geometric indicies calculations

function geometric_indicies(seed, g, selective)
	if selective
		correlations_array = readdlm(string("seed",seed,"correlations_array_selective.csv"), ',', Float64)
	else
		correlations_array = readdlm(string("seed",seed,"correlations_array.csv"), ',', Float64)
	end
	diam_g = diameter(g)
	graph_size = nv(g)
	geodesics_array = zeros(Int, graph_size, graph_size)
	distances = zeros(nv(g)^2 - nv(g), 2)

	for v in vertices(g)
		geodesics_array[v,:] = gdistances(g,v)
	end

	Rs = zeros(diam_g)

	for d=1:diam_g
		R = 0
		for mu in 1:graph_size
			for v in 1:graph_size
				if (mu != v)
					la = geodesics_array[mu,v]<=d ? 1 : -1
					R += la*correlations_array[mu, v]
				end
			end
		end
		R /= ((graph_size*graph_size)-graph_size)
		Rs[d] = R
	end

	if selective
		open(string("seed",seed,"_Rs_selective.txt"),"a") do io
		   println(io,"Rs ",Rs)
		end
	else
		open(string("seed",seed,"_Rs.txt"),"a") do io
		   println(io,"Rs ",Rs)
		end
	end
end

# examples
seed = 1636988619
g = generate_triK5_graph()
cd_scatter(seed, g)

seed = 1636950525
g = generate_triK5_graph()
geometric_indicies(seed, g, true) # selective neurons
geometric_indicies(seed, g, false) # all neurons

ext: jl

max_stars_repo_name: tfburns/BurnsHagaFukai

max_issues_repo_name: tfburns/BurnsHagaFukai

max_forks_repo_name: tfburns/BurnsHagaFukai


text:
using OptimTestProblems
using Base.Test

@test true


text:
# utilities not used in the final version, just for exploration

"""
Read the given columns from `io`, and return them in a vector of
`AutoEnum` structures.
"""
function read_categories(io::IO, columns::Vector{Int},
                         autoenums=[AutoEnum{AbstractString}() for c in columns])
    while !eof(io)
        fields = split(chomp(readline(io)), ";")
        for (column, autoenum) in zip(columns, autoenums)
            tryparse(autoenum, fields[column])
        end
    end
    autoenums
end


text:
module WidgetsBase

using Observables

to_node(x) = Observable(x)
to_node(x::Observable) = x

abstract type AbstractWidget{T} <: Observables.AbstractObservable{T} end

Observables.observe(x::AbstractWidget) = x.value

include("ui_elements.jl")

export Button, Slider, RangeSlider, TextField, NumberInput, Checkbox

end # module


text:
immutable UTF16String <: String
    data::Array{Uint16,1}
end

utf16_is_lead(c::Uint16) = (c & 0xfc00) == 0xd800
utf16_is_trail(c::Uint16) = (c & 0xfc00) == 0xdc00
utf16_is_surrogate(c::Uint16) = (c & 0xf800) == 0xd800
utf16_get_supplementary(lead::Uint16, trail::Uint16) = char((lead-0xd7f7)<<10 + trail)

function endof(s::UTF16String)
    d = s.data
    i = length(d)
    i == 0 && return i
    utf16_is_surrogate(d[i]) ? i-1 : i
end

function next(s::UTF16String, i::Int)
    if !utf16_is_surrogate(s.data[i])
        return char(s.data[i]), i+1
    elseif length(s.data) > i && utf16_is_lead(s.data[i]) && utf16_is_trail(s.data[i+1])
        return utf16_get_supplementary(s.data[i], s.data[i+1]), i+2
    end
    error("invalid UTF-16 character index")
end

# TODO: optmize this
function encode16(s::String)
    buf = Uint16[]
    for c in s
        if c < 0x10000
            push!(buf, uint16(c))
        else
            push!(buf, uint16(0xd7c0 + (c>>10) & 0x3ff))
            push!(buf, uint16(0xdc00 + c & 0x3ff))
        end
    end
    UTF16String(buf)
end

utf16(x) = convert(UTF16String, x)
convert(::Type{UTF16String}, s::UTF16String) = s
convert(::Type{UTF16String}, s::String) = encode16(s)
convert(::Type{UTF8String}, s::UTF16String) =
    sprint(length(s.data), io->for c in s; write(io,c::Char); end)
convert(::Type{Array{Uint16,1}}, s::UTF16String) = s.data
convert(::Type{Array{Uint16}}, s::UTF16String) = s.data

sizeof(s::UTF16String) = sizeof(s.data)
convert{T<:Union(Int16,Uint16)}(::Type{Ptr{T}}, s::UTF16String) =
    convert(Ptr{T}, s.data)

function is_valid_utf16(data::Array{Uint16})
    i = 1
    n = length(data)
    while i < n # check for unpaired surrogates
        if utf16_is_lead(data[i]) && utf16_is_trail(data[i+1])
            i += 2
        elseif utf16_is_surrogate(data[i])
            return false
        else
            i += 1
        end
    end
    return i > n || !utf16_is_surrogate(data[i])
end

is_valid_utf16(s::UTF16String) = is_valid_utf16(s.data)

function convert(::Type{UTF16String}, data::Array{Uint16})
    !is_valid_utf16(data) && throw(ArgumentError("invalid UTF16 data"))
    UTF16String(data)
end

function convert(T::Type{UTF16String}, bytes::Array{Uint8})
    isempty(bytes) && return UTF16String(Uint16[])
    isodd(length(bytes)) && throw(ArgumentError("odd number of bytes"))
    data = reinterpret(Uint16, bytes)
    # check for byte-order mark (BOM):
    if data[1] == 0xfeff        # native byte order
        convert(T, data[2:end])
    elseif data[1] == 0xfffe    # byte-swapped
        convert(T, Uint16[bswap(data[i]) for i=2:length(data)])
    else
        convert(T, copy(data)) # assume native byte order
    end
end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
using Mimi
using Random

export getpage

include("../../src/utils/load_parameters.jl")
include("../../src/utils/mctools.jl")

include("mcs_growth.jl")
include("compute_scc_growth.jl")

include("../../src/components/RCPSSPScenario.jl")
include("../../src/components/extensions/CO2emissions_growth.jl")
include("../../src/components/CO2cycle.jl")
include("../../src/components/extensions/CO2forcing_growth.jl")
include("../../src/components/extensions/CH4emissions_growth.jl")
include("../../src/components/CH4cycle.jl")
include("../../src/components/CH4forcing.jl")
include("../../src/components/extensions/N2Oemissions_growth.jl")
include("../../src/components/N2Ocycle.jl")
include("../../src/components/N2Oforcing.jl")
include("../../src/components/extensions/LGemissions_growth.jl")
include("../../src/components/LGcycle.jl")
include("../../src/components/LGforcing.jl")
include("../../src/components/SulphateForcing.jl")
include("../../src/components/TotalForcing.jl")
include("../../src/components/ClimateTemperature.jl")
include("../../src/components/SeaLevelRise.jl")
include("../../src/components/extensions/GDP_growth_lowpass.jl")
include("../../src/components/MarketDamages.jl")
include("../../src/components/MarketDamagesBurke.jl")
include("../../src/components/NonMarketDamages.jl")
include("../../src/components/Discontinuity.jl")
include("../../src/components/AdaptationCosts.jl")
include("../../src/components/extensions/SLRDamages_growth.jl")
include("../../src/components/AbatementCostParameters.jl")
include("../../src/components/AbatementCosts.jl")
include("../../src/components/TotalAbatementCosts.jl")
include("../../src/components/TotalAdaptationCosts.jl")
include("../../src/components/Population.jl")
include("../../src/components/extensions/EquityWeighting_growth.jl")
include("../../src/components/extensions/PermafrostSiBCASA_growth.jl")
include("../../src/components/extensions/PermafrostJULES_growth.jl")
include("../../src/components/PermafrostTotal.jl")

function buildpage(m::Model, scenario::String, use_permafrost::Bool=true, use_seaice::Bool=true, use_page09damages::Bool=false)

    # add all the components
    scenario = addrcpsspscenario(m, scenario)
    climtemp = addclimatetemperature(m, use_seaice)
    if use_permafrost
        permafrost_sibcasa = add_comp!(m, PermafrostSiBCASA)
        permafrost_jules = add_comp!(m, PermafrostJULES)
        permafrost = add_comp!(m, PermafrostTotal)
    end

    # Socio-Economics
    population = addpopulation(m)
    gdp = add_comp!(m, GDP) # one can change the names per @defcomp to normal names to make the changes throughout this file minimal from the original.

    co2emit = add_comp!(m, co2emissions)
    co2cycle = addco2cycle(m, use_permafrost)
    add_comp!(m, co2forcing)
    ch4emit = add_comp!(m, ch4emissions)
    ch4cycle = addch4cycle(m, use_permafrost)
    add_comp!(m, ch4forcing)
    n2oemit = add_comp!(m, n2oemissions)
    add_comp!(m, n2ocycle)
    add_comp!(m, n2oforcing)
    lgemit = add_comp!(m, LGemissions)
    add_comp!(m, LGcycle)
    add_comp!(m, LGforcing)
    sulfemit = add_comp!(m, SulphateForcing)
    totalforcing = add_comp!(m, TotalForcing)
    add_comp!(m, SeaLevelRise)

    # Abatement Costs
    abatementcostparameters_CO2 = addabatementcostparameters(m, :CO2)
    abatementcostparameters_CH4 = addabatementcostparameters(m, :CH4)
    abatementcostparameters_N2O = addabatementcostparameters(m, :N2O)
    abatementcostparameters_Lin = addabatementcostparameters(m, :Lin)

    set_param!(m, :automult_autonomoustechchange, .65)

    abatementcosts_CO2 = addabatementcosts(m, :CO2)
    abatementcosts_CH4 = addabatementcosts(m, :CH4)
    abatementcosts_N2O = addabatementcosts(m, :N2O)
    abatementcosts_Lin = addabatementcosts(m, :Lin)
    add_comp!(m, TotalAbatementCosts)

    # Adaptation Costs
    adaptationcosts_sealevel = addadaptationcosts_sealevel(m)
    adaptationcosts_economic = addadaptationcosts_economic(m)
    adaptationcosts_noneconomic = addadaptationcosts_noneconomic(m)
    add_comp!(m, TotalAdaptationCosts)

    # Impacts
    slrdamages = addslrdamages(m)
    marketdamages = addmarketdamages(m)
    marketdamagesburke = addmarketdamagesburke(m)
    nonmarketdamages = addnonmarketdamages(m)
    add_comp!(m, Discontinuity)

    # Equity weighting and Total Costs
    equityweighting = add_comp!(m, EquityWeighting)

    # connect parameters together
    connect_param!(m, :ClimateTemperature => :fant_anthroforcing, :TotalForcing => :fant_anthroforcing)

    population[:popgrw_populationgrowth] = scenario[:popgrw_populationgrowth]

    connect_param!(m, :GDP => :pop_population, :Population => :pop_population)
    gdp[:grw_gdpgrowthrate] = scenario[:grw_gdpgrowthrate]

    if use_permafrost
        permafrost_sibcasa[:rt_g] = climtemp[:rt_g_globaltemperature]
        permafrost_jules[:rt_g] = climtemp[:rt_g_globaltemperature]
        permafrost[:perm_sib_ce_co2] = permafrost_sibcasa[:perm_sib_ce_co2]
        permafrost[:perm_sib_e_co2] = permafrost_sibcasa[:perm_sib_e_co2]
        permafrost[:perm_sib_ce_ch4] = permafrost_sibcasa[:perm_sib_ce_ch4]
        permafrost[:perm_jul_ce_co2] = permafrost_jules[:perm_jul_ce_co2]
        permafrost[:perm_jul_e_co2] = permafrost_jules[:perm_jul_e_co2]
        permafrost[:perm_jul_ce_ch4] = permafrost_jules[:perm_jul_ce_ch4]
    end

    co2emit[:er_CO2emissionsgrowth] = scenario[:er_CO2emissionsgrowth]

    # feed counterfactual GDP (for level effects) and actual GDP into emissions components to re-scale scenario emissions
    connect_param!(m, :co2emissions => :gdp_leveleffect, :GDP => :gdp_leveleffect)
    connect_param!(m, :co2emissions => :gdp, :GDP => :gdp)

    connect_param!(m, :CO2Cycle => :e_globalCO2emissions, :co2emissions => :e_globalCO2emissions)
    connect_param!(m, :CO2Cycle => :rt_g_globaltemperature, :ClimateTemperature => :rt_g_globaltemperature)
    if use_permafrost
        co2cycle[:permte_permafrostemissions] = permafrost[:perm_tot_e_co2]
    end

    connect_param!(m, :co2forcing => :c_CO2concentration, :CO2Cycle => :c_CO2concentration)

    ch4emit[:er_CH4emissionsgrowth] = scenario[:er_CH4emissionsgrowth]

    connect_param!(m, :ch4emissions => :gdp_leveleffect, :GDP => :gdp_leveleffect)
    connect_param!(m, :ch4emissions => :gdp, :GDP => :gdp)

    connect_param!(m, :CH4Cycle => :e_globalCH4emissions, :ch4emissions => :e_globalCH4emissions)
    connect_param!(m, :CH4Cycle => :rtl_g0_baselandtemp, :ClimateTemperature => :rtl_g0_baselandtemp)
    connect_param!(m, :CH4Cycle => :rtl_g_landtemperature, :ClimateTemperature => :rtl_g_landtemperature)
    if use_permafrost
        ch4cycle[:permtce_permafrostemissions] = permafrost[:perm_tot_ce_ch4]
    end

    connect_param!(m, :ch4forcing => :c_CH4concentration, :CH4Cycle => :c_CH4concentration)
    connect_param!(m, :ch4forcing => :c_N2Oconcentration, :n2ocycle => :c_N2Oconcentration)

    n2oemit[:er_N2Oemissionsgrowth] = scenario[:er_N2Oemissionsgrowth]

    connect_param!(m, :n2oemissions => :gdp_leveleffect, :GDP => :gdp_leveleffect)
    connect_param!(m, :n2oemissions => :gdp, :GDP => :gdp)

    connect_param!(m, :n2ocycle => :e_globalN2Oemissions, :n2oemissions => :e_globalN2Oemissions)
    connect_param!(m, :n2ocycle => :rtl_g0_baselandtemp, :ClimateTemperature => :rtl_g0_baselandtemp)
    connect_param!(m, :n2ocycle => :rtl_g_landtemperature, :ClimateTemperature => :rtl_g_landtemperature)

    connect_param!(m, :n2oforcing => :c_CH4concentration, :CH4Cycle => :c_CH4concentration)
    connect_param!(m, :n2oforcing => :c_N2Oconcentration, :n2ocycle => :c_N2Oconcentration)

    lgemit[:er_LGemissionsgrowth] = scenario[:er_LGemissionsgrowth]

    connect_param!(m, :LGemissions => :gdp_leveleffect, :GDP => :gdp_leveleffect)
    connect_param!(m, :LGemissions => :gdp, :GDP => :gdp)

    connect_param!(m, :LGcycle => :e_globalLGemissions, :LGemissions => :e_globalLGemissions)
    connect_param!(m, :LGcycle => :rtl_g0_baselandtemp, :ClimateTemperature => :rtl_g0_baselandtemp)
    connect_param!(m, :LGcycle => :rtl_g_landtemperature, :ClimateTemperature => :rtl_g_landtemperature)

    connect_param!(m, :LGforcing => :c_LGconcentration, :LGcycle => :c_LGconcentration)

    sulfemit[:pse_sulphatevsbase] = scenario[:pse_sulphatevsbase]

    connect_param!(m, :TotalForcing => :f_CO2forcing, :co2forcing => :f_CO2forcing)
    connect_param!(m, :TotalForcing => :f_CH4forcing, :ch4forcing => :f_CH4forcing)
    connect_param!(m, :TotalForcing => :f_N2Oforcing, :n2oforcing => :f_N2Oforcing)
    connect_param!(m, :TotalForcing => :f_lineargasforcing, :LGforcing => :f_LGforcing)
    totalforcing[:exf_excessforcing] = scenario[:exf_excessforcing]
    connect_param!(m, :TotalForcing => :fs_sulfateforcing, :SulphateForcing => :fs_sulphateforcing)

    connect_param!(m, :SeaLevelRise => :rt_g_globaltemperature, :ClimateTemperature => :rt_g_globaltemperature)

    if use_page09damages
        connect_param!(m, :GDP => :isat_ImpactinclSaturationandAdaptation, :MarketDamages => :isat_ImpactinclSaturationandAdaptation)
    else
        connect_param!(m, :GDP => :isat_ImpactinclSaturationandAdaptation, :MarketDamagesBurke => :isat_ImpactinclSaturationandAdaptation)
    end

    for allabatement in [
        (:AbatementCostParametersCO2, :AbatementCostsCO2, :er_CO2emissionsgrowth),
        (:AbatementCostParametersCH4, :AbatementCostsCH4, :er_CH4emissionsgrowth),
        (:AbatementCostParametersN2O, :AbatementCostsN2O, :er_N2Oemissionsgrowth),
        (:AbatementCostParametersLin, :AbatementCostsLin, :er_LGemissionsgrowth)]

        abatementcostparameters, abatementcosts, er_parameter = allabatement

        connect_param!(m, abatementcostparameters => :yagg, :GDP => :yagg_periodspan)
        connect_param!(m, abatementcostparameters => :cbe_absoluteemissionreductions, abatementcosts => :cbe_absoluteemissionreductions)

        connect_param!(m, abatementcosts => :zc_zerocostemissions, abatementcostparameters => :zc_zerocostemissions)
        connect_param!(m, abatementcosts => :q0_absolutecutbacksatnegativecost, abatementcostparameters => :q0_absolutecutbacksatnegativecost)
        connect_param!(m, abatementcosts => :blo, abatementcostparameters => :blo)
        connect_param!(m, abatementcosts => :alo, abatementcostparameters => :alo)
        connect_param!(m, abatementcosts => :bhi, abatementcostparameters => :bhi)
        connect_param!(m, abatementcosts => :ahi, abatementcostparameters => :ahi)
        connect_param!(m, abatementcosts => :er_emissionsgrowth, :RCPSSPScenario => er_parameter)

    end

    connect_param!(m, :TotalAbatementCosts => :tc_totalcosts_co2, :AbatementCostsCO2 => :tc_totalcost)
    connect_param!(m, :TotalAbatementCosts => :tc_totalcosts_n2o, :AbatementCostsN2O => :tc_totalcost)
    connect_param!(m, :TotalAbatementCosts => :tc_totalcosts_ch4, :AbatementCostsCH4 => :tc_totalcost)
    connect_param!(m, :TotalAbatementCosts => :tc_totalcosts_linear, :AbatementCostsLin => :tc_totalcost)
    connect_param!(m, :TotalAbatementCosts => :pop_population, :Population => :pop_population)

    connect_param!(m, :AdaptiveCostsEconomic => :gdp, :GDP => :gdp)
    connect_param!(m, :AdaptiveCostsNonEconomic => :gdp, :GDP => :gdp)
    connect_param!(m, :AdaptiveCostsSeaLevel => :gdp, :GDP => :gdp)

    connect_param!(m, :TotalAdaptationCosts => :ac_adaptationcosts_economic, :AdaptiveCostsEconomic => :ac_adaptivecosts)
    connect_param!(m, :TotalAdaptationCosts => :ac_adaptationcosts_noneconomic, :AdaptiveCostsNonEconomic => :ac_adaptivecosts)
    connect_param!(m, :TotalAdaptationCosts => :ac_adaptationcosts_sealevelrise, :AdaptiveCostsSeaLevel => :ac_adaptivecosts)
    connect_param!(m, :TotalAdaptationCosts => :pop_population, :Population => :pop_population)

    connect_param!(m, :SLRDamages => :s_sealevel, :SeaLevelRise => :s_sealevel)
    connect_param!(m, :SLRDamages => :cons_percap_consumption, :GDP => :cons_percap_consumption)
    connect_param!(m, :SLRDamages => :cons_percap_consumption_0, :GDP => :cons_percap_consumption_0)
    connect_param!(m, :SLRDamages => :tct_per_cap_totalcostspercap, :TotalAbatementCosts => :tct_per_cap_totalcostspercap)
    connect_param!(m, :SLRDamages => :act_percap_adaptationcosts, :TotalAdaptationCosts => :act_percap_adaptationcosts)
    connect_param!(m, :SLRDamages => :atl_adjustedtolerablelevelofsealevelrise, :AdaptiveCostsSeaLevel => :atl_adjustedtolerablelevel, ignoreunits=true)
    connect_param!(m, :SLRDamages => :imp_actualreductionSLR, :AdaptiveCostsSeaLevel => :imp_adaptedimpacts)
    connect_param!(m, :SLRDamages => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)
    ###############################################
    # Growth Effects - additional variables and parameters
    ###############################################
    connect_param!(m, :SLRDamages => :cons_percap_consumption_noconvergence, :GDP => :cons_percap_consumption_noconvergence)
    connect_param!(m, :SLRDamages => :cbabsn_pcconsumptionbound_neighbourhood, :GDP => :cbabsn_pcconsumptionbound_neighbourhood)
    connect_param!(m, :SLRDamages => :cbaux1_pcconsumptionbound_auxiliary1, :GDP => :cbaux1_pcconsumptionbound_auxiliary1)
    connect_param!(m, :SLRDamages => :cbaux2_pcconsumptionbound_auxiliary2, :GDP => :cbaux2_pcconsumptionbound_auxiliary2)
    connect_param!(m, :SLRDamages => :cons_percap_consumption_noconvergence, :GDP => :cons_percap_consumption_noconvergence)
    ###############################################

    connect_param!(m, :MarketDamages => :rtl_realizedtemperature, :ClimateTemperature => :rtl_realizedtemperature)
    connect_param!(m, :MarketDamages => :rgdp_per_cap_SLRRemainGDP, :SLRDamages => :rgdp_per_cap_SLRRemainGDP)
    connect_param!(m, :MarketDamages => :rcons_per_cap_SLRRemainConsumption, :SLRDamages => :rcons_per_cap_SLRRemainConsumption)
    connect_param!(m, :MarketDamages => :atl_adjustedtolerableleveloftemprise, :AdaptiveCostsEconomic => :atl_adjustedtolerablelevel, ignoreunits=true) # not required for Burke damages
    connect_param!(m, :MarketDamages => :imp_actualreduction, :AdaptiveCostsEconomic => :imp_adaptedimpacts) # not required for Burke damages
    connect_param!(m, :MarketDamages => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)

    connect_param!(m, :MarketDamagesBurke => :rtl_realizedtemperature, :ClimateTemperature => :rtl_realizedtemperature)
    connect_param!(m, :MarketDamagesBurke => :rgdp_per_cap_SLRRemainGDP, :SLRDamages => :rgdp_per_cap_SLRRemainGDP)
    connect_param!(m, :MarketDamagesBurke => :rcons_per_cap_SLRRemainConsumption, :SLRDamages => :rcons_per_cap_SLRRemainConsumption)
    connect_param!(m, :MarketDamagesBurke => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)

    connect_param!(m, :NonMarketDamages => :rtl_realizedtemperature, :ClimateTemperature => :rtl_realizedtemperature)
    if use_page09damages
        connect_param!(m, :NonMarketDamages => :rgdp_per_cap_MarketRemainGDP, :MarketDamages => :rgdp_per_cap_MarketRemainGDP)
        connect_param!(m, :NonMarketDamages => :rcons_per_cap_MarketRemainConsumption, :MarketDamages => :rcons_per_cap_MarketRemainConsumption)
    else
        connect_param!(m, :NonMarketDamages => :rgdp_per_cap_MarketRemainGDP, :MarketDamagesBurke => :rgdp_per_cap_MarketRemainGDP)
        connect_param!(m, :NonMarketDamages => :rcons_per_cap_MarketRemainConsumption, :MarketDamagesBurke => :rcons_per_cap_MarketRemainConsumption)
    end
    connect_param!(m, :NonMarketDamages => :atl_adjustedtolerableleveloftemprise, :AdaptiveCostsNonEconomic => :atl_adjustedtolerablelevel, ignoreunits=true)
    connect_param!(m, :NonMarketDamages => :imp_actualreduction, :AdaptiveCostsNonEconomic => :imp_adaptedimpacts)
    connect_param!(m, :NonMarketDamages => :isatg_impactfxnsaturation, :GDP => :isatg_impactfxnsaturation)

    connect_param!(m, :Discontinuity => :rgdp_per_cap_NonMarketRemainGDP, :NonMarketDamages => :rgdp_per_cap_NonMarketRemainGDP)
    connect_param!(m, :Discontinuity => :rt_g_globaltemperature, :ClimateTemperature => :rt_g_globaltemperature)
    connect_param!(m, :Discontinuity => :rgdp_per_cap_NonMarketRemainGDP, :NonMarketDamages => :rgdp_per_cap_NonMarketRemainGDP)
    connect_param!(m, :Discontinuity => :rcons_per_cap_NonMarketRemainConsumption, :NonMarketDamages => :rcons_per_cap_NonMarketRemainConsumption)
    connect_param!(m, :Discontinuity => :isatg_saturationmodification, :GDP => :isatg_impactfxnsaturation)

    connect_param!(m, :EquityWeighting => :pop_population, :Population => :pop_population)
    connect_param!(m, :EquityWeighting => :tct_percap_totalcosts_total, :TotalAbatementCosts => :tct_per_cap_totalcostspercap)
    connect_param!(m, :EquityWeighting => :act_adaptationcosts_total, :TotalAdaptationCosts => :act_adaptationcosts_total)
    connect_param!(m, :EquityWeighting => :act_percap_adaptationcosts, :TotalAdaptationCosts => :act_percap_adaptationcosts)
    connect_param!(m, :EquityWeighting => :cons_percap_consumption, :GDP => :cons_percap_consumption)
    connect_param!(m, :EquityWeighting => :cons_percap_consumption_0, :GDP => :cons_percap_consumption_0)
    connect_param!(m, :EquityWeighting => :cons_percap_aftercosts, :SLRDamages => :cons_percap_aftercosts)
    connect_param!(m, :EquityWeighting => :rcons_percap_dis, :Discontinuity => :rcons_per_cap_DiscRemainConsumption)
    connect_param!(m, :EquityWeighting => :yagg_periodspan, :GDP => :yagg_periodspan)
    ###############################################
    # Growth Effects - additional variables and parameters
    ###############################################
    connect_param!(m, :EquityWeighting => :grwnet_realizedgdpgrowth, :GDP => :grwnet_realizedgdpgrowth)
    connect_param!(m, :EquityWeighting => :lgdp_gdploss, :GDP => :lgdp_gdploss)
    ###############################################
    equityweighting[:grw_gdpgrowthrate] = scenario[:grw_gdpgrowthrate]
    equityweighting[:popgrw_populationgrowth] = scenario[:popgrw_populationgrowth]

    return m
end

function initpage(m::Model)
    p = load_parameters(m)
    p["y_year_0"] = 2015.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)
end

function getpage(scenario::String="RCP4.5 & SSP2", use_permafrost::Bool=true, use_seaice::Bool=true, use_page09damages::Bool=false)
    m = Model()
    set_dimension!(m, :time, [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300])
    set_dimension!(m, :region, ["EU", "USA", "OECD","USSR","China","SEAsia","Africa","LatAmerica"])
    set_dimension!(m, :draw, Array(1:5000))

    buildpage(m, scenario, use_permafrost, use_seaice, use_page09damages)

    # next: add vector and panel example
    initpage(m)

    return m
end

get_model = getpage


text:
function build_LOModel(data::Dict{String, Any})

    lom = LaplacianOptModel(data)

    LOpt.variable_LOModel(lom)

    if lom.data["solution_type"] == "exact"
        LOpt.constraint_LOModel(lom)
    end

    LOpt.objective_LOModel(lom)

    return lom
end

function variable_LOModel(lom::LaplacianOptModel)

    LOpt.variable_lifted_W_matrix(lom)
    LOpt.variable_edge_onoff(lom)
    LOpt.variable_algebraic_connectivity(lom)

    return
end

function constraint_LOModel(lom::LaplacianOptModel)

    LOpt.constraint_build_W_var_matrix(lom)

    LOpt.constraint_topology_no_self_loops(lom)
    LOpt.constraint_topology_vertex_cutset(lom)
    LOpt.constraint_topology_total_edges(lom)

    LOpt.constraint_lazycallback_wrapper(lom)

    return
end

function objective_LOModel(lom::LaplacianOptModel)
    LOpt.objective_maximize_algebraic_connectivity(lom)

    return
end

function optimize_LOModel!(lom::LaplacianOptModel; optimizer=nothing)
    if lom.data["relax_integrality"]
        JuMP.relax_integrality(lom.model)
    end

    if JuMP.mode(lom.model) != JuMP.DIRECT && optimizer !== nothing
        if lom.model.moi_backend.state == MOI.Utilities.NO_OPTIMIZER
            JuMP.set_optimizer(lom.model, optimizer)
        else
            Memento.warn(_LOGGER, "Model already contains optimizer, cannot use optimizer specified in `optimize_LOModel!`")
        end
    end

    if JuMP.mode(lom.model) != JuMP.DIRECT && lom.model.moi_backend.state == MOI.Utilities.NO_OPTIMIZER
        Memento.error(_LOGGER, "No optimizer specified in `optimize_LOModel!` or the given JuMP model.")
    end

    start_time = time()

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

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

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

    lom.result = LOpt.build_LOModel_result(lom, solve_time)

    return lom.result
end

function run_LOpt_model(params::Dict{String, Any}, lom_optimizer::MOI.OptimizerWithAttributes; visualize_solution = false, visualizing_tool = "graphviz")

    data = LOpt.get_data(params)

    model_lopt  = LOpt.build_LOModel(data)

    result_lopt = LOpt.optimize_LOModel!(model_lopt, optimizer = lom_optimizer)

    if visualize_solution
        LOpt.visualize_solution(result_lopt, data, visualizing_tool = visualizing_tool)
    end

    return result_lopt
end

function run_MaxSpanTree_model(params::Dict{String, Any}, lom_optimizer::MOI.OptimizerWithAttributes; visualize_solution = false, visualizing_tool = "graphviz", lazy_callback = false)
    data = LOpt.get_data(params)

    model_mst = LOpt.build_MaxSpanTree_model(data, lazy_callback)

    result_mst = LOpt.optimize_LOModel!(model_mst, optimizer = lom_optimizer)

    if visualize_solution
        LOpt.visualize_solution(result_mst, data, visualizing_tool = visualizing_tool)
    end

    return result_mst

end

function build_MaxSpanTree_model(data::Dict{String, Any}, lazy_callback::Bool)

    m_mst = LaplacianOptModel(data)

    LOpt.variable_MaxSpanTree_model(m_mst, lazy_callback)

    LOpt.constraint_MaxSpanTree_model(m_mst, lazy_callback)

    LOpt.objective_MaxSpanTree_model(m_mst)

    return m_mst
end

function objective_MaxSpanTree_model(lom::LaplacianOptModel)

    LOpt.objective_maximize_spanning_tree_cost(lom)

    return
end

function variable_MaxSpanTree_model(lom::LaplacianOptModel, lazy_callback::Bool)

    LOpt.variable_edge_onoff(lom)

    if !lazy_callback
        LOpt.variable_multi_commodity_flow(lom)
    end

    return
end

function constraint_MaxSpanTree_model(lom::LaplacianOptModel, lazy_callback::Bool)

    LOpt.constraint_topology_no_self_loops(lom)
    LOpt.constraint_topology_vertex_cutset(lom)
    LOpt.constraint_topology_total_edges(lom)

    if lazy_callback
        LOpt.constraint_lazycallback_wrapper(lom)
    else
        LOpt.constraint_topology_multi_commodity_flow(lom)
    end

    return
end


text:
# Constants

const baseDir = "/Users/lutz/Documents/projects/p2018/migrant_wage_gains";

# Filter criteria
const birthYearMin = 1900
const birthYearMax = 1990

# Country indicators
const homeCountry = 0
const unknownCountry = -1

# --------------


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.SimulationEngine(ModiaMath/SimulationEngine/_module.jl)
#
import DASSL
using DASSL

mutable struct IntegratorData
    model::ModiaMath.AbstractSimulationModel    # Model specific data structure
    simulationState::ModiaMath.SimulationState
    statistics::ModiaMath.SimulationStatistics   # Simulation statistics
    z::Vector{Float64}          # Vector of event indicators (zero crossings). If one of z[i] passes
                                # zero, that is beforeEvent(z[i])*z[i] < 0, an event is triggered
                                # (allocated during instanciation according to nz).
    zDir::Vector{Int32}          # zDir[i] =  0: Root is reported for both crossing directions
                                #         =  1: Root is reported when crossing from negative to positive direction
    eweight::Vector{Float64}    # Weight vector
    last_t::Float64              # Time instant of the last call to idasol_f (residue function)
    last_norm_y::Float64         # Norm of y of the last call to idasol_f (residue function)
    hcur::MVector{1,Float64}     # Current step size argument for IDAGetCurrentStep
    order::MVector{1,Int32}       # Current order for IDAGetCurrentOrder
    r::Vector{Float64} # Residues vector used in sol_f function

    y::Vector{Float64}    # Julia vector of y wrapping y_N_Vector vector
    yp::Vector{Float64}   # Julia vector of yp wrapping yp_N_Vector vector
    fulljac::Union{Matrix{Float64}, Nothing}  # Julia matrix of fulljac wrapping fulljac_N_xxx

    function IntegratorData(model::ModiaMath.AbstractSimulationModel, simulationState::ModiaMath.SimulationState)
        ny = simulationState.nx
        nz = simulationState.nz
        new(model, simulationState, simulationState.statistics,
          ones(nz), fill(0, 0), zeros(ny), typemin(Float64), 0.0,
          MVector{1,Float64}(0.0), MVector{1,Int32}(0), zeros(ny))
    end
end

include("UserData.jl")
#=
function updateStatistics!(sol_mem, stat::ModiaMath.SimulationStatistics)
    info = 0
    info = sol_mem.steps
    stat.nSteps        += info

    info = sol_mem.jac_evals
    stat.nJac          += info

    info = sol_mem.fails
    stat.nErrTestFails += info

end
=#

"""
    result = closeTimePoints(t1,t2,epsilon)

The function returns true if the two time points t1 and t2 are close together (with respect to epsilon)
"""
closeTimePoints(t1::Float64, t2::Float64, epsilon::Float64) = abs(t1 - t2) / max(abs(t2), 10 * epsilon) <= 100 * epsilon

"""
    ModiaMath.simulate!(simulationModel; log=false, startTime=NaN, stopTime=NaN,
                                                    tolerance=NaN, interval=NaN)

Simulates a DAE `simulationModel` that is defined with package `Modia`, package `Modia3D` or
with the `ModiaMath.@component` macro. The DAE is mathematically described as
implicit, index 1 DAE with events (containing an ODE or a semi-explicit
index 1 DAE as special cases). For details see [`ModiaMath.StructureOfDAE`](@ref).

During continuous integration, the DAE is solved with the
[Sundials](https://computation.llnl.gov/projects/sundials) IDA solver
(accessed via the Julia [Sundials](https://github.com/JuliaDiffEq/Sundials.jl) interface package).

Input arguments

- `simulationModel::ModiaMath.AbstractSimulationModel`: Model struct (generated with Modia, Modia3D or `ModiaMath.@component`).
- `log::Bool`: = true, if logging is enabled, otherwise it is disabled.
- `startTime::Float64`: Start time of the simulation in [s].
                        If startTime=NaN, the default startTime is used that is defined by the `simulationModel`.
- `stopTime::Float64`: Stop time of the simulation in [s].
                       If stopTime=NaN, the default stopTime is used that is defined by the `simulationModel`.
- `tolerance::Float64`: The relative tolerance for the integration.
  The absolute tolerance is computed as `0.1*tolerance*nominal(variable)` where
  `nominal(variable)` is the nominal value of the variable.
   If tolerance=NaN, the default tolerance is used that is defined by the `simulationModel`.
- `interval::Float64`: Output interval for results in [s]. If events occur, the event time instants are
  additionally added to the result.
  If interval=NaN, the default interval is used that is defined by the `simulationModel`.
"""
function simulate!(model::ModiaMath.AbstractSimulationModel;
                   tolerance=NaN,
                   startTime=NaN,
                   stopTime=NaN,
                   interval=NaN,
                   log::Bool=false,
                   KLUorderingChoice::Int=1)

    use_fulljac = false

    sim           = model.simulationState
    sim.model     = model
    sim.tolerance = isnan(tolerance) ? sim.defaultTolerance : convert(Float64, tolerance)
    sim.startTime = isnan(startTime) ? sim.defaultStartTime : convert(Float64, startTime)
    sim.stopTime  = isnan(stopTime)  ? sim.defaultStopTime  : convert(Float64, stopTime)
    sim.interval  = isnan(interval)  ? (isnan(sim.defaultInterval) ? (sim.stopTime - sim.startTime)/500.0 : sim.defaultInterval)  : convert(Float64, interval)
    sim.time      = sim.startTime
    @assert(sim.stopTime >= sim.startTime)
    @assert(sim.interval >= 0.0)
    @assert(sim.tolerance > 0.0)

    stopTime2::Float64  = sim.stopTime
    startTime2::Float64 = sim.startTime
    tolerance::Float64  = sim.tolerance
    interval::Float64   = sim.interval
    logger              = sim.logger
    ModiaMath.setLog!(logger, log)

    if ModiaMath.isLogInfos(logger)
        println("... ModiaMath.simulate! (version ", ModiaMath.Version, " ", ModiaMath.Date, ") to simulate model: ", sim.name)
    end

    # Start timing measure
    cpuStart::UInt64 = time_ns()
    cpuLast::UInt64  = cpuStart
    cpuStartIntegration::UInt64 = cpuStart

    # Init variables
    nt::Int           = round(Int64, 1 + (stopTime2 - startTime2) / interval)
    t0::Float64       = startTime2
    tEnd::Float64     = stopTime2
    tReached::Float64 = t0
    epsilon::Float64  = eps()

    # Initialize simulation model and store result after initialization
    statistics = sim.statistics
    ModiaMath.reInitializeStatistics!(statistics, t0, stopTime2, interval, tolerance)
    statistics.h0 = interval
    statistics.hMin = Inf
    statistics.hMax = 0
    if ModiaMath.isLogInfos(logger)
        println("      Initialization at time = ", t0, " s")
    end

    init = DAE.initialize!(model, sim, t0, nt, tolerance)
    @assert(length(init.y0)  == sim.nx)
    @assert(length(init.yp0) == sim.nx)
    @assert(tolerance > 0.0)

    if init.terminate || tEnd <= t0
        statistics.h0   = 0.0
        statistics.hMin = 0.0
        @goto TerminateSimulation
    end

    if ModiaMath.isLogInfos(logger)
        println("      Simulation started")
    end

    maxTime::Float64          = init.maxTime
    nextEventTime::Float64    = init.nextEventTime
    integrateToEvent::Bool = init.integrateToEvent

    if nextEventTime < t0
        nextEventTime = typemax(Float64)
    end

    # Initialize auxiliary variables
    simModel = IntegratorData(model, sim)
    simModel.hcur[1] = statistics.h0
    stateEvent::Bool = false
    timeEvent::Bool  = false
    flag = false
    ny = sim.nx
    nz = sim.nz

    y  = copy(init.y0)
    yp = copy(init.yp0)
    simModel.y  = y
    simModel.yp = yp
    tolAbs = tolerance * 5e-4 #magic number passes tests
    if use_fulljac
        fulljac = zeros(sim.nx, sim.nx)
        simModel.fulljac = fulljac
    else
        simModel.fulljac = nothing
    end

    eventInfo = DAE.EventInfo()
    tret    = [0.0]
    jc::Int = 0  # communication point counter
    if ModiaMath.isLogInfos(logger) && nz > 0
        rootInfo = fill(0, nz)
    end

    #will be used later to let solver operate with vectors of more presice than Float64 reals
    y_Vector  = copy(y)
    yp_Vector = copy(yp)
    relTol =  0.0001 * tolerance * init.y_nominal
    for i in 1:ny
        if !init.y_errorControl[i]
            tolAbs[i] = 1e5 * init.y_nominal[i]   # switch tolerance control off
        end
    end
    # Run solution
    # Initialize solvers
    tReached = t0
    solver = dassl()
    #solver = IDA()

    y0 = y_Vector
    yp0 = yp_Vector
    if use_fulljac
        # lsJacFn
        jac_fun = sol_fulljac
    end

    # Sundials.IDASetId(mem, y_states_nvector)
    #if sim.sparse
    #   nnz_jac = nnz(sim.jac)
    #   Sundials.IDAKLU(mem, ny, nnz_jac)
    #   Sundials.IDAKLUSetOrdering(mem, KLUorderingChoice)
    #   Sundials.IDASlsSetSparseJacFn(mem);
    #else
    #A  = Sundials.SUNDenseMatrix(ny, ny)
    #LS = Sundials.SUNDenseLinearSolver(y, A)
    #Sundials.IDADlsSetLinearSolver(mem, LS, A)
    #end

    # Initialize zero crossing function, if required
    hasZeroCross = nz > 0
    if hasZeroCross
        ModiaMath.DAE.getEventIndicators!(simModel.model, sim, t0, y0, yp0, simModel.z)
        simModel.zDir = zeros(nz)
        sim.zDir = simModel.zDir
    end

    # Initialize event variables
    tStop::Float64 = tEnd
    tNext = t0
    tc    = t0
    nTimeEvents = 0
    restart     = ModiaMath.NoRestart

    # Check for event at time = t0. If yes, trigger event and store result after event
    if closeTimePoints(t0, nextEventTime, epsilon)
        tReached = t0
        if ModiaMath.isLogEvents(logger)
            println("\n      Time event at time = $tReached s")
        end

        DAE.reset!(eventInfo)
        DAE.processEvent!(model, sim, tReached, y_Vector, yp_Vector, eventInfo)
        restart          = eventInfo.restart
        maxTime          = eventInfo.maxTime
        nextEventTime    = eventInfo.nextEventTime
        integrateToEvent = eventInfo.integrateToEvent

        statistics.nTimeEvents += 1
        if ModiaMath.isLogEvents(logger)
            println("        restart = ", restart)
        end

        # Reinitialize DAE if required by model
        if restart == ModiaMath.Restart
            statistics.nRestartEvents += 1
            simModel.y = y
            simModel.yp = yp
            DAE.reinitialize!(simModel.model, sim, tReached)
        end
    end

    # Loop over communication points
    cpuStartIntegration = time_ns()

    flag = false

    while tReached < tEnd && restart != ModiaMath.Terminate  #------------------ while --------------
        #println("entered")
        if tReached >= tc
            # Last communication point reached -> increment communication point
            jc += 1
            tc = startTime2 + jc * interval
        end

        if closeTimePoints(nextEventTime, tc, epsilon)
            nextEventTime = tc
        end

        if nextEventTime <= tc # next tStop is a time event
            event = true
            tNext = nextEventTime
            tStop = integrateToEvent ? nextEventTime : min(tEnd, maxTime)
        else        # next tStop is a communication point (no event)
            event = false
            tNext = tc
            tStop = min(tEnd, maxTime)
        end

        # Integrate in direction of tStop
        #println("... tReached = ", tReached, ", tNext = ", tNext, ", tStop = ", tStop,
        #        ", y = ", y[1], ", yp = ", yp[1])
        #initialize derivatives
        r =  simModel.r
        tspan = (tReached, tNext)
        resprob!(r, du, u, p, t) = sol_f!(simModel, sim, t, u, du, r, simModel.hcur[1])
        differential_vars = [true,true,false]
        prob = DAEProblem{true}(resprob!, yp, y, tspan, differential_vars=differential_vars)
        sol = Vector{Any}
        try #some problems require higher tolerance than passed, some fail with this hish tolerance
            sol = solve(prob, solver, reltol = relTol, abstol = tolAbs, initstep = simModel.hcur[1] )
        catch
            tAbs = tolAbs * 100
            rTol = relTol * 100
            sol = solve(prob, solver, reltol = rTol, abstol = tAbs, initstep = simModel.hcur[1] )
        end

        endInd = length(sol.t)
        simModel.yp = sol.du[end]
        if hasZeroCross

            len = length(sol.t)
            for el in 2:len
                ModiaMath.DAE.getEventIndicators!(simModel.model, sim, sol.t[el-1], sol.u[el-1], sol.du[el-1], simModel.z)
                old_z = copy(simModel.z)
                #println("old_z = $old_z")
                ModiaMath.DAE.getEventIndicators!(simModel.model, sim, sol.t[el], sol.u[el], sol.du[el], simModel.z)
                zs = old_z.*(simModel.z)
                z = copy(simModel.z)
                #println("new_z = $z, zs = $zs")
                #println("Event ind = $z")
                if (tReached!=t0)
                    flag = any(x->x<0, zs) && !flag
                    #println("flag = $flag")
                end
                if flag
                endInd = el
                #println("before event at all: y =", sol.u, ", \n yp = ", sol.du)
                    break
                end
            end
        end

        #update statistics
        nSteps = length(sol.t)-1
        if nSteps>0
            steps = zeros(nSteps)
            for it in 1:nSteps
                step = sol.t[it+1] - sol.t[it]
                steps[it] = step
            end
            simModel.hcur[1] = steps[end]
            statistics.h0 = steps[end]

            sim.time             = tReached
            simModel.last_t      = tReached
            simModel.last_norm_y = norm(simModel.y, Inf)
            #println("stat = ", statistics.h0)
            statistics.hMin = min(statistics.hMin, minimum(steps))
            statistics.hMax = max(statistics.hMax, maximum(steps))
        end

        if flag
            #binary search
            ind_l = endInd-1
            ind_r = endInd
            #println("indices are ind_l = $ind_l and ind_r = $ind_r")
            yMin = sol.u[ind_l]
            ypMin = sol.du[ind_l]
            tMin = sol.t[ind_l]
            yMax = sol.u[ind_r]
            ypMax = sol.du[ind_r]
            tMax = sol.t[ind_r]
            z = copy(simModel.z)
            while !any(x->abs(x)<tolAbs, z) #binsearch by values approx
            #while  false #binsearch by values approx
                tspan = (tMin, tMax)
                yMid = (yMax + yMin)/2
                ypMid = (ypMax + ypMin)/2
                tMid = (tMax + tMin)/2
                ModiaMath.DAE.getEventIndicators!(simModel.model, sim, tMid, yMid, ypMid, simModel.z)
                z = copy(simModel.z)
                if any(x->x<0, old_z.*(simModel.z))
                    yMax = yMid
                    ypMax = ypMid
                    tMax = tMid
                else
                    yMin = yMid
                    ypMin = ypMid
                    tMin = tMid
                end
                #println("y internal = $yMid, yp = $ypMid,  z = $z, t = $tMid")
            end  #end binsearch by values approx
            y = yMax
            tReached = tMax
            yp = ypMax
            simModel.yp = yp
        else #-- if not flag
            y = sol.u[end]
            yp = sol.du[end]
            tReached = sol.t[end]
        end #-- end if flag --

        #println("y new = $y, yp = $yp,  tol = $tolerance, t = $tReached")
        simModel.y = y
        simModel.yp = yp

        sim.time             = tReached
        simModel.last_t      = tReached
        simModel.last_norm_y = norm(simModel.y, Inf)

        stateEvent = flag
        timeEvent  = event && closeTimePoints(tNext, tReached, epsilon)
        isEvent    = timeEvent || stateEvent
        # Store result point at tReached
        DAE.computeAndStoreResult!(model, sim, tReached, y, yp)

        # If event at tReached, handle the event
        if isEvent
            if ModiaMath.isLogEvents(logger)
                if timeEvent && stateEvent
                    print("\n      Time and state (zero-crossing) event at time = $tReached s")
                elseif timeEvent
                    println("\n      Time event at time = $tReached s")
                elseif stateEvent
                    print("\n      State event (zero-crossing) at time = $tReached s")
                end
                #println("\n y = $y, yp=$yp")
                if stateEvent
                   # Print information about the root
                    ModiaMath.DAE.getEventIndicators!(simModel.model, sim, tReached, y, yp, simModel.z)
                    rootInfo = simModel.z
                    print(" (")
                    firstRoot = true
                    for iroot in eachindex(rootInfo)
                        if rootInfo[iroot] != 0
                            if firstRoot
                                firstRoot = false
                            else
                                print(", ")
                            end
                            print("z[", iroot, "]", rootInfo[iroot] > 0 ? " > 0" : " < 0")
                        end
                    end
                    print(")\n")
                end
            end
            old_z = simModel.z
            # Trigger event
            DAE.reset!(eventInfo)
            DAE.processEvent!(model, sim, tReached, y, yp, eventInfo)
            #println("\n  event = $eventInfo")
            #println("\n after event y = $y, yp=$yp, t = $tReached")

            restart          = eventInfo.restart
            maxTime          = eventInfo.maxTime
            nextEventTime    = eventInfo.nextEventTime
            integrateToEvent = eventInfo.integrateToEvent

            if timeEvent
                statistics.nTimeEvents += 1
            end

            if stateEvent
                statistics.nStateEvents += 1
            end

            if restart == ModiaMath.Restart || restart == ModiaMath.FullRestart
                statistics.nRestartEvents += 1
            end

            if ModiaMath.isLogEvents(logger)
                println("        restart = ", restart)
            end

            # Handle restart flag
            if restart == ModiaMath.Restart
                # Restart (dimensions do not change)
                simModel.y = y
                simModel.yp = yp
                DAE.reinitialize!(simModel.model, sim, tReached)
                #y = simModel.y
                #yp = simModel.yp
                ##Sundials.__IDAReInit(mem, tReached, y_N_Vector, yp_N_Vector);
            elseif restart == ModiaMath.FullRestart
                error("FullRestart not yet supported")
            end
        end

        # print every 5 s the time in order that the user sees the progress
        if ModiaMath.isLogProgress(logger)
            cpuNew = time_ns()
            if (cpuNew - cpuLast) * 1e-9 > 5.0
                cpuLast = cpuNew
                @printf("      progress: integrated up to time = %.2g s\n", tReached)
            end
        end
    end #---------------------- end while ------------------------------------

    # Finalize statistics for simulation run

    @label TerminateSimulation
    if ModiaMath.isLogInfos(logger)
        println("\n      Simulation is terminated at time = ", tEnd, " s")
    end

    result = DAE.terminate!(model, sim, tReached, y, yp)
    statistics.cpuTimeInitialization = (cpuStartIntegration - cpuStart) * 1e-9
    statistics.cpuTimeIntegration    = (time_ns() - cpuStartIntegration) * 1e-9

    if ModiaMath.isLogStatistics(logger)
        println("\n      Statistics (get help with ?ModiaMath.SimulationStatistics):")
        display(statistics)
    end
    return result
end


text:
using Test

@testset "Test for countZeros" begin include("test_countZeros.jl") end
@testset "Test for subdivision algorithm" begin include("test_subdivision.jl") end


text:
export PPOTrajectory

using MacroTools

struct PPOTrajectory{T<:CircularCompactSARTSATrajectory,P,names,types} <:
       AbstractTrajectory{names,types}
    trajectory::T
    action_log_prob::P
end

function PPOTrajectory(;
    capacity,
    action_log_prob_size = (),
    action_log_prob_type = Float32,
    kw...,
)
    t = CircularCompactSARTSATrajectory(; capacity = capacity, kw...)
    p = CircularArrayBuffer{action_log_prob_type}(action_log_prob_size..., capacity + 1)
    names = typeof(t).parameters[1]
    types = typeof(t).parameters[2]
    PPOTrajectory{
        typeof(t),
        typeof(p),
        (
            :state,
            :action,
            :action_log_prob,
            :reward,
            :terminal,
            :next_state,
            :next_action,
            :next_action_log_prob,
        ),
        Tuple{
            types.parameters[1:2]...,
            frame_type(p),
            types.parameters[3:end]...,
            frame_type(p),
        },
    }(
        t,
        p,
    )
end

MacroTools.@forward PPOTrajectory.trajectory Base.length, Base.isempty, RLCore.isfull

function RLCore.get_trace(t::PPOTrajectory, s::Symbol)
    if s == :action_log_prob
        select_last_dim(
            t.action_log_prob,
            1:(nframes(t.action_log_prob) > 1 ? nframes(t.action_log_prob) - 1 :
               nframes(t.action_log_prob)),
        )
    elseif s == :next_action_log_prob
        select_last_dim(t.action_log_prob, 2:nframes(t.action_log_prob))
    else
        get_trace(t.trajectory, s)
    end
end

Base.getindex(t::PPOTrajectory, s::Symbol) =
    s == :action_log_prob ? t.action_log_prob : t.trajectory[s]

function Base.getindex(p::PPOTrajectory, i::Int)
    s, a, r, t, s′, a′ = p.trajectory[i]
    (
        state = s,
        action = a,
        action_log_prob = select_last_dim(p.action_log_prob, i),
        reward = r,
        terminal = t,
        next_state = s′,
        next_action = a′,
        next_action_log_prob = select_last_dim(p.action_log_prob, i + 1),
    )
end

function Base.empty!(b::PPOTrajectory)
    empty!(b.action_log_prob)
    empty!(b.trajectory)
end

function Base.push!(b::PPOTrajectory, kv::Pair{Symbol})
    k, v = kv
    if k == :action_log_prob || k == :next_action_log_prob
        push!(b.action_log_prob, v)
    else
        push!(b.trajectory, kv)
    end
end

function Base.pop!(t::PPOTrajectory, s::Symbol)
    if s == :action_log_prob || s == :next_action_log_prob
        pop!(t.action_log_prob)
    else
        pop!(t.trajectory, s)
    end
end

function Base.pop!(t::PPOTrajectory)
    (pop!(t.trajectory)..., action_log_prob = pop!(t.action_log_prob))
end


text:
module Defs

export @vassert
export vifthen

macro vassert(cond)
end

@inline function vifthen{T}(cond::Bool, x::T, y::T)
    return cond ? x : y
end

end


text:
export get_info_from_DB
import Random

"""
    function get_info_from_DB(db, horns, stokes)

This function can be used to take information about specific horns
from the STRIP instrument database.

Given:
- The database (db = Stripeline.InstrumentDB())
- A vector containg the names of the horns of interest (e.g. ["Y0",
  "I1", "R6"])
- A string indicating the Stokes parameter object of the simulation
  (`"Q"` or `"U"`)

It returns 5 arrays, containing for each horn:
- the 3D vector containing the orientation of the horn in the sky
- the unique ID of the polarimeter associated to the horn
- the estimated bandwidth, in Hz
- the estimated noise temperature, in K
- the estimated knee frequency, in Hz

For those polarimeters without valid measurements of the knee
frequency, the function randomly picks a value from a vector
containing all the fknee measurements available in the database.

"""
function get_info_from_DB(db, horns, stokes)

    @assert stokes in ["Q", "U"]

    orientations=[]
    polarimeters = Array{Int64}(undef, length(horns))
    β_hz = Array{Float64}(undef,length(horns))
    tnoise_k = Array{Float64}(undef,length(horns))
    fknee_hz = Array{Float64}(undef,length(horns))
    slopes = Array{Float64}(undef,length(horns))

    measured_fknees = Float64[]
    measured_slopes = Float64[]

    for i in 1:length(horns)
        polarimeters[i] = db.focalplane[horns[i]].polid
        if(db.detectors[polarimeters[i]].spectrum.fknee_q_hz !=0)
            if stokes == "Q"
                fknee = db.detectors[polarimeters[i]].spectrum.fknee_q_hz
                slope = db.detectors[polarimeters[i]].spectrum.slope_q
            else stokes == "U"
                fknee = db.detectors[polarimeters[i]].spectrum.fknee_u_hz
                slope = db.detectors[polarimeters[i]].spectrum.slope_u
            end
            push!(measured_fknees, fknee)
            push!(measured_slopes, clamp(slope, 0.0, 2.0))

        end
    end
    rng = Random.MersenneTwister(1234)

    for i in 1:length(horns)
        append!(orientations, [db.focalplane[horns[i]].orientation])
        β_hz[i] = db.detectors[polarimeters[i]].bandshape.bandwidth_hz
        tnoise_k[i] = db.detectors[polarimeters[i]].tnoise.tnoise_k

        if(db.detectors[polarimeters[i]].spectrum.fknee_q_hz ==0)
            fknee_hz[i] = rand(rng, measured_fknees)
            slopes[i] = rand(rng, measured_slopes)
        else
            if stokes == "Q"
                fknee_hz[i] = db.detectors[polarimeters[i]].spectrum.fknee_q_hz
                slope = db.detectors[polarimeters[i]].spectrum.slope_q
            else
                fknee_hz[i] = db.detectors[polarimeters[i]].spectrum.fknee_u_hz
                slope = db.detectors[polarimeters[i]].spectrum.slope_u
            end
            #correct for too high slopes, since noise generator works till slope=2
            if slope > 2
                slopes[i] = 2.0
            else
                slopes[i] = slope
            end
        end
    end
    return orientations, polarimeters, β_hz, tnoise_k, fknee_hz, slopes
end


text:
@testset "MinimumAmong2MILP: $(fct_type), dimension $(array_dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], array_dim in [2, 3], T in [Int, Float64]
    dim = 1 + array_dim
    mock = MOIU.MockOptimizer(MILPModel{T}())
    model = COIB.MinimumAmong2MILP{T}(mock)

    if T == Int
        @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    end
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.MinimumAmong,
    )

    if T == Int
        x, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:dim])
    elseif T == Float64
        x = MOI.add_variables(model, dim)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x)
    else
        @assert false
    end

    @test_throws AssertionError MOI.add_constraint(model, fct, CP.MinimumAmong(array_dim))

    for i in 1:array_dim
        MOI.add_constraint(model, x[1 + i], MOI.GreaterThan(zero(T)))
        MOI.add_constraint(model, x[1 + i], MOI.LessThan(one(T)))
    end
    c = MOI.add_constraint(model, fct, CP.MinimumAmong(array_dim))

    for i in 1:dim
        @test MOI.is_valid(model, x[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.MinimumAmong}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.MinimumAmong) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,)]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == array_dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == array_dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == array_dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == array_dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == 1

        @test MOI.get(bridge, MOI.ListOfVariableIndices()) == bridge.vars
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}())) == Set(bridge.vars_bin)
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}())) == Set(bridge.cons_gt)
        @test Set(MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}())) == Set(bridge.cons_lt)
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == [bridge.con_choose_one]
    end

    @testset "Binary variables" begin
        @test length(bridge.vars) == array_dim
        for i in 1:array_dim
            @test MOI.is_valid(model, bridge.vars[i])
        end

        @test MOI.is_valid(model, bridge.con_choose_one)
        s = MOI.get(model, MOI.ConstraintSet(), bridge.con_choose_one)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_choose_one)

        @test typeof(s) == MOI.EqualTo{T}
        @test s.value == one(T)
        @test length(f.terms) == array_dim
        @test f.constant == zero(T)

        for i in 1:array_dim
            t = f.terms[i]
            @test t.coefficient === one(T)
            @test t.variable == bridge.vars[i]
        end
    end

    @testset "Constraints: greater than" begin
        @test length(bridge.cons_gt) == array_dim
        for i in 1:array_dim
            @test MOI.is_valid(model, bridge.cons_gt[i])
            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_gt[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_gt[i])

            @test typeof(s) == MOI.GreaterThan{T}
            @test s.lower == zero(T)
            @test length(f.terms) == 3
            @test f.constant == -one(T)

            t1 = f.terms[1]
            @test t1.coefficient === one(T)
            @test t1.variable == x[1]

            t2 = f.terms[2]
            @test t2.coefficient === -one(T)
            @test t2.variable == x[i + 1]

            t3 = f.terms[3]
            @test t3.coefficient === one(T)
            @test t3.variable == bridge.vars[i]
        end
    end

    @testset "Constraints: less than" begin
        @test length(bridge.cons_lt) == array_dim
        for i in 1:array_dim
            @test MOI.is_valid(model, bridge.cons_lt[i])
            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_lt[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_lt[i])

            @test typeof(s) == MOI.LessThan{T}
            @test s.upper == zero(T)
            @test length(f.terms) == 2
            @test f.constant == zero(T)

            t1 = f.terms[1]
            @test t1.coefficient === one(T)
            @test t1.variable == x[1]

            t2 = f.terms[2]
            @test t2.coefficient === -one(T)
            @test t2.variable == x[i + 1]
        end
    end
end


text:
# Exports
export satellite_plan_baseline
export satellite_plan_resource_baseline

"""
Binary search method to find the first opportunity with an index greater than
the target simulation time (Epoch or elapsed time)
"""
function find_next_opp_index(problem::SatPlanningProblem, target::Union{Epoch, Real})
    # Set search indices
    l = 1
    r = length(problem.opportunities)
    m = 0

    # Perform binary search to get first opportunity with start time
    # greater or equal to the target value
    while l < r
        m = floor(Int, (l+r)/2)
        # println("l: $l, r: $r, avg: $((l+r)/2), m: $m, t: $(problem.opportunities[m].t_start)")

        if problem.opportunities[m].t_start < target
            l = m + 1
        else
            r = m
        end
    end

    return m
end

function baseline_step(problem::SatPlanningProblem, state::SatMDPState, action::Opportunity)
    # If done return same state, just done
    if typeof(action) == Done
        return SatMDPState(time=state.time,
                last_action=Done(t_start=state.time),
                last_cdo_action=state.last_cdo_action,
                request_ids=copy(state.request_ids),
                power=state.power,
                data=state.data)
    end

    # Current state time
    time0 = state.time
    time  = state.time

    # Store last viable action
    last_cdo_action = state.last_cdo_action

    # Resource
    power = state.power
    data  = state.data

    # Copy observed requests
    request_ids = copy(state.request_ids)

    # Finish state
    done = false

    # Advance time to next action
    time = action.t_start

    # Update data generate
    power_generated = 0.0
    data_generated  = problem.spacecraft[1].datagen_backorbit*(time - time0)

    if typeof(action) == Noop
        # Do nothing if Noop
    elseif typeof(action) == Sunpoint
        # Charge for duration of sunpoint
        power_generated += problem.spacecraft[1].powergen_sunpoint*(time - time0)

    elseif typeof(action) == Collect
        # Update last action
        last_cdo_action = action

        collect_generation = action.duration * problem.spacecraft[1].datagen_image

        if data + collect_generation < 1.0
            # Only perform collect if we have capacilty
            data_generated  += collect_generation
            power_generated += action.duration * problem.spacecraft[1].powergen_image

            push!(request_ids, action.location.id)
        end

    elseif typeof(action) == Contact
        # Update last action
        last_cdo_action = action

        data_generated  += action.duration * problem.spacecraft[1].datagen_downlink
        power_generated += action.duration * problem.spacecraft[1].powergen_downlink
    else
        throw(ErrorException("Action $(string(action)) is not of a known type."))
    end

    # Update power and data numbers
    power += power_generated
    data  += data_generated

    # Ensure resource limits are enforced
    if power > 1.0 power = 1.0 end

    if power < 0.0
        power = 0.0
        done = true
    end

    if data > 1.0 data = 1.0 end
    if data < 0.0 data = 0.0 end

    return SatMDPState(time=time,
            last_action=action,
            last_cdo_action=last_cdo_action,
            request_ids=request_ids,
            power=power,
            data=data)
end

function next_feasible_action(problem::SatPlanningProblem, state::SatMDPState)

    # Find index of next possible action
    idx_start = find_next_opp_index(problem::SatPlanningProblem, state.time)

    if idx_start > length(problem.opportunities)
        return nothing
    end

    # Ensure index is for next action after current state
    while problem.opportunities[idx_start].t_start <= state.time
        idx_start += 1

        if idx_start > length(problem.opportunities)
            return nothing
        end
    end

    # Return all opportunities inside planning horizon
    fopps = problem.opportunities[idx_start:end]

    # List of all possible actions
    actions = Opportunity[]

    # Populate list of possible actions
    for opp in fopps
        # Skip Collects until contact if high data
        if typeof(opp) == Collect && state.data >= 0.5
            continue
        end

        # Valuate transition to see if valid
        valid = true

        # Check to see if each transition constraint is valid
        for constraint in problem.constraints
            if valid == false
                break
            end

            valid = valid && constraint(state.last_cdo_action, opp)
        end

        # If valid transition add to edges
        if valid == true
            return opp
        end
    end

    return nothing
end

function satellite_plan_baseline(problem::SatPlanningProblem)

    # Get Initial State
    init_opp = problem.opportunities[1]
    state = SatMDPState(time=init_opp.t_start, last_action=init_opp, last_cdo_action=init_opp)

    states = SatMDPState[state]
    plan = Opportunity[problem.opportunities[1]]
    requests = Request[]

    while next_feasible_action(problem, state) != nothing
        action = next_feasible_action(problem, state)

        # Advance state
        state = baseline_step(problem, state, action)

        push!(states, state)
        push!(plan, action)
    end

    reward = 0.0

    for action in plan
        if typeof(action) == Collect
            reward += action.reward
        end
    end

    return plan, reward
end

function satellite_plan_resource_baseline(problem::SatPlanningProblem; allow_repeats::Bool=false)

    # Get Initial State
    init_opp = problem.opportunities[1]
    state = SatMDPState(time=init_opp.t_start, last_action=init_opp, last_cdo_action=init_opp)
    action = state.last_cdo_action

    states = SatMDPState[state]
    plan = Opportunity[problem.opportunities[1]]
    requests = Request[]

    while next_feasible_action(problem, state) != nothing
        # Get next feasible action
        action = next_feasible_action(problem, state)

        # If able to take collect or downlink do it
        power_generated = 0.0
        data_generated  = problem.spacecraft[1].datagen_backorbit*(action.t_start - state.time)

        if typeof(action) == Collect
            # Resources consumed by action
            data_generated  += action.duration * problem.spacecraft[1].datagen_image
            power_generated += action.duration * problem.spacecraft[1].powergen_image
        elseif typeof(action) == Contact
            # Update data generation
            data_generated  += action.duration * problem.spacecraft[1].datagen_downlink
            power_generated += action.duration * problem.spacecraft[1].powergen_downlink
        end

        if 0.0 < state.power + power_generated && state.data + data_generated < 1.0
            # Take action if able (e.g. do nothing)
        else
            # If not able sunpoint
            action = Sunpoint(t_start=action.t_start)
        end

        # println("State: $state")
        # println("Action: $action")

        # Advance state
        state = baseline_step(problem, state, action)

        push!(states, state)
        push!(plan, action)
    end

    reward = 0.0

    for action in plan
        if typeof(action) == Collect
            reward += action.reward
        end
    end

    return states, plan, reward
end


text:
module ExactDiagonalization

using Reexport: @reexport

include("FED.jl")

@reexport using .FED

end # module


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
using Pkg, DrWatson

# ╔═╡ e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, StatsBase, DataFrames
	Random.seed!(1)
end

# ╔═╡ ac3ceabe-1ee4-11eb-3935-313e96aafd62
md"## Listing 2.7"

# ╔═╡ c4bef1a4-1ee9-11eb-0b58-8be033f6fbad
begin
	A = Set(['a','e','i','o','u'])
	B = Set(['x','y','z'])
	omega = 'a':'z'
	N = 10^6
end

# ╔═╡ 23c4c304-1eea-11eb-2a57-d9b2e99a3253
begin
	df = DataFrame()
	for _ in 1:5
    	mcEst1 = sum([in(sample(omega),A) || in(sample(omega),B) for _ in 1:N])/N
    	mcEst2 = sum([in(sample(omega),union(A,B)) for _ in 1:N])/N
    	append!(df, DataFrame(:mcEst1 => mcEst1, :mcEst2 => mcEst2))
	end
	df
end

# ╔═╡ 6285dea2-1ee5-11eb-2b6a-f9906c0d3351
md"## End of listing 2.7"

# ╔═╡ Cell order:
# ╟─ac3ceabe-1ee4-11eb-3935-313e96aafd62
# ╠═de7a5e8a-1ee4-11eb-2bbf-0f0eb3e489ec
# ╠═e4eeb13a-1ee4-11eb-0113-91f9a9c3b659
# ╠═c4bef1a4-1ee9-11eb-0b58-8be033f6fbad
# ╠═23c4c304-1eea-11eb-2a57-d9b2e99a3253
# ╟─6285dea2-1ee5-11eb-2b6a-f9906c0d3351


text:
@testset "Boolean Fairness Metrics" begin
    ft = @load_toyfairtensor
    dp = DemographicParity()
    @test dp(ft) == false
    @test dp.C == 3
    A = zeros(3, 8)
    A[1, [1, 3, 5, 7]] = [4, 4, -3, -3]
    A[2, [1, 3, 5, 7]] = [3, 3, -3, -3]
    @test dp.A == A
end


text:
module Rust

export compile_code, compile_function

import Base.Libdl: dlext, dlopen, dlsym

include(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
include(joinpath(dirname(@__FILE__), "AST.jl"))

import .AST

global file_counter = -1
global generated_file_dir = ""

function compile_file_num()
    global file_counter
    file_counter += 1
    return file_counter
end

function generation_dir()
    global generated_file_dir

    if !isdir(generated_file_dir)
        generated_file_dir = mktempdir()
    end

    return generated_file_dir
end

rust_file_name(crate_name) = "$crate_name.rs"
rust_lib_name(crate_name) = "lib$crate_name.$dlext"

unique_crate_name(name="") = "$name$(compile_file_num())"

const RUST_FILE_BOILERPLATE = """
    extern crate libc;
    use libc::*;
    use std::ops::*;
"""

const RUST_FUNCTION_BOILERPLATE = "#[no_mangle]\npub extern """

function compile_code(rust_code)
    crate_name = unique_crate_name()
    dir = generation_dir()
    rust_path = joinpath(dir, rust_file_name(crate_name))

    open(rust_path, "w") do fp
        println(fp, rust_code)
    end

    lib_path = compile_file(crate_name, dir, rust_path)
    return dlopen(lib_path)
end

function compile_function(func::Function, args::ANY)
    func_name, code = AST.generate_rust(func, args)

    io = IOBuffer()
    println(io, RUST_FILE_BOILERPLATE)
    println(io)
    println(io, RUST_FUNCTION_BOILERPLATE, code)

    lib = compile_code(String(io))
    func = dlsym(lib, func_name)
end

function compile_file(crate_name, dir, rust_path)
    run(`rustc $rust_path --crate-name $crate_name --crate-type dylib -C opt-level=3 --out-dir $dir --emit=dep-info,link --extern libc=$LIBC_RLIB`)

    lib_path = joinpath(dir, rust_lib_name(crate_name))
    if !isfile(lib_path)
        error("Compilation succeeded but Rust.jl failed to generate the library in the correct place")
    end

    return lib_path
end

end # module


text:
include(joinpath("docs", "make.jl"))

cd(@__DIR__)

# fix all links to prepend them with `FranklinTemplates.jl/`
for (root, _, files) ∈ walkdir(joinpath(@__DIR__, "docs", "build"))
    for file ∈ files
        endswith(file, ".html") || continue
        path  = joinpath(root, file)
        htmls = read(path, String)
        htmls = replace(htmls, "href=\"/" => "href=\"/FranklinTemplates.jl/")
        htmls = replace(htmls, "src=\"/" => "src=\"/FranklinTemplates.jl/")
        write(path, htmls)
    end
end

using NodeJS

const JS_GHP = """
    var ghpages = require('gh-pages');
    ghpages.publish('docs/build', function(err) {});
    """
run(`$(NodeJS.nodejs_cmd()) -e $JS_GHP`)


text:
import OrdinaryDiffEq: ODEProblem, solve, Tsit5
import Zygote
using DiffEqSensitivity, Test

dynamics = (x, _p, _t) -> x

function loss(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0), params)
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = InterpolatingAdjoint(autojacvec=ZygoteVJP()))
    sum(Array(rollout)[:, end])
end

function loss2(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0), params)
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()))
    sum(Array(rollout)[:, end])
end

function loss3(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0), params)
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = InterpolatingAdjoint(autojacvec=TrackerVJP()))
    sum(Array(rollout)[:, end])
end

function loss4(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0))
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = InterpolatingAdjoint(autojacvec=ZygoteVJP()))
    sum(Array(rollout)[:, end])
end

function loss5(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0))
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = InterpolatingAdjoint(autojacvec=EnzymeVJP()))
    sum(Array(rollout)[:, end])
end

function loss6(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0))
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = BacksolveAdjoint(autojacvec=ZygoteVJP()))
    sum(Array(rollout)[:, end])
end

function loss7(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0))
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = QuadratureAdjoint(autojacvec=ZygoteVJP()))
    sum(Array(rollout)[:, end])
end

function loss8(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0))
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = QuadratureAdjoint(autojacvec=ReverseDiffVJP()))
    sum(Array(rollout)[:, end])
end

function loss9(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0))
    rollout = solve(problem, Tsit5(), u0 = u0, p = params, sensealg = QuadratureAdjoint(autojacvec=EnzymeVJP()))
    sum(Array(rollout)[:, end])
end

function loss10(params)
    u0 = zeros(2)
    problem = ODEProblem(dynamics, u0, (0.0, 1.0))
    rollout = solve(problem, Tsit5(), u0 = u0, p = params)
    sum(Array(rollout)[:, end])
end

@test Zygote.gradient(dynamics, 0.0, nothing, nothing) == (1.0,nothing,nothing)

@test Zygote.gradient(loss, nothing)[1] === nothing
@test_broken Zygote.gradient(loss2, nothing)
@test_broken Zygote.gradient(loss3, nothing)
@test Zygote.gradient(loss4, nothing)[1] === nothing
@test Zygote.gradient(loss5, nothing)[1] === nothing
@test Zygote.gradient(loss6, nothing)[1] === nothing
@test Zygote.gradient(loss7, nothing)[1] === nothing
@test Zygote.gradient(loss8, nothing)[1] === nothing
@test Zygote.gradient(loss9, nothing)[1] === nothing
@test Zygote.gradient(loss10, nothing)[1] === nothing

@test Zygote.gradient(loss, zeros(123))[1] == zeros(123)
@test Zygote.gradient(loss2, zeros(123))[1] == zeros(123)
@test Zygote.gradient(loss3, zeros(123))[1] == zeros(123)
@test Zygote.gradient(loss4, zeros(123))[1] == zeros(123)
@test Zygote.gradient(loss5, zeros(123))[1] == zeros(123)
@test Zygote.gradient(loss6, zeros(123))[1] == zeros(123)
@test_broken Zygote.gradient(loss7, zeros(123))[1] == zeros(123)
@test Zygote.gradient(loss8, zeros(123))[1] == zeros(123)
@test Zygote.gradient(loss9, zeros(123))[1] == zeros(123)
@test Zygote.gradient(loss10, zeros(123))[1] == zeros(123)


text:
using Test
using CLIMA

let
  for N = 1:10
    for FloatType in (Float64, Float32)
      (ξ, ω) = CLIMA.Mesh.Elements.lglpoints(FloatType, N)
      I∫ = CLIMA.Mesh.Grids.indefinite_integral_interpolation_matrix(ξ, ω)
      for n = 1:N
        if N == 1
          @test sum(abs.(I∫ * ξ)) < 10*eps(FloatType)
        else
          @test I∫ * ξ.^n ≈ (ξ.^(n+1) .- (-1).^(n+1)) / (n+1)
        end
      end
    end
  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
#####################################################
### From Julia 0.6 to Julia 1.0
### --- SIMPLEXN - parallel version ---
### Authors: Fabio Fatelli, Emanuele Loprevite
#####################################################

"""
	plarExtrude1(model::Tuple{Array{Array{T,1},1},Array{Array{Int64,1},1}}, pattern::Array{Int64,1}) where T<:Real
		::Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}

Algorithm for multidimensional extrusion of a simplicial complex.
Can be applied to 0-, 1-, 2-, ... simplicial models, to get a 1-, 2-, 3-, ... model.
The pattern `Array` is used to specify how to decompose the added dimension.

A `model` is a LAR model, i.e. a pair (vertices,cells) to be extruded, whereas pattern is an array of `Int64`, to be used as lateral measures of the *extruded* model. `pattern` elements are assumed as either *solid* or *empty* measures, according to their (+/-) sign.

# Example
```julia
julia> V = [[0,0], [1,0], [2,0], [0,1], [1,1], [2,1], [0,2], [1,2], [2,2]];

julia> FV = [[1,2,4],[2,3,5],[3,5,6],[4,5,7],[5,7,8],[6,8,9]];

julia> pattern = repeat([1,2,-3],outer=4);

julia> model = (V,FV);

julia> W,FW = plarExtrude1(model, pattern);
```
"""
# Generation of the output model vertices in a multiple extrusion of a LAR model
using SharedArrays

@everywhere function plarExtrude1(model::Tuple{Array{Array{T,1},1},Array{Array{Int64,1},1}}, pattern::Array{Int64,1}) where T<:Real
	V, FV = model
	d, m = length(FV[1]), length(pattern)
	dd1 = d*(d+1)
	coords = cumsum(append!([0],abs.(pattern))) # built-in function cumsum
	outVertices = @spawn [vcat(v,z) for z in coords for v in V]
	offset, outcells, rangelimit = length(V), SharedArray{Int64}(m,dd1*length(FV)), d*m
	@sync @distributed for j in 1:length(FV)
		tube = [v+k*offset for k in 0:m for v in FV[j]]
		celltube = Int64[]
		celltube = @sync @distributed (append!) for k in 1:rangelimit
			tube[k:k+d]
		end
		outcells[:,(j-1)*dd1+1:(j-1)*dd1+dd1] = reshape(celltube,dd1,m)'
	end
	p = convert(SharedArray,findall(x->x>0,pattern))
	cellGroups = SharedArray{Int64}(length(p),size(outcells)[2])
	@sync @distributed for k in 1:length(p)
		cellGroups[k,:] = outcells[p[k],:]
	end
	cellGroupsL = vec(cellGroups')
	outCellGroups = Array{Int64,1}[]
	outCellGroups = @distributed (append!) for k in 1:d+1:length(cellGroupsL)
			[cellGroupsL[k:k+d]]
		end
	return fetch(outVertices), outCellGroups
end

"""
	plarSimplexGrid1(shape::Array{Int64,1})::Tuple{Array{Array{Int64,1},1},Array{Array{Int64,1},1}}

Generate a simplicial complex decomposition of a cubical grid of ``d``-cuboids, where ``d`` is the length of `shape` array. Vertices (0-cells) of the grid have `Int64` coordinates.

# Example
```julia
julia> plarSimplexGrid1([0]) # 0-dimensional complex
# output
(Array{Int64,1}[[0]], Array{Int64,1}[])

julia> V,EV = plarSimplexGrid1([1]) # 1-dimensional complex
# output
(Array{Int64,1}[[0], [1]], Array{Int64,1}[[0, 1]])

julia> V,FV = plarSimplexGrid1([1,1]) # 2-dimensional complex
# output
(Array{Int64,1}[[0, 0], [1, 0], [0, 1], [1, 1]], Array{Int64,1}[[0, 1, 2], [1, 2, 3]])

julia> V,CV = plarSimplexGrid1([10,10,1]) # 3-dimensional complex
# output
(Array{Int64,1}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [6, 0, 0], [7, 0, 0], [8, 0, 0], [9, 0, 0]  …  [1, 10, 1], [2, 10, 1], [3, 10, 1], [4, 10, 1], [5, 10, 1], [6, 10, 1], [7, 10, 1], [8, 10, 1], [9, 10, 1], [10, 10, 1]], Array{Int64,1}[[0, 1, 11, 121], [1, 11, 121, 122], [11, 121, 122, 132], [1, 11, 12, 122], [11, 12, 122, 132], [12, 122, 132, 133], [1, 2, 12, 122], [2, 12, 122, 123], [12, 122, 123, 133], [2, 12, 13, 123]  …  [118, 228, 229, 239], [108, 118, 119, 229], [118, 119, 229, 239], [119, 229, 239, 240], [108, 109, 119, 229], [109, 119, 229, 230], [119, 229, 230, 240], [109, 119, 120, 230], [119, 120, 230, 240], [120, 230, 240, 241]])

julia> V,HV = plarSimplexGrid1([1,1,1,1]) # 4-dim cellular complex from the 4D simplex
# output
(Array{Int64,1}[[0, 0, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0], [1, 1, 0, 0], [0, 0, 1, 0], [1, 0, 1, 0], [0, 1, 1, 0], [1, 1, 1, 0], [0, 0, 0, 1], [1, 0, 0, 1], [0, 1, 0, 1], [1, 1, 0, 1], [0, 0, 1, 1], [1, 0, 1, 1], [0, 1, 1, 1], [1, 1, 1, 1]], Array{Int64,1}[[0, 1, 2, 4, 8], [1, 2, 4, 8, 9], [2, 4, 8, 9, 10], [4, 8, 9, 10, 12], [1, 2, 4, 5, 9], [2, 4, 5, 9, 10], [4, 5, 9, 10, 12], [5, 9, 10, 12, 13], [2, 4, 5, 6, 10], [4, 5, 6, 10, 12]  …  [3, 5, 9, 10, 11], [5, 9, 10, 11, 13], [2, 3, 5, 6, 10], [3, 5, 6, 10, 11], [5, 6, 10, 11, 13], [6, 10, 11, 13, 14], [3, 5, 6, 7, 11], [5, 6, 7, 11, 13], [6, 7, 11, 13, 14], [7, 11, 13, 14, 15]])
```
"""
# Generation of simplicial grids of any dimension and shape
@everywhere function plarSimplexGrid1(shape::Array{Int64,1})
	model = [Int64[]],[[0]] # the empty simplicial model
	for item in shape # no parallel
		model = plarExtrude1(model,repeat([1],item))
	end
	return model
end

"""
	plarSimplexFacets(simplices::Array{Array{Int64,1},1})::Array{Array{Int64,1},1}

Compute the `(d-1)`-skeleton (set of `facets`) of a simplicial `d`-complex.

# Example
```julia
julia> V,FV = plarSimplexGrid1([1,1])
# output
(Array{Int64,1}[[0, 0], [1, 0], [0, 1], [1, 1]], Array{Int64,1}[[0, 1, 2], [1, 2, 3]])

julia> W,CW = plarExtrude1((V,FV), [1])
# output
(Array{Int64,1}[[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0], [0, 0, 1], [1, 0, 1], [0, 1, 1], [1, 1, 1]], Array{Int64,1}[[0, 1, 2, 4], [1, 2, 4, 5], [2, 4, 5, 6], [1, 2, 3, 5], [2, 3, 5, 6], [3, 5, 6, 7]])

julia> FW = plarSimplexFacets(CW)
18-element Array{Array{Int64,1},1}:
 [0, 1, 2]
 [0, 1, 4]
 [0, 2, 4]
 [1, 2, 3]
 [1, 2, 4]
 [1, 2, 5]
 [1, 3, 5]
 [1, 4, 5]
 [2, 3, 5]
 [2, 3, 6]
 [2, 4, 5]
 [2, 4, 6]
 [2, 5, 6]
 [3, 5, 6]
 [3, 5, 7]
 [3, 6, 7]
 [4, 5, 6]
 [5, 6, 7]
```
"""
# Extraction of non-oriented (d-1)-facets of d-dimensional simplices
@everywhere using Combinatorics # for combinations() function

@everywhere function plarSimplexFacets(simplices::Array{Array{Int64,1},1})
	out = Array{Int64,1}[]
	d = length(simplices[1])
	out = @distributed (append!) for simplex in simplices
			collect(combinations(simplex,d-1))
		end
	return sort!(unique(out),lt=isless)
end

"""
	pquads2tria(model::Tuple{Array{Array{T,1},1},Array{Array{Int64,1},1}}) where T<:Real
		::Tuple{Array{Array{Float64,1},1},Array{Array{Int64,1},1}}

Give the conversion of a LAR boundary representation (B-Rep), i.e. a LAR model `V, FV` made of 2D faces, usually quads but also general polygons, into a LAR model `verts, triangles` made by triangles.

# Example
```julia
julia> model = plarSimplexGrid1([2,2,2])
# output
(Array{Int64,1}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [0, 1, 0], [1, 1, 0], [2, 1, 0], [0, 2, 0], [1, 2, 0], [2, 2, 0], [0, 0, 1]  …  [2, 2, 1], [0, 0, 2], [1, 0, 2], [2, 0, 2], [0, 1, 2], [1, 1, 2], [2, 1, 2], [0, 2, 2], [1, 2, 2], [2, 2, 2]], Array{Int64,1}[[0, 1, 3, 9], [1, 3, 9, 10], [3, 9, 10, 12], [1, 3, 4, 10], [3, 4, 10, 12], [4, 10, 12, 13], [1, 2, 4, 10], [2, 4, 10, 11], [4, 10, 11, 13], [2, 4, 5, 11]  …  [15, 21, 22, 24], [13, 15, 16, 22], [15, 16, 22, 24], [16, 22, 24, 25], [13, 14, 16, 22], [14, 16, 22, 23], [16, 22, 23, 25], [14, 16, 17, 23], [16, 17, 23, 25], [17, 23, 25, 26]])

julia> pquads2tria(model)
# output
(Array{Float64,1}[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0], [2.0, 1.0, 0.0], [0.0, 2.0, 0.0], [1.0, 2.0, 0.0], [2.0, 2.0, 0.0], [0.0, 0.0, 1.0]  …  [0.25, 1.5, 1.75], [0.75, 1.5, 1.25], [0.5, 1.75, 1.5], [0.75, 1.75, 1.75], [1.25, 1.25, 1.25], [1.5, 1.25, 1.5], [1.25, 1.5, 1.75], [1.75, 1.5, 1.25], [1.5, 1.75, 1.5], [1.75, 1.75, 1.75]], Array{Int64,1}[[27, 3, 9], [27, 9, 0], [27, 0, 1], [27, 1, 3], [28, 9, 10], [28, 10, 1], [28, 1, 3], [28, 3, 9], [29, 10, 12], [29, 12, 3]  …  [72, 16, 17], [72, 17, 23], [73, 23, 25], [73, 25, 16], [73, 16, 17], [73, 17, 23], [74, 26, 25], [74, 25, 17], [74, 17, 23], [74, 23, 26]])
```
"""
# Transformation to triangles by sorting circularly the vertices of faces
using LinearAlgebra

@everywhere function pquads2tria(model::Tuple{Array{Array{T,1},1},Array{Array{Int64,1},1}}) where T<:Real
	V, FV = model
	if typeof(V) != Array{Array{Float64,1},1}
		V = convert(Array{Array{Float64,1},1},V)
	end
	out = Array{Int64,1}[]
	nverts = length(V)-1
	for face in FV # no parallel
		arr = [V[v+1] for v in face]
		centroid = sum(arr)/length(arr)
		append!(V,[centroid])
		nverts += 1
		v1, v2 = V[face[1]+1]-centroid, V[face[2]+1]-centroid
		v3 = cross(v1,v2)
		if norm(v3) < 1/(10^3)
			v1, v2 = V[face[1]+1]-centroid, V[face[3]+1]-centroid
			v3 = cross(v1,v2)
		end
		transf = inv(copy(hcat(v1,v2,v3)'))
		verts = [(V[v+1]'*transf)'[1:end-1] for v in face]
		tcentroid = sum(verts)/length(verts)
		tverts = pmap(x->x-tcentroid,verts)
		iterator = collect(zip(tverts,face))
		rverts = [[atan(reverse(iterator[i][1])...),iterator[i][2]] for i in 1:length(iterator)]
		rvertsS = sort(rverts,lt=(x,y)->isless(x[1],y[1]))
		ord = [pair[2] for pair in rvertsS]
		append!(ord,ord[1])
		edges = [[i[2],ord[i[1]+1]] for i in enumerate(ord[1:end-1])]
		triangles = pmap(x->prepend!(x,nverts),edges)
		append!(out,triangles)
	end
	return V, out
end

max_stars_repo_name: EmaLoprevite/Simplexn

max_issues_repo_name: EmaLoprevite/Simplexn

max_forks_repo_name: EmaLoprevite/Simplexn


text:
using SatisfiabilityInterface
using Test

@testset "SatisfiabilityInterface.jl" begin

    include("graph_colouring.jl")
    include("sudoku.jl")
    # include("MOI_wrapper.jl")
    # include("JuMP.jl")

end

max_stars_repo_name: dpsanders/SatisfiabilityInterface

max_issues_repo_name: dpsanders/SatisfiabilityInterface

max_forks_repo_name: dpsanders/SatisfiabilityInterface


text:
# What if we turn all fields into independent arrays?

struct TripPaymentColumnarData
    vendor_id::Vector{Int}
    tpep_pickup_datetime::Vector{String}
    tpep_dropoff_datetime::Vector{String}
    passenger_count::Vector{Int}
    trip_distance::Vector{Float64}
    fare_amount::Vector{Float64}
    extra::Vector{Float64}
    mta_tax::Vector{Float64}
    tip_amount::Vector{Float64}
    tolls_amount::Vector{Float64}
    improvement_surcharge::Vector{Float64}
    total_amount::Vector{Float64}
end

# Reformat data
columar_records = TripPaymentColumnarData(
    [r.vendor_id for r in records],
    [r.tpep_pickup_datetime for r in records],
    [r.tpep_dropoff_datetime for r in records],
    [r.passenger_count for r in records],
    [r.trip_distance for r in records],
    [r.fare_amount for r in records],
    [r.extra for r in records],
    [r.mta_tax for r in records],
    [r.tip_amount for r in records],
    [r.tolls_amount for r in records],
    [r.improvement_surcharge for r in records],
    [r.total_amount for r in records]
);

@btime mean(columar_records.fare_amount);
#=
julia> @btime mean(columar_records.fare_amount);
  27.170 μs (1 allocation: 16 bytes)
=#


text:
using MediaCLI
using Test

@testset "MediaCLI.jl" begin
    include("discovery.jl")
end


text:
module PandorasBoxDustSpriteColorController

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/dustSpriteColorController" DustSpriteColorController(x::Integer, y::Integer, eyeTexture::String="danger/dustcreature/eyes", eyeColor::String="Red", borderColor::String="Red,Green,Blue")

const placements = Ahorn.PlacementDict(
    "Dust Sprite Color Controller (Pandora's Box)" => Ahorn.EntityPlacement(
        DustSpriteColorController
    )
)

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

Ahorn.editingOptions(entity::DustSpriteColorController) = Dict{String, Any}(
    "eyeColor" => colors
)

function Ahorn.selection(entity::DustSpriteColorController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 10, y - 10, 20, 20)
end

tintColor = (137, 251, 255, 255) ./ 255

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DustSpriteColorController, room::Maple.Room)
    Ahorn.Cairo.save(ctx)

    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1);

    Ahorn.drawCircle(ctx, 0, 0, 10, (1.0, 1.0, 1.0, 1.0))

    Ahorn.Cairo.restore(ctx)

    Ahorn.drawSprite(ctx, "danger/dustcreature/base00", 0, 0, tint=tintColor)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
module Verilog

  #support for -v suffixed verilog-style ranges.
  #support for python-style range() indexing
  include("verilog_ranges.jl")

  include("Wire.jl")
  include("wire-display.jl")

  include("booleanops.jl")
  include("arithmeticops.jl")

  #the magical verilog macro
  include("verilog_macro.jl")

  #verigen - where verilog generation happens
  include("Verigen.jl")
  #verilog generation passes these objects around.
  include("WireObject.jl")

  #codegen contains resources to create .v files and verilate them.
  if is_linux()
    include("codegen.jl")
    include("verilator-adapter.c.jl")
  end

end # module


text:
# using Revise
# using Pkg; Pkg.activate(".")
using HyperModularity
using Test

## Reading in a dataset
dataset = "contact-high-school-classes"
maxsize = 5
minsize = 2
return_labels = true
H, L = read_hypergraph_data(dataset,maxsize,minsize,return_labels)

## Run graph Louvain on reduced graph

# weighted = true means that when clique expandin hyperedge e
# weight all edges by 1/(|e|-1)
# weighted = false means they are weighted by 1, but the
# resulting graph is still integer weighted
weighted = true

# binary = true means that after clique expansion, all edges
# are unweighted. If true, "weighted" flag doesn't matter
binary = false
A = CliqueExpansion(H,weighted,binary)

# Run modularity on clique expansion
gamma = 1.0
Z = CliqueExpansionModularity(H,gamma;weighted=true,randflag=false,binary=false,clusterpenalty=0.0,maxits=1000)

# Compute the MLE for the resolution parameter based on Z
(ωᵢ, ωₒ) = computeDyadicResolutionParameter(H, Z; mode = 0)
γ = (ωᵢ - ωₒ)/(log(ωᵢ) - log(ωₒ))

# log-likelihood
loglike = dyadicLogLikelihood(H, Z, ωᵢ, ωₒ)

# Above steps can be iteratively repeated to reach global max for log-likelihood
# Good place to test log likelihood computations

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg