text:
module IssueReporter

using MetadataTools, URIParser, GitHub, DocStringExtensions

"""
$(SIGNATURES)

Takes the name of a METADATA Julia package and returns the associated git URL.

#Examples
```jldoctest
julia> IssueReporter.package_uri("Trie")
"git://github.com/JuliaArchive/Trie.jl.git"
````
"""
function package_uri(package_name::String)
  MetadataTools.get_pkg(package_name).url
end

"""
$(SIGNATURES)

Checks if the required GitHub authentication token is defined.
"""
function token_is_defined()
  if ! haskey(ENV, "GITHUB_ACCESS_TOKEN")
    secrets_path = joinpath(Pkg.dir("IssueReporter"), "secrets.jl")
    isfile(secrets_path) && include(secrets_path)
    haskey(ENV, "GITHUB_ACCESS_TOKEN") || return false
  end

  true
end

"""
$(SIGNATURES)

Returns the configured GitHub authentication token, if defined -- or throws an error otherwise.
"""
function token()
  token_is_defined() && return ENV["GITHUB_ACCESS_TOKEN"]
  error("""ENV["GITHUB_ACCESS_TOKEN"] is not set -- please make sure it's passed as a command line argument or defined in the `secrets.jl` file.""")
end

"""
$(SIGNATURES)

Accepts a GitHub authentication token and stores it in the `secrets.jl` configuration file.
The `secrets.jl` file is added to `.gitignore` so it won't be accidentally commited.
"""
function configure(token::String)
  open(joinpath(Pkg.dir("IssueReporter"), "secrets.jl"), "w") do io
    write(io, """ENV["GITHUB_ACCESS_TOKEN"] = "$token" """)
  end
end

"""
$(SIGNATURES)

Performs GitHub authentication and returns the OAuth2 object, required by further GitHub API calls.
"""
function github_auth()
  token() |> GitHub.authenticate
end

"""
$(SIGNATURES)

Converts a registered Julia package name to the corresponding GitHub "username/repo_name" string.

#Examples
```jldoctest
julia> IssueReporter.repo_id("Trie")
"JuliaArchive/Trie.jl"
````
"""
function repo_id(package_name::String)
  pkg_url = package_uri(package_name) |> URIParser.parse_url
  repo_info = endswith(pkg_url.path, ".git") ? replace(pkg_url.path, r".git$", "") : pkg_url.path

  repo_info[2:end]
end

"""
$(SIGNATURES)

Creates a new GitHub issue with the title `title` and the content `body` onto the repo corresponding to the registered package called `package_name`.
"""
function report(package_name::String, title::String, body::String)
  # GitHub.create_issue(repo_id(package_name), auth = github_auth(), params = Dict(:title => title, body => body)) # <-- commented out so we don't spam real repos
  GitHub.create_issue("essenciary/julia-by-example-test-repo", auth = github_auth(), params = Dict(:title => title, body => body))
end

end # module


text:
include("main.jl");
include("SHEMS_optimizer.jl");
include("SHEMS_optimizer_seco.jl");
include("SHEMS_optimizer_sesu.jl");

function yearly_SHEMS(h_start=1, h_end=8760, objective=1, case=1, costfactor=1.0, outputflag=1, bc_violations=79)

    # Initialize technical setup according to case______________________
    # set_SHEMS_parameters(h_start, h_end, h_predict, h_control, rolling_flag, costfactor)
    sh, hp, fh, hw, b, m = set_SHEMS_parameters(h_start, h_end, (h_end-h_start)+1, (h_end-h_start)+1,
        false, case, costfactor, outputflag);

    if objective==1   # minimize costs
        results  = SHEMS_optimizer(sh, hp, fh, hw, b, m);
    elseif objective==2  # maximize self-consumption
        results  = SHEMS_optimizer_seco(sh, hp, fh, hw, b, m, bc_violations);
    elseif objective==3  # maximize self-sufficiency
        results  = SHEMS_optimizer_sesu(sh, hp, fh, hw, b, m, bc_violations);
    end

    # write to results folder___________________________________________________
    write_to_results_file(hcat(results, ones(size(results,1))*m.h_predict), m, objective, case, costfactor)
    return nothing
end

function roll_SHEMS(h_start, h_end, h_predict, h_control, case=1, costfactor=1.0, outputflag=0)

    # Initialize technical setup__________________________________________________
    # set_SHEMS_parameters(h_start, h_end, h_predict, h_control, rolling_flag, costfactor, outputflag)
    sh, hp, fh, hw, b, m = set_SHEMS_parameters(h_start, h_end, h_predict, h_control, true, costfactor,
        case, outputflag);

    # intitial run_______________________________________________________________
    sh.soc_b, sh.soc_fh, sh.soc_hw, results_new  = SHEMS_optimizer(sh, hp, fh, hw, b, m);
    results = results_new;
    sh.h_start = m.h_start + m.h_control;

    # loop runs for rest of horizon______________________________________________
    for h = sh.h_start:m.h_control:(m.h_end-m.h_predict)
        sh.h_start = h;
        sh.soc_b, sh.soc_fh, sh.soc_hw, results_new  = SHEMS_optimizer(sh, hp, fh, hw, b, m);
        results = vcat(results, results_new);
    end

    # write to results folder____________________________________________________
    write_to_results_file(hcat(results, ones(size(results,1))*m.h_predict), m, 1, case, costfactor)
    return nothing
end

function set_SHEMS_parameters(h_start, h_end, h_predict, h_control, rolling_flag, case=1, costfactor=1.0, outputflag=0)

    # Initialize technical setup________________________________________________
    # Model_SHEMS(h_start, h_end, h_predict, h_control, big, rolling_flag, solver, mip_gap, output_flag, presolve_flag)
    m = Model_SHEMS(h_start, h_end,  h_predict, h_control, 60, rolling_flag, "Cbc", 0.005f0, outputflag, -1);
    # HeatPump(eta, rate_max)
    hp = HeatPump(1.0f0, 3.0f0);
    # ThermalStorage(eta, volume, loss, t_supply, soc_min, soc_max)
    fh = ThermalStorage(1.0f0, 10.0f0, 0.045f0, 30.0f0, 20.0f0, 22.0f0);
    hw = ThermalStorage(1.0f0, 180.0f0, 0.035f0, 45.0f0, 20.0f0, 180.0f0);
    if case==1 # base case
        # Battery(eta, soc_min, soc_max, rate_max, loss)
        b = Battery(0.95f0, 0.0f0, 13.5f0, 3.3f0, 0.00003f0);
        # SHEMS(costfactor, p_buy, p_sell, soc_b, soc_fh, soc_hw, h_start)
        sh = SHEMS(costfactor, 0.3f0, 0.1f0, 13.5f0, 22.0f0, 180.0f0, h_start);
    elseif case==2 # no battery
        # set soc_max, rate_max to zero for no battery
        b = Battery(0.95f0, 0.0f0, 0.0f0, 0.0f0, 0.00003f0);
        # soc_b zero for no battery
        sh = SHEMS(costfactor, 0.3f0, 0.1f0, 0.0f0, 22.0f0, 180.0f0, h_start);
    elseif case==3 # no grid feed-in compensation
        # Battery(eta, soc_min, soc_max, rate_max, loss)
        b = Battery(0.95f0, 0.0f0, 13.5f0, 3.3f0, 0.00003f0);
        # set p_sell to zero for no feedin tariff
        sh = SHEMS(costfactor, 0.3f0, 0.0f0, 13.5f0, 22.0f0, 180.0f0, h_start);
    elseif case==3 # no battery and no grid feed-in compensation
        # Battery(eta, soc_min, soc_max, rate_max, loss)
        b = Battery(0.95f0, 0.0f0, 0.0f0, 0.0f0, 0.00003f0);
        # set p_sell to zero for no feedin tariff
        sh = SHEMS(costfactor, 0.3f0, 0.0f0, 0.0f0, 22.0f0, 180.0f0, h_start);
    end

    return sh, hp, fh, hw, b, m
end

function write_to_results_file(results, m, objective=1, case=1, costfactor=1.0)
    date=200531;
    CSV.write("results/$(date)_results_$(m.h_predict)_$(m.h_control)_$(m.h_start)-$(m.h_end)_
        $(objective)_$(case)_$(costfactor).csv", DataFrame(results), header=["Temp_FH", "Vol_HW",
            "Soc_B", "V_HW_plus", "V_HW_minus", "T_FH_plus", "T_FH_minus", "profits", "COP_FH",
            "COP_HW","PV_DE", "B_DE", "GR_DE", "PV_B", "PV_GR", "PV_HP","GR_HP", "B_HP", "HP_FH", "HP_HW",
            "month", "day", "hour", "horizon"]);
    return nothing
end

function COPcalc(ts, t_outside)
    # Calculate coefficients of performance for every time period (1:h_predict)
    return cop = max.((5.8*ones(size(t_outside,1))) -(1. /14) * abs.((ts.t_supply*ones(size(t_outside,1)))
            -t_outside), 0);
end

ext: jl


text:
facts("Coordinates") do

    context("Create") do
        context("Brain Vision") do
            bv = BrainVision(0, 0, 0)
            context("Show") do
                show(bv)
            end
        end
        context("Talairach") do
            tal = Talairach(68.3, -26.9, 8.3)
            context("Show") do
                show(tal)
            end
        end
        context("SPM") do
            spm = SPM(68.3, -26.9, 8.3)
            context("Show") do
                show(spm)
            end
        end
        context("Unknown") do
            uk = UnknownCoordinate(1.2, 2, 3.2)
            context("Show") do
                show(uk)
            end
        end
    end

    context("Convert") do
        context("MNI -> Talairach") do

            # Values from table IV in Lancaster et al 2007
            # TODO Find better points to translate

            mni = SPM(73.7, -26.0, 7.0)
            tal = Talairach(68.3, -26.9, 8.3)
            @fact euclidean(convert(Talairach, mni), tal) --> roughly(0; atol=2)

            # As an electrode

            # Test as an electrode
            e = Electrode("test", SPM(73.7, -26.0, 7.0), Dict())
            e = conv_spm_mni2tal(e)

            @fact e.coordinate.x --> roughly(tal.x, 1.5)
            @fact e.coordinate.y --> roughly(tal.y, 1.5)
            @fact e.coordinate.z --> roughly(tal.z, 1.5)
            @fact isa(e, EEG.Sensor) --> true
            @fact isa(e, EEG.Electrode) --> true
            @fact isa(e.coordinate, EEG.Talairach) --> true

        end

        context("BrainVision -> Talairach") do

            # TODO this just tests it runs, need to check values

            bv = BrainVision(0, 0, 0)
            tal = Talairach(128, 128, 128)
            @fact euclidean(convert(Talairach, bv), tal) --> roughly(0; atol=1.5)

        end
    end

    context("Distances") do

        @fact euclidean(Talairach(0, 0, 0), Talairach(1, 1, 1)) --> sqrt(3)

        v = [0, 0, 0]
        @fact euclidean(Talairach(1, 1, 1), v) --> sqrt(3)
        @fact euclidean(v, Talairach(1, 1, 1)) --> sqrt(3)

    end
end


text:
using Distributions, StatsBase, ProgressMeter

export draw_sample_mode!, draw_sample_mode, draw_sample_gamma!, draw_sample_gamma
export draw_sample!, draw_sample

function double_iex_multinomial_edit_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::Union{SimIexInsertDelete,SimIexSplitMerge},
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )

    N = length(S_curr)
    dist = posterior.dist
    V = posterior.V
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_in_lb, K_in_ub = (K_inner.l, K_inner.u)
    data = posterior.data
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    aux_mcmc = mcmc.aux_mcmc

    δ = rand(1:mcmc.ν_ed)  # Number of edits to enact
    rem_edits = δ # Remaining edits to allocate
    j = 0 # Keeps track how many interaction have been edited
    log_ratio = 0.0

    for i = 1:N

        # If at end we just assign all remaining edits to final interaction
        if i == N
            δ_tmp = rem_edits
        # Otherwise we sample the number of edits via rescaled Binomial
        else
            p = 1/(N-i+1)
            δ_tmp = rand(Binomial(rem_edits, p)) # Number of edits to perform on ith interaction
        end

        # println("   Index $i getting $δ_tmp edits...")
        # If we sampled zero edits we skip to next iteration
        if δ_tmp == 0
            continue
        else
            j += 1 # Increment j
            # Make edits ....
            @inbounds n = length(S_curr[i])
            d = rand(0:min(n-K_in_lb, δ_tmp))
            m = n + δ_tmp - 2*d

            # Catch invalid proposals
            if (m > K_in_ub)
                # Here we just reject the proposal
                for i in 1:N
                    @inbounds copy!(S_prop[i], S_curr[i])
                end
                return 0, suff_stat_curr
            end

            ind_del = view(mcmc.ind_del, 1:d)
            ind_add = view(mcmc.ind_add, 1:(δ_tmp-d))
            vals_del = view(mcmc.vals, 1:d)

            # Sample indexing info and new entries (all in-place)
            StatsBase.seqsample_a!(1:n, ind_del)
            StatsBase.seqsample_a!(1:m, ind_add)

            # *** HERE IS DIFFERENT FROM MODEL SAMPLER ***
            # The delete_insert_informed() function does the sampling + editing
            @inbounds log_ratio += delete_insert_informed!(
                S_prop[i],
                ind_del, ind_add, vals_del,
                P)

            @inbounds mcmc.ind_update[j] = i # Store which interaction was updated

            # Add to log_ratio
            # log_prod_term += log(b - a + 1) - log(ub(m, δ_tmp) - lb(m, δ_tmp, model) +1)
            log_ratio += log(min(n-K_in_lb, δ_tmp)+1) - log(min(m-K_in_lb, δ_tmp)+1)

        end

        # Update rem_edits
        rem_edits -= δ_tmp

        # If no more left terminate
        if rem_edits == 0
            break
        end

    end

    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
    )

    if aux_init_at_prev
        tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = γ_curr * (
        sum(x -> dist(x, S_prop), aux_data)
        - sum(x -> dist(x, S_curr), aux_data)
    )

    suff_stat_prop = sum(x -> dist(x, S_prop), data)
    log_lik_ratio = γ_curr * (
        suff_stat_curr - suff_stat_prop
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    log_multinom_term = log_multinomial_ratio(S_curr, S_prop)

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio + log_multinom_term
    # println("Edit alloc proposing...")
    # println("$(S_prop)")

    # @show log_lik_ratio, aux_log_lik_ratio, log_multinom_term, log_ratio, log_α
    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        for i in view(mcmc.ind_update, 1:j)
            @inbounds copy!(S_curr[i], S_prop[i])
        end
        return 1, suff_stat_prop
    else
        for i in view(mcmc.ind_update, 1:j)
            @inbounds copy!(S_prop[i], S_curr[i])
        end
        return 0, suff_stat_curr
    end
end

function double_iex_flip_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::Union{SimIexInsertDelete,SimIexInsertDeleteProportional,SimIexSplitMerge},
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )

    dist = posterior.dist
    V = posterior.V
    data = posterior.data
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)

    aux_mcmc = mcmc.aux_mcmc
    lengths = length.(S_curr)

    δ = rand(1:min(mcmc.ν_ed,sum(lengths)))
    # @show δ, lengths
    alloc = rand_restr_bins(lengths, δ)
    ind = mcmc.ind_add
    log_ratio = 0.0

    for (key,val) in pairs(alloc)
        ind_flip = view(ind, 1:val)
        @inbounds StatsBase.seqsample_a!(1:lengths[key], ind_flip)
        log_ratio += flip_informed_excl!(
            S_prop[key],
            ind_flip,
            P
        )
    end

    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
    )

    if aux_init_at_prev
        tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    log_multinom_term = log_multinomial_ratio(S_curr, S_prop)
    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio + log_multinom_term

    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        for i in keys(alloc)
            @inbounds copy!(S_curr[i], S_prop[i])
        end
        return 1, suff_stat_prop
    else
        for i in keys(alloc)
            @inbounds copy!(S_prop[i], S_curr[i])
        end
        return 0, suff_stat_curr
    end

end

function double_iex_trans_dim_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::Union{SimIexInsertDelete,SimIexInsertDeleteProportional},
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )

    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_out_lb, K_out_ub = (K_outer.l, K_outer.u)
    K_in_ub = K_inner.u
    data = posterior.data
    dist = posterior.dist
    V = posterior.V
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    ν_td = mcmc.ν_td
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    log_ratio = 0.0

    # Enact insertion / deletion
    N = length(S_curr)
    ε = rand(1:ν_td)
    d = rand(0:min(ν_td, N))
    a = ε - d
    # Catch invalid proposal (outside dimension bounds)
    M = N - d + a  # Number of paths in proposal
    if (M < K_out_lb) | (M > K_out_ub)
        return 0, suff_stat_curr
    end
    if d > 0
        log_ratio += imcmc_multi_delete_prop_sample!(
            S_curr, S_prop,
            mcmc,
            d,
            V
        )
    end
    if a > 0
        log_ratio += imcmc_multi_insert_prop_sample!(
            S_curr, S_prop,
            mcmc,
            a,
            V, K_in_ub
        )
    end

    log_ratio += log(min(ν_td, N) + 1) - log(min(ν_td, M) + 1)

    # Now do accept-reject step (**THIS IS WHERE WE DIFFER FROM MODEL SAMPLER***)
    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
    )

    if aux_init_at_prev
        tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    log_multinom_term = log_multinomial_ratio(S_curr, S_prop)

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio + log_multinom_term

    # Note that we copy interactions between S_prop (resp. S_curr) and prop_pointers (resp .curr_pointers) by hand.
    if log(rand()) < log_α
        # Do deletions to S_curr
        ind = view(mcmc.ind_td_del, 1:d)
        for i in Iterators.reverse(ind) # If we didnt do reverse would have to update indices
            @inbounds tmp = popat!(S_curr, i)
            pushfirst!(curr_pointers, tmp)
        end
        # Do insertions
        ind = view(mcmc.ind_td_add, 1:a)
        for i in ind
            tmp = popfirst!(curr_pointers)
            copy!(tmp, S_prop[i])
            insert!(S_curr, i, tmp)
        end
        return 1, suff_stat_prop
    else
        # Remove insertions and return to prop_pointers
        ind = view(mcmc.ind_td_add, 1:a)
        for i in Iterators.reverse(ind)
            @inbounds tmp = popat!(S_prop, i)
            pushfirst!(prop_pointers, tmp)
        end
        # Re-insert deletions
        ind = view(mcmc.ind_td_del, 1:d)
        for i in ind
            @inbounds tmp = popfirst!(prop_pointers)
            copy!(tmp, S_curr[i])
            insert!(S_prop, i, tmp)
        end
        return 0, suff_stat_curr
    end

end

function double_iex_trans_dim_informed_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::T,
    p_ins::Categorical,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    ) where {T<:SimPosteriorSampler}

    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_out_lb, K_out_ub = (K_outer.l, K_outer.u)
    K_in_ub = K_inner.u
    data = posterior.data
    dist = posterior.dist
    V = posterior.V
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    ν_td = mcmc.ν_td
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    log_ratio = 0.0

    # Enact insertion / deletion
    N = length(S_curr)
    ε = rand(1:ν_td)
    d = rand(0:min(ν_td, N))
    a = ε - d
    # Catch invalid proposal (outside dimension bounds)
    M = N - d + a  # Number of paths in proposal
    if (M < K_out_lb) | (M > K_out_ub)
        return 0, suff_stat_curr
    end
    if d > 0
        log_ratio += imcmc_multi_delete_prop_sample_informed!(
            S_curr, S_prop,
            mcmc,
            d,
            p_ins
        )
    end
    if a > 0
        log_ratio += imcmc_multi_insert_prop_sample_informed!(
            S_curr, S_prop,
            mcmc,
            a,
            p_ins,
            K_in_ub
        )
    end

    log_ratio += log(min(ν_td, N) + 1) - log(min(ν_td, M) + 1)
    # Now do accept-reject step (**THIS IS WHERE WE DIFFER FROM MODEL SAMPLER***)
    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
    )

    if aux_init_at_prev
        tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_suff_stat_curr = mapreduce(x -> dist(x, S_curr), + , aux_data)
    aux_suff_stat_prop = mapreduce(x -> dist(x, S_prop), +, aux_data)
    aux_log_lik_ratio = -γ_curr * (
        aux_suff_stat_curr - aux_suff_stat_prop
    )

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    log_multinom_term = log_multinomial_ratio(S_curr, S_prop)

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio + log_multinom_term

    # println("Acceptance prob: $(exp(log_α))")
    # @show is_insert, suff_stat_curr, suff_stat_prop, log_lik_ratio, aux_suff_stat_curr, aux_suff_stat_prop, aux_log_lik_ratio, log_ratio, log_multinom_term, log_prior_ratio, log_α, exp(log_α)
    # Note that we copy interactions between S_prop (resp. S_curr) and prop_pointers (resp .curr_pointers) by hand.
    if log(rand()) < log_α
        # Do deletions to S_curr
        ind = view(mcmc.ind_td_del, 1:d)
        for i in Iterators.reverse(ind) # If we didnt do reverse would have to update indices
            @inbounds tmp = popat!(S_curr, i)
            pushfirst!(curr_pointers, tmp)
        end
        # Do insertions
        ind = view(mcmc.ind_td_add, 1:a)
        for i in ind
            tmp = popfirst!(curr_pointers)
            copy!(tmp, S_prop[i])
            insert!(S_curr, i, tmp)
        end
        return 1, suff_stat_prop
    else
        # Remove insertions and return to prop_pointers
        ind = view(mcmc.ind_td_add, 1:a)
        for i in Iterators.reverse(ind)
            @inbounds tmp = popat!(S_prop, i)
            pushfirst!(prop_pointers, tmp)
        end
        # Re-insert deletions
        ind = view(mcmc.ind_td_del, 1:d)
        for i in ind
            @inbounds tmp = popfirst!(prop_pointers)
            copy!(tmp, S_curr[i])
            insert!(S_prop, i, tmp)
        end
        return 0, suff_stat_curr
    end

end

function draw_sample_mode!(
    sample_out::Union{InteractionSequenceSample{Int}, SubArray},
    mcmc::Union{SimIexInsertDelete,SimIexInsertDeleteProportional},
    posterior::SimPosterior,
    γ_fixed::Float64;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    if loading_bar
        iter = Progress(
            length(sample_out) * lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for γ = $(γ_fixed) and n = $(posterior.sample_size) (mode conditional)....")  # Loading bar. Minimum update interval: 1 second
    end

    # Define aliases for pointers to the storage of current vals and proposals
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    β = mcmc.β
    aux_mcmc = mcmc.aux_mcmc

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for i in 1:length(S_init)
        migrate!(S_curr, curr_pointers, i, 1)
        migrate!(S_prop, prop_pointers, i, 1)
        copy!(S_curr[i], S_init[i])
        copy!(S_prop[i], S_init[i])
    end

    γ_curr = γ_fixed

    sample_count = 1 # Keeps which sample to be stored we are working to get
    i = 0 # Keeps track all samples (included lags and burn_ins)

    tr_dim_count = 0
    tr_dim_acc_count = 0
    ed_count = 0
    ed_acc_count = 0
    flp_count = 0
    flp_acc_count = 0

    aux_data = [[Int[]] for i in 1:posterior.sample_size]
    # Initialise the aux_data
    aux_model = SIM(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer)
    draw_sample!(aux_data, aux_mcmc, aux_model, burn_in=10000)

    # Evaluate sufficient statistic
    suff_stat_curr = mapreduce(
        x -> posterior.dist(S_curr, x),
        +,
        posterior.data
    )
    suff_stats = Float64[suff_stat_curr] # Storage for all sufficient stats (for diagnostics)
    P, vmap, vmap_inv = get_informed_proposal_matrix(posterior, mcmc.α)
    p_ins = get_informed_insertion_dist(posterior, mcmc.α)
    while sample_count ≤ length(sample_out)
        i += 1
        # Store value
        if (i > burn_in) & (((i - burn_in - 1) % lag)==0)
            @inbounds sample_out[sample_count] = deepcopy(S_curr)
            push!(suff_stats, suff_stat_curr)
            sample_count += 1
        end
        # W.P. do update move (accept-reject done internally by function call)
        if rand() < β
            if rand() < 0.5
                was_acc, suff_stat_curr = double_iex_multinomial_edit_accept_reject!(
                    S_curr, S_prop,
                    posterior, γ_curr,
                    mcmc, P,
                    aux_data,
                    suff_stat_curr,
                    aux_init_at_prev
                )
                ed_acc_count += was_acc
                ed_count += 1
            else
                was_acc, suff_stat_curr = double_iex_flip_accept_reject!(
                    S_curr, S_prop,
                    posterior, γ_curr,
                    mcmc, P,
                    aux_data,
                    suff_stat_curr,
                    aux_init_at_prev
                )
                flp_acc_count += was_acc
                flp_count += 1
            end
        # Else do trans-dim move. We will do accept-reject move here
        else
            was_acc, suff_stat_curr = double_iex_trans_dim_informed_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc,
                p_ins,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            tr_dim_acc_count += was_acc
            tr_dim_count += 1
        end
        if loading_bar
            next!(iter)
        end
    end
    for i in 1:length(S_curr)
        migrate!(curr_pointers, S_curr, 1, 1)
        migrate!(prop_pointers, S_prop, 1, 1)
    end
    return (
        ed_count, ed_acc_count,
        flp_count, flp_acc_count,
        tr_dim_count, tr_dim_acc_count,
        suff_stats
    )
end

function draw_sample_mode(
    mcmc::T,
    posterior::SimPosterior,
    γ_fixed::Float64;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    ) where {T<:SimPosteriorSampler}

    sample_out = Vector{InteractionSequence{Int}}(undef, desired_samples)
    draw_sample_mode!(
        sample_out,
        mcmc, posterior,
        γ_fixed,
        burn_in=burn_in, lag=lag, S_init=S_init,
        loading_bar=loading_bar,
        aux_init_at_prev=aux_init_at_prev
        )
    return sample_out

end

function (mcmc::T where {T<:SimPosteriorSampler})(
    posterior::SimPosterior,
    γ_fixed::Float64;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )
    sample_out = Vector{InteractionSequence{Int}}(undef, desired_samples)

    (
        edit_count, edit_acc_count,
        flip_count, flip_acc_count,
        trans_dim_count, trans_dim_acc_count,
        suff_stats
        ) = draw_sample_mode!(
            sample_out,
            mcmc,
            posterior, γ_fixed,
            burn_in=burn_in,
            lag=lag,
            S_init=S_init,
            loading_bar=loading_bar,
            aux_init_at_prev=aux_init_at_prev
            )

    p_measures = Dict(
        "Proportion Update Moves" => (edit_count+flip_count)/(edit_count+flip_count+trans_dim_count),
        "Edit Alloc Move Acceptance Probability" => edit_acc_count / edit_count,
        "Flip Alloc Acceptance Probability" => flip_acc_count / flip_count,
        "Trans-Dimensional Move Acceptance Probability" => trans_dim_acc_count / trans_dim_count
        )
    output = SimPosteriorModeConditionalMcmcOutput(
            γ_fixed,
            sample_out,
            posterior,
            suff_stats,
            p_measures
            )

    return output

end

# Dispersion Conditional
# ----------------------

function draw_sample_gamma!(
    sample_out::Union{Vector{Float64}, SubArray},
    mcmc::T,
    posterior::SimPosterior,
    S_fixed::InteractionSequence{Int};
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    γ_init::Float64=4.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    ) where {T<:SimPosteriorSampler}

    if loading_bar
        iter = Progress(
            length(sample_out) * lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for n = $(posterior.sample_size) (dispersion conditional)....")  # Loading bar. Minimum update interval: 1 second
    end

    # Define aliases for pointers to the storage of current vals and proposals
    ε = mcmc.ε
    aux_mcmc = mcmc.aux_mcmc

    acc_count = 0
    i = 1 # Which iteration we are on
    sample_count = 1  # Which sample we are working to get

    S_curr = deepcopy(S_fixed)
    γ_curr = γ_init
    aux_data = [[Int[]] for i in 1:posterior.sample_size]

    # Evaluate sufficient statistic
    suff_stat = mapreduce(
        x -> posterior.dist(S_curr, x),
        +,
        posterior.data
        )

    # Initialise the aux_data
    aux_model = SIM(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer
    )

    draw_sample!(aux_data, aux_mcmc, aux_model, burn_in=10000)

    while sample_count ≤ length(sample_out)
        # Store value
        if (i > burn_in) & (((i-1) % lag)==0)
            @inbounds sample_out[sample_count] = γ_curr
            sample_count += 1
        end

        γ_prop = rand_reflect(γ_curr, ε, 0.0, Inf)

        aux_model = SIM(
            S_curr, γ_prop,
            posterior.dist,
            posterior.V,
            posterior.K_inner, posterior.K_outer
        )
        if aux_init_at_prev
            tmp = deepcopy(aux_data[end])
            draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
        else
            draw_sample!(aux_data, aux_mcmc, aux_model)
        end

        # Accept reject

        log_lik_ratio = (γ_curr - γ_prop) * suff_stat
        aux_log_lik_ratio = (γ_prop - γ_curr) * sum_of_dists(aux_data, S_curr, posterior.dist)

        log_α = (
            logpdf(posterior.γ_prior, γ_prop)
            - logpdf(posterior.γ_prior, γ_curr)
            + log_lik_ratio + aux_log_lik_ratio
        )
        if log(rand()) < log_α
            γ_curr = γ_prop
            acc_count += 1
        end
        if loading_bar
            next!(iter)
        end
        i += 1

    end
    return acc_count
end

function draw_sample_gamma(
    mcmc::T,
    posterior::SimPosterior,
    S_fixed::InteractionSequence{Int};
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    γ_init::Float64,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    ) where {T<:SimPosteriorSampler}

    sample_out = Vector{Float64}(undef, desired_samples)
    draw_sample_gamme!(
        sample_out,
        mcmc, posterior,
        S_fixed,
        burn_in=burn_in, lag=lag, γ_init=γ_init,
        loading_bar=loading_bar,
        aux_init_at_prev=aux_init_at_prev
        )
    return sample_out

end

function (mcmc::T where {T<:SimPosteriorSampler})(
    posterior::SimPosterior,
    S_fixed::InteractionSequence{Int};
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    γ_init::Float64=5.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    sample_out = Vector{Float64}(undef, desired_samples)

    acc_count = draw_sample_gamma!(
            sample_out,
            mcmc,
            posterior, S_fixed,
            burn_in=burn_in,
            lag=lag,
            γ_init=γ_init,
            loading_bar=loading_bar,
            aux_init_at_prev=aux_init_at_prev
            )

    p_measures = Dict(
            "Acceptance Probability" => acc_count/desired_samples
        )

    output = SimPosteriorDispersionConditionalMcmcOutput(
            S_fixed,
            sample_out,
            posterior.γ_prior,
            posterior.data,
            p_measures
            )

    return output

end

# Joint Distribution
# ------------------

# Mode accept-reject

function accept_reject_mode!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::Union{SimIexInsertDelete,SimIexInsertDeleteProportional},
    P::CumCondProbMatrix,
    p_ins::Categorical,
    aux_data::InteractionSequenceSample{Int},
    acc_count::Vector{Int},
    count::Vector{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool=false
    )

    β = mcmc.β
    if rand() < β
        if rand() < 0.5
            was_accepted, suff_stat_curr = double_iex_multinomial_edit_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc, P,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            acc_count[1] += was_accepted
            count[1] += 1
        else
            was_accepted, suff_stat_curr = double_iex_flip_accept_reject!(
                S_curr, S_prop,
                posterior, γ_curr,
                mcmc, P,
                aux_data,
                suff_stat_curr,
                aux_init_at_prev
            )
            acc_count[2] += was_accepted
            count[2] += 1
        end
    else
        was_accepted, suff_stat_curr = double_iex_trans_dim_informed_accept_reject!(
            S_curr, S_prop,
            posterior, γ_curr,
            mcmc, p_ins,
            aux_data,
            suff_stat_curr,
            aux_init_at_prev
        )
        acc_count[3] += was_accepted
        count[3] += 1
    end
    return suff_stat_curr
end

# Gamma accept-reject

function accept_reject_gamma!(
    γ_curr::Float64,
    S_curr::InteractionSequence{Int},
    posterior::SimPosterior,
    mcmc::T,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool=false
    ) where {T<:SimPosteriorSampler}

    ε = mcmc.ε
    aux_mcmc = mcmc.aux_mcmc
    dist = posterior.dist

    γ_prop = rand_reflect(γ_curr, ε, 0.0, Inf)

    aux_model = SIM(
        S_curr, γ_prop,
        posterior.dist,
        posterior.V,
        posterior.K_inner, posterior.K_outer
    )
    if aux_init_at_prev
        @inbounds tmp = deepcopy(aux_data[end])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end
    # Accept reject

    log_lik_ratio = (γ_curr - γ_prop) * suff_stat_curr

    aux_log_lik_ratio = (γ_prop - γ_curr) * sum(dist(x,S_curr) for x in aux_data)

    log_α = (
        logpdf(posterior.γ_prior, γ_prop)
        - logpdf(posterior.γ_prior, γ_curr)
        + log_lik_ratio + aux_log_lik_ratio
    )
    if log(rand()) < log_α
        return γ_prop, 1
    else
        return γ_curr, 0
    end
end

function draw_sample!(
    sample_out_S::Union{InteractionSequenceSample{Int},SubArray},
    sample_out_gamma::Union{Vector{Float64},SubArray},
    mcmc::T,
    posterior::SimPosterior;
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    ) where {T<:SimPosteriorSampler}

    if loading_bar
        iter = Progress(
            length(sample_out_S) * lag + burn_in, # How many iters
            1,  # At which granularity to update loading bar
            "Chain for n = $(posterior.sample_size) (joint)....")  # Loading bar. Minimum update interval: 1 second
    end

    # Define aliases for pointers to the storage of current vals and proposals
    curr_pointers = mcmc.curr_pointers
    prop_pointers = mcmc.prop_pointers
    aux_mcmc = mcmc.aux_mcmc

    S_curr = InteractionSequence{Int}()
    S_prop = InteractionSequence{Int}()
    for i in 1:length(S_init)
        migrate!(S_curr, curr_pointers, i, 1)
        migrate!(S_prop, prop_pointers, i, 1)
        copy!(S_curr[i], S_init[i])
        copy!(S_prop[i], S_init[i])
    end
    γ_curr = copy(γ_init)

    sample_count = 1 # Keeps which sample to be stored we are working to get
    i = 1 # Keeps track all samples (included lags and burn_ins)

    acc_count = [0,0,0]
    count = [0,0,0]
    γ_acc_count = 0

    aux_data = [[Int[]] for i in 1:posterior.sample_size]
    # Initialise the aux_data
    aux_model = SIM(
        S_curr, γ_curr,
        posterior.dist,
        posterior.V,
        posterior.K_inner,
        posterior.K_outer)
    draw_sample!(aux_data, aux_mcmc, aux_model, burn_in=10000)
    # Initialise sufficient statistic
    suff_stat_curr = mapreduce(
        x -> posterior.dist(S_curr, x),
        +,
        posterior.data
    )
    suff_stats = Float64[suff_stat_curr]
    # Get informed proposal matrix
    P, vmap, vmap_inv = get_informed_proposal_matrix(posterior, mcmc.α)
    p_ins = get_informed_insertion_dist(posterior, mcmc.α)
    while sample_count ≤ length(sample_out_S)
        # Store values
        if (i > burn_in) & (((i-1) % lag)==0)
            @inbounds sample_out_S[sample_count] = deepcopy(S_curr)
            @inbounds sample_out_gamma[sample_count] = copy(γ_curr)
            push!(suff_stats, suff_stat_curr)
            sample_count += 1
        end

        # Update mode
        # -----------
        suff_stat_curr = accept_reject_mode!(
            S_curr, S_prop,
            posterior, γ_curr,
            mcmc, P, p_ins,
            aux_data,
            acc_count, count,
            suff_stat_curr,
            aux_init_at_prev
        )
        # Update gamma
        # ------------
        γ_curr, tmp =  accept_reject_gamma!(
            γ_curr,
            S_curr,
            posterior,
            mcmc,
            aux_data,
            suff_stat_curr,
            aux_init_at_prev
        )
        γ_acc_count += tmp
        if loading_bar
            next!(iter)
        end
        i += 1
    end

    for i in 1:length(S_curr)
        migrate!(curr_pointers, S_curr, 1, 1)
        migrate!(prop_pointers, S_prop, 1, 1)
    end

    ed_acc_prob = acc_count[1] / count[1]
    flip_acc_prob = acc_count[2] / count[2]
    td_acc_prob = acc_count[3] / count[3]
    γ_acc_prob = γ_acc_count / sum(count)
    return ed_acc_prob, flip_acc_prob, td_acc_prob, γ_acc_prob, suff_stats
end

function draw_sample(
    mcmc::Union{SimIexInsertDelete,SimIexInsertDeleteProportional},
    posterior::SimPosterior;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    sample_out_S = Vector{InteractionSequence{Int}}(undef, desired_samples)
    sample_out_gamma = Vector{Float64}(undef, desired_samples)

    draw_sample!(
        sample_out_S,
        sample_out_gamma,
        mcmc,
        posterior,
        burn_in=burn_in, lag=lag,
        S_init=S_init, γ_init=γ_init,
        loading_bar=loading_bar,
        aux_init_at_prev=aux_init_at_prev
    )

    return (S=sample_out_S, gamma=sample_out_gamma)
end

function (mcmc::T where {T<:SimPosteriorSampler})(
    posterior::SimPosterior;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    S_init::InteractionSequence{Int}=sample_frechet_mean(posterior.data, posterior.dist),
    γ_init::Float64=5.0,
    loading_bar::Bool=true,
    aux_init_at_prev::Bool=false
    )

    sample_out_S = Vector{InteractionSequence{Int}}(undef, desired_samples)
    sample_out_gamma = Vector{Float64}(undef, desired_samples)

    ed_acc_prob, flip_acc_prob, td_acc_prob, γ_acc_prob, suff_stats = draw_sample!(
        sample_out_S,
        sample_out_gamma,
        mcmc,
        posterior,
        burn_in=burn_in, lag=lag,
        S_init=S_init, γ_init=γ_init,
        loading_bar=loading_bar,
        aux_init_at_prev=aux_init_at_prev
    )

    p_measures = Dict(
        "Dipsersion acceptance probability" => γ_acc_prob,
        "Edit allocation acceptance probability" => ed_acc_prob,
        "Flip allocation acceptance probability" => flip_acc_prob,
        "Trans-dimensional acceptance probability" => td_acc_prob
        )

    return SimPosteriorMcmcOutput(
        sample_out_S,
        sample_out_gamma,
        posterior,
        suff_stats,
        p_measures
    )

end


text:
# Julia wrapper for header: /usr/include/scip/cons_sos2.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrSOS2(scip)
    ccall((:SCIPincludeConshdlrSOS2, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsSOS2(scip, cons, name, nvars, vars, weights, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsSOS2, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, weights, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicSOS2(scip, cons, name, nvars, vars, weights)
    ccall((:SCIPcreateConsBasicSOS2, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}), scip, cons, name, nvars, vars, weights)
end

function SCIPaddVarSOS2(scip, cons, var, weight)
    ccall((:SCIPaddVarSOS2, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, weight)
end

function SCIPappendVarSOS2(scip, cons, var)
    ccall((:SCIPappendVarSOS2, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}), scip, cons, var)
end

function SCIPgetNVarsSOS2(scip, cons)
    ccall((:SCIPgetNVarsSOS2, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsSOS2(scip, cons)
    ccall((:SCIPgetVarsSOS2, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetWeightsSOS2(scip, cons)
    ccall((:SCIPgetWeightsSOS2, libscip), Ptr{Cdouble}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end


text:
using VideoIO
using ImageFeatures
using Images
using ImageDraw
# Construct a AVInput object to access the video and audio streams in a video container
# io = VideoIO.open(video_file)
io = VideoIO.openvideo("test.mp4") # for testing purposes

function posindex(x)
	return floor(Int, x%600),ceil(Int, x/600)
end

function firstlast(imgdash)
	dia1pos = 0
	dia2pos = 0
	changed = true
	for (id,i) in enumerate(imgdash)

		if i == 1 && changed == true
			dia1pos = posindex(id)
			changed = false
		elseif i == 1
			dia2pos = posindex(id)
		else
			continue
		end
	end
	return dia1pos, dia2pos
end

function centerradius(dia1, dia2)
	center = (floor(Int, (dia1[1]+dia2[1])/2 ),floor(Int, (dia1[2]+dia2[2])/2))
	radius = floor(Int,sqrt((dia1[1]-dia2[1])^2 + (dia1[2]-dia2[2])^2) /2)
	return center, radius
end

for i in 1: 144
	img = load("./lotofdata/screenshot_$(i).jpeg")
    img = convert(Array{Float64}, Gray.(img)) .> 0.5
    first, last = firstlast(img)
    center, radius = centerradius(first, last)
	img = draw!(RGB.(img), Polygon(RectanglePoints(Point(center[2]-radius-5, center[1]-radius-5), Point(center[2]+radius+5, center[1]+radius+5))), RGB(0,1,0))
	img = draw(img, Polygon(RectanglePoints(Point(10, 10), Point(100, 100))), RGB{N0f8}(1))
    save("./lotofdata/results/$(i).png", img)
	break
end

# for (id, i) in enumerate(io)
#      # img = load("./lotofdata/screenshot_$(i).jpeg")
#     img = convert(Array{Float64}, Gray.(i)) .> 0.5
#     first, last = firstlast(img)
#     center, radius = centerradius(first, last)
#     img = draw(RGB.(img), Polygon(RectanglePoints(Point(center[2]-radius-5, center[1]-radius-5), Point(center[2]+radius+5, center[1]+radius+5))), RGB(0,1,0))
#     save("./images/$(id).png", img)
# end

# close(io)

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
using DashBootstrapComponents

input_group = dbc_inputgroup([
    dbc_inputgroupaddon(
        dbc_button("Random name", id="input-group-button", n_clicks=0),
        addon_type="prepend",
    ),
    dbc_input(id="input-group-button-input", placeholder="name"),
]);

callback!(
    app,
    Output("input-group-button-input", "value"),
    Input("input-group-button", "n_clicks"),
) do n_clicks
    if n_clicks > 0
        names = ["Arthur Dent", "Ford Prefect", "Trillian Astra"]
        which = n_clicks % length(names)
        return names[which + 1]
    else
        return ""
    end
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
#! /usr/bin/env julia
A = zeros(Int, 1000, 1000);

# #1 @ 249,597: 20x15
parser = r"@ (\d+),(\d+): (\d+)x(\d+)"

function getClaimArea(horiz, vert, width, height)
  horiz += 1
  vert += 1
  view(A, horiz:(horiz+width-1), vert:(vert+height-1))
end

function paintA!(horiz, vert, width, height)
  b = getClaimArea(horiz, vert, width, height)

  # faster than b .+= 1
  map!(x -> x + 1, b, b)
end

claims = map(readlines()) do l
  m = match(parser, l)
  map(s -> parse(Int, s), m.captures)
end

foreach(c -> paintA!(c...), claims)

claimIdx = findfirst(claims) do claim
  claimArea = getClaimArea(claim...)
  all(claimArea .== 1)
end

# display(A)
println("part1:", count(x -> x > 1, A))
println("part2:", claimIdx)

max_stars_repo_name: paolobueno/advent-of-julia

max_issues_repo_name: paolobueno/advent-of-julia

max_forks_repo_name: paolobueno/advent-of-julia


text:
struct ODEProblem{F,T<:Tuple{Number,Number},U<:AbstractVector,P<:AbstractVector}
    f::F
    tspan::T
    u0::U
    θ::P
end

abstract type ODESolver end

struct Euler{T} <: ODESolver
    dt::T
end

function (solver::Euler)(prob::ODEProblem, u, t)
    f, θ, dt  = prob.f, prob.θ, solver.dt
    (u + dt*f(u,θ), t+dt)
end

function solve(prob::ODEProblem, solver::ODESolver)
    t = prob.tspan[1]; u = prob.u0
    us = [u]; ts = [t]
    while t < prob.tspan[2]
        (u,t) = solver(prob, u, t)
        push!(us,u)
        push!(ts,t)
    end
    ts, reduce(hcat,us)
end

# Define & Solve ODE

function lotkavolterra(x,θ)
    α, β, γ, δ = θ
    x₁, x₂ = x

    dx₁ = α*x₁ - β*x₁*x₂
    dx₂ = δ*x₁*x₂ - γ*x₂

    [dx₁, dx₂]
end

θ = [0.1,0.2,0.3,0.2]
u0 = [1.0,1.0]
tspan = (0.,100.)
prob = ODEProblem(lotkavolterra,tspan,u0,θ)


text: println("Julia says, Hello")


text:
@testset "OneHotArray" begin
  idx = rand(5:10)
  ot = OneHot(10)
  o = OneHot{10}(idx)
  ind = rand(1:10, 5, 5)
  oa = OneHotArray{10}(ind)

  @test PrimitiveOneHot.onehotsize(oa) == 10
  @test OneHotArray(15, oa) |> PrimitiveOneHot.onehotsize == 15

  @test size(oa) == (10, 5, 5)
  @test oa[:, 5, 5] == oa.onehots[5, 5]
  @test oa[:, 5, :] == OneHotArray(oa.onehots[5, :])
  @test oa[ind[1], 1, 1]
  @test oa[5, :, 2] == collect(oa)[5, :, 2]
  @test oa[5, :, :] == collect(oa)[5, :, :]

  @test reinterpret(Int32, oa) == Int32.(ind)

  @test vcat(o, o) == vcat(collect(o), collect(o))
  @test hcat(o, o).onehots == [o,o]
  @test cat(o, o; dims=1) == vcat(collect(o), collect(o))
  @test cat(o, o; dims=2).onehots == [o,o]

  @test vcat(oa, oa) == vcat(collect(oa), collect(oa))
  @test hcat(oa, oa).onehots == vcat(oa.onehots, oa.onehots)
  @test cat(oa, oa; dims=1) == vcat(collect(oa), collect(oa))
  @test cat(oa, oa; dims=2).onehots == vcat(oa.onehots, oa.onehots)
  @test cat(oa, oa; dims=3).onehots == cat(oa.onehots, oa.onehots; dims=2)

  @test reshape(oa, 10, 25) isa OneHotArray
  @test reshape(oa, 10, :) isa OneHotArray
  @test reshape(oa, :, 25) isa OneHotArray
  @test reshape(oa, 50, :) isa Base.ReshapedArray{Bool, 2}
  @test reshape(oa, 5, 10, 5) isa Base.ReshapedArray{Bool, 3}

end


text:
# code for testing our package
using PSVrelief, Test

@testset "steam" begin
    # test case steam
    # From API 520, Part I, 2000 (section 3.7.2)
    # W = 69615 kg/h
    # P = 11032*1.1 + 101.3 = 12236 kPa
    # saturated
    # API gets 1100 mm2, and this is with a Napier factor of 51.5, we are using 51.45
    @test (abs(PSVrelief.PSVsteamSize(69615.0, 12236.0, "Sat") - 1100.0) < 2.0)
end

@testset "vapour" begin
    # test case vapour
    # From API 520, Part I, 2013 (section 5.6.3.2.3)
    # W = 24270.0 kg/h
    # P = 670.0 kPa
    # MW = 51
    # Z = 0.9
    # k = Cp/Cv = 1.11, or C = 328.
    # T = 348 - 273
    # API gets 3698 mm2
    # W, P, Tcelcius, MW, k, Z
    @test (abs(PSVrelief.PSVvaporSize(24270.0, 670.0, 75.0, 51.0, 1.11, 0.9) - 3698) < 2.0)
end

@testset "liquid1" begin
    # test case liquid1
    # from API 520, Part I, 2013 (section 5.8.1.3)
    # flow = 6814 litres/min
    # density = 900 kg/m3
    # P = 1896 kpag (incudes 10%)
    # Pback = 345 kPag
    # viscosity 2000 SSU (API uses archaic units, the table gives 431.7 cSt)
    # API gets 3180 mm2 with a backpressure correction of 0.97, or 3085 without the backpressure correciton
    # PSVliquidSize(W, P, Pback, d, mu)
    # PSVliquidSize(900.0*6814*60.0/1000.0, 1896.0, 345.0, 900.0, 431.7*0.900)
    @test (abs(PSVrelief.PSVliquidSize(900.0*6814*60.0/1000.0, 1896.0, 345.0, 900.0, 431.7*0.900) - 3085.0) < 10.0)
end

@testset "liquid2" begin
    # test case liquid2
    # from API 520, Part I, 2013 (section 5.8.1.3)
    # flow = 6814 litres/min
    # density = 900 kg/m3
    # P = 1896 kpag (incudes 10%)
    # Pback = 345 kPag
    # viscosity 2000 SSU (API uses archaic units, the table gives 431.7 cSt)
    # we will flip around the sizing case to get the rating case
    # API gets 3180 mm2 with a backpressure correction of 0.97, or 3085 without the backpressure correciton
    # PSVliquidSize(W, P, Pback, d, mu)
    targetFlow = 900.0*6814*60.0/1000.0
    targetArea = PSVrelief.PSVliquidSize(targetFlow, 1896.0, 345.0, 900.0, 431.7*0.900)
    @test (abs(PSVrelief.PSVliquidRate(targetArea, 1896.0, 345.0, 900.0, 431.7*0.900) - targetFlow) < 1.0)
end

max_stars_repo_name: kevindorma/PSVrelief

max_issues_repo_name: kevindorma/PSVrelief

max_forks_repo_name: kevindorma/PSVrelief


text:
export BrickletIndustrialDualAnalogInVoltageCallbackThreshold
struct BrickletIndustrialDualAnalogInVoltageCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletIndustrialDualAnalogInCalibration
struct BrickletIndustrialDualAnalogInCalibration
    offset::Vector{Integer}
    gain::Vector{Integer}
end

export BrickletIndustrialDualAnalogInIdentity
struct BrickletIndustrialDualAnalogInIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialDualAnalogIn
"""
Measures two DC voltages between -35V and +35V with 24bit resolution each
"""
mutable struct BrickletIndustrialDualAnalogIn <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialDualAnalogIn(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_dual_analog_in")
        deviceInternal = package.BrickletIndustrialDualAnalogIn(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage for the given channel.

If you want to get the voltage periodically, it is recommended to use the
:cb:`Voltage` callback and set the period with
:func:`Set Voltage Callback Period`.
"""
function get_voltage(device::BrickletIndustrialDualAnalogIn, channel)
    return device.deviceInternal.get_voltage(channel)
end

export set_voltage_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Voltage` callback is triggered
periodically for the given channel. A value of 0 turns the callback off.

The :cb:`Voltage` callback is only triggered if the voltage has changed since the
last triggering.
"""
function set_voltage_callback_period(device::BrickletIndustrialDualAnalogIn, channel, period)
    device.deviceInternal.set_voltage_callback_period(channel, period)
end

export get_voltage_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Voltage Callback Period`.
"""
function get_voltage_callback_period(device::BrickletIndustrialDualAnalogIn, channel)
    return device.deviceInternal.get_voltage_callback_period(channel)
end

export set_voltage_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Voltage Reached` callback for the given
channel.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the voltage is *outside* the min and max values"
 "'i'",    "Callback is triggered when the voltage is *inside* the min and max values"
 "'<'",    "Callback is triggered when the voltage is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the voltage is greater than the min value (max is ignored)"
"""
function set_voltage_callback_threshold(device::BrickletIndustrialDualAnalogIn, channel, option, min, max)
    device.deviceInternal.set_voltage_callback_threshold(channel, option, min, max)
end

export get_voltage_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Voltage Callback Threshold`.
"""
function get_voltage_callback_threshold(device::BrickletIndustrialDualAnalogIn, channel)
    return device.deviceInternal.get_voltage_callback_threshold(channel)
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callback

* :cb:`Voltage Reached`

is triggered, if the threshold

* :func:`Set Voltage Callback Threshold`

keeps being reached.
"""
function set_debounce_period(device::BrickletIndustrialDualAnalogIn, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletIndustrialDualAnalogIn)
    return device.deviceInternal.get_debounce_period()
end

export set_sample_rate
"""
    $(SIGNATURES)

Sets the sample rate. The sample rate can be between 1 sample per second
and 976 samples per second. Decreasing the sample rate will also decrease the
noise on the data.
"""
function set_sample_rate(device::BrickletIndustrialDualAnalogIn, rate)
    device.deviceInternal.set_sample_rate(rate)
end

export get_sample_rate
"""
    $(SIGNATURES)

Returns the sample rate as set by :func:`Set Sample Rate`.
"""
function get_sample_rate(device::BrickletIndustrialDualAnalogIn)
    return device.deviceInternal.get_sample_rate()
end

export set_calibration
"""
    $(SIGNATURES)

Sets offset and gain of MCP3911 internal calibration registers.

See MCP3911 datasheet 7.7 and 7.8. The Industrial Dual Analog In Bricklet
is already factory calibrated by Tinkerforge. It should not be necessary
for you to use this function
"""
function set_calibration(device::BrickletIndustrialDualAnalogIn, offset, gain)
    device.deviceInternal.set_calibration(offset, gain)
end

export get_calibration
"""
    $(SIGNATURES)

Returns the calibration as set by :func:`Set Calibration`.
"""
function get_calibration(device::BrickletIndustrialDualAnalogIn)
    return device.deviceInternal.get_calibration()
end

export get_adc_values
"""
    $(SIGNATURES)

Returns the ADC values as given by the MCP3911 IC. This function
is needed for proper calibration, see :func:`Set Calibration`.
"""
function get_adc_values(device::BrickletIndustrialDualAnalogIn)
    return device.deviceInternal.get_adc_values()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialDualAnalogIn)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialDualAnalogIn, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
@testset "Tests: types" begin
    gate = QCO.GateData("H_1", 2)
    @test gate.isreal
    gate = QCO.GateData("S_3", 3)
    @test !gate.isreal

    # Add more tests for U3 and R gates.
end


text:
using KernelAbstractions

@kernel function compute_vorticity!(vorticity, fx, fy, k)
    i, j = @index(Global, NTuple)
    vorticity[i, j, k] = fx[i, j-1] - fx[i, j] - fy[i-1, j] + fy[i, j]
end

#=
@inline δx⁻(i, j, f::AbstractArray) = @inbounds f[i, j] - f[i-1, j]
@inline δx⁺(i, j, f::AbstractArray) = @inbounds f[i+1, j] - f[i, j]

@inline δy(i, j, f::AbstractArray) = @inbounds f[i, j] - f[i, j-1]

@inline δx(i, j, f::Function, args...) = f(i, j, args...) - f(i-1, j, args...)
@inline δy(i, j, f::Function, args...) = f(i, j, args...) - f(i, j-1, args...)

@inline δ²x(i, j, f, args...) = δx(i, j, δx, f, args...)

fx(i, j, args...) = #
fy(i, j, args...) = #

@kernel function compute_vorticity!(vorticity, k, args...)
    i, j = @index(Global, NTuple)
    vorticity[i, j, k] = fx(i, j-1, args...) - fx(i, j, args...) - fy(i-1, j, args...) + fy(i, j, args...)
end
=#


text:
"""
Daylight generator table for the
relative spectral power distribution (SPD) SD(λ)
of a D series illuminant

λ : wavelength in nm
S0: mean
S1: 1st char. vect.
S2: 2nd char. vect.

of the CCT dependent daylight SPDs
"""
const daylight_generator_table=[
       300    0.04    0.02   0.00;
       305    3.02    2.26   1.00;
       310    6.00    4.50   2.00;
       315   17.80   13.45   3.00;
       320   29.60   22.40   4.00;
       325   42.45   32.20   6.25;
       330   55.30   42.00   8.50;
       335   56.30   41.30   8.15;
       340   57.30   40.60   7.80;
       345   59.55   41.10   7.25;
       350   61.80   41.60   6.70;
       355   61.65   39.80   6.00;
       360   61.50   38.00   5.30;
       365   65.15   40.20   5.70;
       370   68.80   42.40   6.10;
       375   66.10   40.45   4.55;
       380   63.40   38.50   3.00;
       385   64.60   36.75   2.10;
       390   65.80   35.00   1.20;
       395   80.30   39.20   0.05;
       400   94.80   43.40  -1.10;
       405   99.80   44.85  -0.80;
       410  104.80   46.30  -0.50;
       415  105.35   45.10  -0.60;
       420  105.90   43.90  -0.70;
       425  101.35   40.50  -0.95;
       430   96.80   37.10  -1.20;
       435  105.35   36.90  -1.90;
       440  113.90   36.70  -2.60;
       445  119.75   36.30  -2.75;
       450  125.60   35.90  -2.90;
       455  125.55   34.25  -2.85;
       460  125.50   32.60  -2.80;
       465  123.40   30.25  -2.70;
       470  121.30   27.90  -2.60;
       475  121.30   26.10  -2.60;
       480  121.30   24.30  -2.60;
       485  117.40   22.20  -2.20;
       490  113.50   20.10  -1.80;
       495  113.30   18.15  -1.65;
       500  113.10   16.20  -1.50;
       505  111.95   14.70  -1.40;
       510  110.80   13.20  -1.30;
       515  108.65   10.90  -1.25;
       520  106.50    8.60  -1.20;
       525  107.65    7.35  -1.10;
       530  108.80    6.10  -1.00;
       535  107.05    5.15  -0.75;
       540  105.30    4.20  -0.50;
       545  104.85    3.05  -0.40;
       550  104.40    1.90  -0.30;
       555  102.20    0.95  -0.15;
       560  100.00    0.00   0.00;
       565   98.00   -0.80   0.10;
       570   96.00   -1.60   0.20;
       575   95.55   -2.55   0.35;
       580   95.10   -3.50   0.50;
       585   92.10   -3.50   1.30;
       590   89.10   -3.50   2.10;
       595   89.80   -4.65   2.65;
       600   90.50   -5.80   3.20;
       605   90.40   -6.50   3.65;
       610   90.30   -7.20   4.10;
       615   89.35   -7.90   4.40;
       620   88.40   -8.60   4.70;
       625   86.20   -9.05   4.90;
       630   84.00   -9.50   5.10;
       635   84.55  -10.20   5.90;
       640   85.10  -10.90   6.70;
       645   83.50  -10.80   7.00;
       650   81.90  -10.70   7.30;
       655   82.25  -11.35   7.95;
       660   82.60  -12.00   8.60;
       665   83.75  -13.00   9.20;
       670   84.90  -14.00   9.80;
       675   83.10  -13.80  10.00;
       680   81.30  -13.60  10.20;
       685   76.60  -12.80   9.25;
       690   71.90  -12.00   8.30;
       695   73.10  -12.65   8.95;
       700   74.30  -13.30   9.60;
       705   75.35  -13.10   9.05;
       710   76.40  -12.90   8.50;
       715   69.85  -11.75   7.75;
       720   63.30  -10.60   7.00;
       725   67.50  -11.10   7.30;
       730   71.70  -11.60   7.60;
       735   74.35  -11.90   7.80;
       740   77.00  -12.20   8.00;
       745   71.10  -11.20   7.35;
       750   65.20  -10.20   6.70;
       755   56.45   -9.00   5.95;
       760   47.70   -7.80   5.20;
       765   58.15   -9.50   6.30;
       770   68.60  -11.20   7.40;
       775   66.80  -10.80   7.10;
       780   65.00  -10.40   6.80;
       785   65.50  -10.50   6.90;
       790   66.00  -10.60   7.00;
       795   63.50  -10.15   6.70;
       800   61.00   -9.70   6.40;
       805   57.15   -9.00   5.95;
       810   53.30   -8.30   5.50;
       815   56.10   -8.80   5.80;
       820   58.90   -9.30   6.10;
       825   60.40   -9.55   6.30;
       830   61.90   -9.80   6.50]

# Source: D. R. W. Hunt, “Measuring Colour”, third edition
# Appendix 5.2, CIE Standard Illuminants
#      λ        B        C
const b_c_series_table=[
       300     0.00    0.00;
       305     0.00    0.00;
       310     0.00    0.00;
       315     0.00    0.00;
       320     0.02    0.01;
       325     0.26    0.20;
       330     0.50    0.40;
       335     1.45    1.55;
       340     2.40    2.70;
       345     5.60    7.00;
       350     9.60   12.90;
       365    12.40   17.20;
       370    15.20   21.40;
       375    18.80   27.50;
       380    22.40   33.00;
       385    26.85   39.92;
       390    31.30   47.40;
       395    36.18   55.17;
       400    41.30   63.30;
       405    46.62   71.81;
       410    52.10   80.60;
       415    57.70   89.53;
       420    63.20   98.10;
       425    68.37  105.80;
       430    73.10  112.40;
       435    77.31  117.75;
       440    80.80  121.50;
       445    83.44  123.45;
       450    85.40  124.00;
       455    86.88  123.60;
       460    88.30  123.10;
       465    90.08  123.30;
       470    92.00  123.80;
       475    93.75  124.09;
       480    95.20  123.90;
       485    96.23  122.92;
       490    96.50  120.70;
       495    95.71  116.90;
       500    94.20  112.10;
       505    92.37  106.98;
       510    90.70  102.30;
       515    89.65   98.81;
       520    89.50   96.90;
       525    90.43   96.78;
       530    92.20   98.00;
       535    94.46   99.94;
       540    96.90  102.10;
       545    99.16  103.95;
       550   101.00  105.20;
       555   102.20  105.67;
       560   102.80  105.30;
       565   102.92  104.11;
       570   102.60  102.30;
       575   101.90  100.15;
       580   101.00   97.80;
       585   100.07   95.43;
       590    99.20   93.20;
       595    98.44   91.22;
       600    98.00   89.70;
       605    98.08   88.83;
       610    98.50   88.40;
       615    99.06   88.19;
       620    99.70   88.10;
       625   100.36   88.06;
       630   101.00   88.00;
       635   101.56   87.86;
       640   102.20   87.80;
       645   103.05   87.99;
       650   103.90   88.20;
       655   104.59   88.20;
       660   105.00   87.90;
       665   105.08   87.22;
       670   104.90   86.30;
       675   104.55   85.30;
       680   103.90   84.00;
       685   102.84   82.21;
       690   101.60   80.20;
       695   100.38   78.24;
       700    99.10   76.30;
       705    97.70   74.36;
       710    96.20   72.40;
       715    94.60   70.40;
       720    92.90   68.30;
       725    91.10   66.30;
       730    89.40   65.40;
       735    88.00   62.80;
       740    86.90   61.50;
       745    85.90   60.20;
       750    85.20   59.20;
       755    84.80   58.50;
       760    84.70   58.10;
       765    84.90   58.00;
       770    85.40   58.20;
       775    86.10   58.50;
       780    87.00   59.10;
       785     0.00   78.91;
       790     0.00   79.55;
       795     0.00   76.48;
       800     0.00   73.40;
       805     0.00   68.66;
       810     0.00   63.92;
       815     0.00   67.35;
       820     0.00   70.78;
       825     0.00   72.61;
       830     0.00   74.44]

# Source: D. R. W. Hunt, “Measuring Colour”, Third edition
# Appendix 5.3, “Spectral power distributions of representative fluorescent lamps”
# checked for correctness
const f_series_table=[
       380     1.87   1.18   0.82   0.57   1.87   1.05   2.56   1.21   0.90   1.11   0.91   0.96;
       385     2.36   1.48   1.02   0.70   2.35   1.31   3.18   1.50   1.12   0.80   0.63   0.64;
       390     2.94   1.84   1.26   0.87   2.92   1.63   3.84   1.81   1.36   0.62   0.46   0.45;
       395     3.47   2.15   1.44   0.98   3.45   1.90   4.53   2.13   1.60   0.57   0.37   0.33;
       400     5.17   3.44   2.57   2.01   5.10   3.11   6.15   3.17   2.59   1.48   1.29   1.19;
       405    19.49  15.69  14.36  13.75  18.91  14.80  19.37  13.08  12.80  12.16  12.68  12.48;
       410     6.13   3.85   2.70   1.95   6.00   3.43   7.37   3.83   3.05   2.12   1.59   1.12;
       415     6.24   3.74   2.45   1.59   6.11   3.30   7.05   3.45   2.56   2.70   1.79   0.94;
       420     7.01   4.19   2.73   1.76   6.85   3.68   7.71   3.86   2.86   3.74   2.46   1.08;
       425     7.79   4.62   3.00   1.93   7.58   4.07   8.41   4.42   3.30   5.14   3.38   1.37;
       430     8.56   2.06   3.28   2.10   8.31   4.45   9.15   5.09   3.82   6.75   4.49   1.78;
       435    43.67  34.98  31.85  30.28  40.76  32.61  44.14  34.10  32.62  34.39  33.94  29.05;
       440    16.94  11.81   9.47   8.03  16.06  10.74  17.52  12.42  10.77  14.86  12.13   7.90;
       445    10.72   6.27   4.02   2.55  10.32   5.48  11.35   7.68   5.84  10.40   6.95   2.65;
       450    11.35   6.63   4.25   2.70  10.91   5.78  12.00   8.60   6.57  10.76   7.19   2.71;
       455    11.89   6.93   4.44   2.82  11.40   6.03  12.58   9.46   7.25  10.67   7.12   2.65;
       460    12.37   7.19   4.59   2.91  11.83   6.25  13.08  10.24   7.86  10.11   6.72   2.49;
       465    12.75   7.40   4.72   2.99  12.17   6.41  13.45  10.84   8.35   9.27   6.13   2.33;
       470    13.00   7.54   4.80   3.04  12.40   6.52  13.71  11.33   8.75   8.29   5.46   2.10;
       475    13.15   7.62   4.86   3.08  12.54   6.58  13.88  11.71   9.06   7.29   4.79   1.91;
       480    13.23   7.65   4.87   3.09  12.58   6.59  13.95  11.98   9.31   7.91   5.66   3.01;
       485    13.17   7.62   4.85   3.09  12.52   6.56  13.93  12.17   9.48  16.64  14.29  10.83;
       490    13.13   7.62   4.88   3.14  12.47   6.56  13.82  12.28   9.61  16.73  14.96  11.88;
       495    12.85   7.45   4.77   3.06  12.20   6.42   3.64  12.32   9.68   10.44  8.97   6.88;
       500    12.52   7.28   4.67   3.00  11.89   6.28  13.43  12.35   9.74   5.94   4.72   3.43;
       505    12.20   7.15   4.62   2.98  11.61   6.20  13.25  12.44   9.88   3.34   2.33   1.49;
       510    11.83   7.05   4.62   3.01  11.33   6.19  13.08  12.55  10.04   2.35   1.47   0.92;
       515    11.50   7.04   4.73   3.14  11.10   6.30  12.93  12.68  10.26   1.88   1.10   0.71;
       520    11.22   7.16   4.99   3.41  10.96   6.60  12.78  12.77  10.48   1.59   0.89   0.60;
       525    11.05   7.47   5.48   3.90  10.97   7.12  12.60  12.72  10.63   1.47   0.83   0.63;
       530    11.03   8.04   6.25   4.69  11.16   7.94  12.44  12.60  10.78   1.80   1.18   1.10;
       535    11.18   8.88   7.34   5.81  11.54   9.07  12.33  12.43  10.96   5.71   4.90   4.56;
       540    11.53  10.01   8.78   7.32  12.12  10.49  12.26  12.22  11.18  40.98  39.59  34.40;
       545    27.74  24.88  23.82  22.59  27.78  25.22  29.52  28.96  27.71  73.69  72.84  65.40;
       550    17.05  16.64  16.14  15.11  17.73  17.46  17.05  16.51  16.29  33.61  32.61  29.48;
       555    13.55  14.59  14.59  13.88  14.47  15.63  12.44  11.79  12.28   8.24   7.52   7.16;
       560    14.33  16.16  16.63  16.33  15.20  17.22  12.58  11.76  12.74   3.38   2.83   3.08;
       565    15.01  17.56  18.49  18.68  15.77  18.53  12.72  11.77  13.21   2.47   1.96   2.47;
       570    15.52  18.62  19.95  20.64  16.10  19.43  12.83  11.84  13.65   2.14   1.67   2.27;
       575    18.29  21.47  23.11  24.28  18.54  21.97  15.46  14.61  16.57   4.86   4.43   5.09;
       580    19.55  22.79  24.69  26.26  19.50  23.01  16.75  16.11  18.14  11.45  11.28  11.96;
       585    15.48  19.29  21.41  23.28  15.39  19.41  12.83  12.34  14.55  14.79  14.76  15.32;
       590    14.91  18.66  20.85  22.94  14.64  18.56  12.67  12.53  14.65  12.16  12.73  14.27;
       595    14.15  17.73  19.93  22.14  13.72  17.42  12.45  12.72  14.66   8.97   9.74  11.86;
       600    13.22  16.54  18.67  20.91  12.69  16.09  12.19  12.92  14.61   6.52   7.33   9.28;
       605    12.19  15.21  17.22  19.43  11.57  14.64  11.89  13.12  14.50   8.31   9.72  12.31;
       610    11.12  13.80  15.65  17.74  10.45  13.15  11.60  13.34  14.39  44.12  55.27  68.53;
       615    10.03  12.36  14.04  16.00   9.35  11.68  11.35  13.61  14.40  34.55  42.58  53.02;
       620     8.95  10.95  12.45  14.42   8.29  10.25  11.12  13.87  14.47  12.09  13.18  14.67;
       625     7.96   9.65  10.95  12.56   7.32   8.95  10.95  14.07  14.62  12.15  13.16  14.38;
       630     7.02   8.40   9.51  10.93   6.41   7.74  10.76  14.20  14.72  10.52  12.26  14.71;
       635     6.20   7.32   8.27   9.52   5.63   6.69  10.42  14.16  14.55   4.43   5.11   6.46;
       640     5.42   6.31   7.11   8.18   4.90   5.71  10.11  14.13  14.40   1.95   2.07   2.57;
       645     4.73   5.43   6.09   7.01   4.26   4.87  10.04  14.34  14.58   2.19   2.34   2.75;
       650     4.15   4.68   5.22   6.00   3.72   4.16  10.02  14.50  14.88   3.19   3.58   4.18;
       655     3.64   4.02   4.45   5.11   3.25   3.55  10.11  14.46  15.51   2.77   3.01   3.44;
       660     3.20   3.45   3.80   4.36   2.83   3.02   9.87  14.00  15.47   2.29   2.48   2.81;
       665     2.81   2.96   3.23   3.69   2.49   2.57   8.65  12.58  13.20   2.00   2.14   2.42;
       670     2.47   2.55   2.75   3.13   2.19   2.20   7.27  10.99  10.57   1.52   1.54   1.64;
       675     2.18   2.19   2.33   2.64   1.93   1.87   6.44   9.98   9.18   1.35   1.33   1.36;
       680     1.93   1.89   1.99   2.24   1.71   1.60   5.83   9.22   8.25   1.47   1.46   1.49;
       685     1.72   1.64   1.70   1.91   1.52   1.37   5.41   8.62   7.57   1.79   1.94   2.14;
       690     1.67   1.53   1.55   1.70   1.48   1.29   5.04   8.07   7.03   1.74   2.00   2.34;
       695     1.43   1.27   1.27   1.39   1.26   1.05   4.57   7.39   6.35   1.02   1.20   1.42;
       700     1.29   1.10   1.09   1.18   1.13   0.91   4.12   6.71   5.72   1.14   1.35   1.61;
       705     1.19   0.99   0.96   1.03   1.05   0.81   3.77   6.16   5.25   3.32   4.10   5.04;
       710     1.08   0.88   0.83   0.88   0.96   0.71   3.46   5.63   4.80   4.49   5.58   6.98;
       715     0.96   0.76   0.71   0.74   0.85   0.61   3.08   5.03   4.29   2.05   2.51   3.19;
       720     0.88   0.68   0.62   0.64   0.78   0.54   2.73   4.46   3.80   0.49   0.57   0.71;
       725     0.81   0.61   0.54   0.54   0.72   0.48   2.47   4.02   3.43   0.24   0.27   0.30;
       730     0.77   0.56   0.49   0.49   0.68   0.44   2.25   3.66   3.12   0.21   0.23   0.26;
       735     0.75   0.54   0.46   0.46   0.67   0.43   2.06   3.36   2.86   0.21   0.21   0.23;
       740     0.73   0.51   0.43   0.42   0.65   0.40   1.90   3.09   2.64   0.24   0.24   0.28;
       745     0.68   0.47   0.39   0.37   0.61   0.37   1.75   2.85   2.43   0.24   0.24   0.28;
       750     0.69   0.47   0.39   0.37   0.62   0.38   1.62   2.65   2.26   0.21   0.20   0.21;
       755     0.64   0.43   0.35   0.33   0.59   0.35   1.54   2.51   2.14   0.17   0.24   0.17;
       760     0.68   0.46   0.38   0.35   0.62   0.39   1.45   2.37   2.02   0.21   0.32   0.21;
       765     0.69   0.47   0.39   0.36   0.64   0.41   1.32   2.15   1.83   0.22   0.26   0.19;
       770     0.61   0.40   0.33   0.31   0.55   0.33   1.17   1.89   1.61   0.17   0.16   0.15;
       775     0.52   0.33   0.28   0.26   0.47   0.26   0.99   1.61   1.38   0.12   0.12   0.10;
       780     0.43   0.27   0.21   0.19   0.40   0.21   0.81   1.32   1.12   0.09   0.09   0.05]

# Source: D. R. W. Hunt, “Measuring Colour”, Third edition
# Date from: Appendix 5.6, “Gas discharge lamps”
# relative SPDs, relative to HPS at 560 nm
# columns:
# Low pressure sodium, High pressure sodium, Mercury MB, Mercury MBF, Mercury MBTF, Mercury HMI, Xenon
# checked for correctness
const gas_discharge_table=[
       380     0.1     5.55     5.14    4.53     4.73   116.39    93.03;
       385     0.0     6.11     4.28    3.91     4.01   114.92    94.59;
       390     0.0     7.25     7.49    6.68     6.40   115.62    96.33;
       395     0.1     8.02     5.26    4.81     4.85   107.07   100.56;
       400     0.0     9.22     7.71    7.73     6.77   108.19   102.81;
       405     0.2    11.37   121.05  107.15   102.94   123.74   100.40;
       410     0.0    11.84    61.57   54.92    54.53   125.42   100.84;
       415     0.0    13.33     5.73    6.08     5.97   123.88   101.45;
       420     0.0    15.39     3.67    4.33     4.19   125.56   102.57;
       425     0.1    16.95     4.29    4.97     4.82   116.39   101.94;
       430     0.1    19.30     6.26    7.15     6.75    89.57   101.29;
       435     0.0    23.70   168.75  146.23   144.91   122.69   101.54;
       440     0.0    25.45   140.57  123.51   131.86   120.03   103.74;
       445     0.0    23.36     5.84    6.82     7.00    71.29   103.67;
       450     0.0    38.22     3.86    4.66     4.96    65.93   110.30;
       455     0.0    27.95     3.22    3.79     4.28    73.25   112.78;
       460     0.0    19.25     2.84    3.46     4.08    80.53   116.52;
       465     0.1    54.71     2.72    3.27     4.05    72.48   129.56;
       470     0.1    57.56     2.74    3.49     4.34    66.65   141.07;
       475     0.0    25.63     3.00    3.70     4.68    67.82   126.45;
       480     0.0     9.94     2.77    3.58     4.68    66.60   115.94;
       485     0.1     7.81     3.10    3.94     5.20    59.40   118.42;
       490     0.0    12.68     5.24    6.04     7.27    60.43   111.75;
       495     0.3    73.42     6.11    6.39     8.03    62.92   113.67;
       500     0.7   150.21     2.88    3.73     5.59    64.03   105.17;
       505     0.0    22.10     2.58    3.49     5.55    68.89   103.88;
       510     0.0    13.02     2.39    3.57     5.73    57.96   102.90;
       515     0.2    39.88     2.58    3.75     6.16    57.28   102.71;
       520     0.1    18.38     2.28    3.67     6.24    54.24   102.29;
       525     0.0    11.27     2.43    3.88     6.70    50.74   101.90;
       530     0.1    11.97     2.53    4.38     7.29    52.25   101.54;
       535     0.1    12.91     3.67    6.74     9.43    53.26   101.47;
       540     0.1    15.01     4.29    9.93    12.21    56.32   101.12;
       545     0.1    24.05   195.18  198.10   195.59   114.85   101.43;
       550     0.1    60.75   181.63  171.97   178.84   119.19   100.98;
       555     0.2    86.90     6.93   10.57    14.46    59.24   100.75;
       560     0.2   100.00     3.74    7.55    11.60    56.91   100.44;
       565     2.1   274.28     3.50    6.61    11.26    62.96   100.28;
       570     8.1   638.70     3.79    6.76    11.92    65.18   100.19;
       575     1.3   317.64   100.00  100.00   100.00   100.00   100.00;
       580     1.4   468.03   292.87  273.65   265.84   158.96    99.86;
       585   131.8   524.92    39.62   42.59    47.39    81.58   100.47;
       590  1000.0   157.99     4.23   15.05    19.07    82.42   100.33;
       595   150.6   521.00     3.38   27.64    28.62    72.41    99.04;
       600     2.0   663.20     2.74   13.65    18.67    74.44    97.17;
       605     1.3   505.88     2.97   10.75    16.68    62.19    96.65;
       610     1.3   372.87     2.78   24.34    26.49    61.75    96.84;
       615     3.5   372.50     3.10   83.81    69.19    61.00    98.50;
       620     1.9   272.75     2.65  149.80   119.79    60.76   100.21;
       625     0.6   188.67     3.51   58.07    53.70    62.91    99.91;
       630     0.6   160.93     2.67   17.99    23.74    55.78    97.80;
       635     0.8   141.89     2.91   10.51    18.90    50.61    96.84;
       640     1.8   127.07     2.66    8.09    18.23    52.97    98.92;
       645     0.6   115.00     2.76    7.86    18.56    49.83    99.18;
       650     0.7   104.10     2.66   10.98    21.14    46.82   101.99;
       655     0.4    98.04     2.78   10.76    21.58    47.96    98.78;
       660     0.3    91.55     2.68    7.25    19.31    56.13    97.14;
       665     0.2    85.00     2.81    6.29    18.98    53.38    98.43;
       670     0.5    92.90     2.96    6.06    19.18    67.89   100.30;
       675     0.2    90.94     3.18    5.91    19.68    73.25   101.29;
       680     0.0    73.37     2.75    5.18    19.43    47.32   101.97;
       685     0.0    56.00     2.95    5.48    20.20    50.08   109.55;
       690     0.1    45.52     6.51    9.83    24.27    44.96   110.88;
       695     0.5    41.23     5.19   18.35    30.79    56.34   102.08;
       700     0.2    38.28     2.82   46.48    51.26    69.68    94.49;
       705     0.6    37.39     3.17   35.63    43.56    36.16    93.05;
       710     0.1    35.86     4.29   15.82    29.82    28.41    98.10;
       715     0.0    35.22     3.22    5.92    23.00    27.15   106.20;
       720     0.1    33.84     2.76    4.61    22.22    28.21    98.24;
       725     0.2    33.52     2.97    4.48    22.86    34.95    95.69;
       730     0.0    32.60     2.74    4.11    22.84    34.94   101.73;
       735     0.0    32.01     2.97    4.23    23.63    27.68   108.38;
       740     0.2    32.51     2.71    3.92    23.41    24.19   102.62;
       745     0.0    32.19     2.97    4.08    24.35    22.61   100.28;
       750     0.4    33.72     2.63    3.84    24.30    21.79    97.33;
       755     0.2    34.04     3.00    4.13    26.32    30.47    97.45;
       760     0.1    31.90     2.67    3.67    25.18    33.75   101.71;
       765    23.6    54.53     3.10    3.99    26.40    36.16   137.21;
       770    55.8    81.51     2.99    3.79    26.22    23.39   105.22;
       775    11.1    40.67     3.86    4.61    27.96    26.55    84.41;
       780     0.0    32.49     2.64    3.41    27.08    24.38    80.55]


text:
module CommunalHelperConnectedDreamBlock

using ..Ahorn, Maple
using Ahorn.CommunalHelper
using Ahorn.CommunalHelperEntityPresets: CustomDreamBlockData

@mapdefdata Entity "CommunalHelper/ConnectedDreamBlock" ConnectedDreamBlock CustomDreamBlockData

const placements = Ahorn.PlacementDict(
    "Connected Dream Block (Normal) (Communal Helper)" => Ahorn.EntityPlacement(
        ConnectedDreamBlock,
        "rectangle",
    ),
    "Connected Dream Block (Feather Mode) (Communal Helper)" => Ahorn.EntityPlacement(
        ConnectedDreamBlock,
        "rectangle",
        Dict{String,Any}(
            "featherMode" => true,
        ),
    ),
    "Connected Dream Block (Normal, One Use) (Communal Helper)" => Ahorn.EntityPlacement(
        ConnectedDreamBlock,
        "rectangle",
        Dict{String,Any}(
            "oneUse" => true,
        ),
    ),
    "Connected Dream Block (Feather Mode, One Use) (Communal Helper)" => Ahorn.EntityPlacement(
        ConnectedDreamBlock,
        "rectangle",
        Dict{String,Any}(
            "featherMode" => true,
            "oneUse" => true,
        ),
    ),
)

Ahorn.minimumSize(entity::ConnectedDreamBlock) = 8, 8
Ahorn.resizable(entity::ConnectedDreamBlock) = true, true

Ahorn.selection(entity::ConnectedDreamBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ConnectedDreamBlock)
    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    renderDreamBlock(ctx, 0, 0, width, height, entity.data)
end

end


text:
"""
    sum_of_sines(domain, c, ω, ϕ)

Create initial condition function `u₀` and an antiderivative `U₀`.
The function is a sum of sines of amplitudes `c`, frequencies `ω` and phase-shifts `ϕ`.
"""
function sum_of_sines(domain, c, ω, ϕ)
    φ(ω, ϕ, x) = sin(ω * x - ϕ)
    Φ(ω, ϕ, x) = ω ≈ 0 ? -sin(ϕ) * x : -cos(ω * x - ϕ) / ω
    u₀(x) = sum(c .* φ.(ω, ϕ, x))
    U₀(x) = sum(c .* Φ.(ω, ϕ, x))
    u₀, U₀
end


text:
#=
This file contains all of the functions to perform the integration on triangles.
=#

mutable struct Quadrature1d
    wg::Array{Float64, 1}
    xg::Array{Float64, 1}
end

mutable struct Quadrature2d
    wg::Array{Float64, 1}
    xg::Array{Float64, 2}
end

Quadrature1d() = Quadrature1d(zeros(Float64, 0), zeros(Float64, 0))
Quadrature2d() = Quadrature2d(zeros(Float64, 0), zeros(Float64, 0, 0))

function quadGL1d(ng::Int64)
    #= Hard-coded quadrature rules on [-1,1] =#
    if ng == 1
        xg = [0.0]; wg = [2.0];
    elseif ng == 2
        xg = [-0.577350269189625764509148780502, 0.577350269189625764509148780502]; wg = [1.0, 1.0]
    elseif ng == 3
        xg = [-0.774596669241483377035853079956, 0.0, 0.774596669241483377035853079956]; wg = [5. /9., 8. /9., 5. /9.]
    elseif ng == 4
        xg = [-0.861136311594052575223946488893, -0.339981043584856264802665759103, 0.339981043584856264802665759103, 0.861136311594052575223946488893]; wg = [0.347854845137453857373063949222, 0.652145154862546142626936050778, 0.652145154862546142626936050778,0.347854845137453857373063949222];
    elseif ng == 5
        xg = [-0.906179845938663992797626878299, -0.538469310105683091036314420700, 0.0, 0.538469310105683091036314420700, 0.906179845938663992797626878299]; wg = [0.236926885056189087514264040720, 0.478628670499366468041291514836, 0.568888888888888888888888888889, 0.478628670499366468041291514836, 0.236926885056189087514264040720];
    else
        xg = []; wg = [];
    end
    return Quadrature1d(wg, xg)
end

function quadGL1d(ng::Int64, a::Float64, b::Float64)
    #= Quadrature rule on [a,b] =#
    ref = quadGL1d(ng)
    cw = (b - a)/2.
    cx = (a + b)/2.
    for i = 1:ng
        ref.wg[i] *= cw
        ref.xg[i] =  cx + cw*ref.xg[i]
    end
    return ref
end

function quadGL1d(ref::Quadrature1d, a::Float64, b::Float64)
    #= Same as before, but from an already-existing rule =#
    cw = (b - a)/2.
    cx = (a + b)/2.
    xg = copy(ref.xg)
    wg = copy(ref.wg)
    for i = 1:length(ref.wg)
        wg[i] *= cw
        xg[i] =  cx + cw*xg[i]
    end
    return Quadrature1d(wg, xg)
end

function quadGL2d(ng::Int64)
    #= Hard-coded 2d quadrature rules on {(0,0),(1,0),(0,1)} =#
    if ng == 1
        xg = reshape([1/3.; 1/3.], (2, 1)); wg = [0.5]
    elseif ng == 3
        xg = [2/3 1/6 1/6; 1/6 2/3 1/6]; wg = [1/6., 1/6., 1/6.];
    elseif ng == 4
        xg = [1/3 0.6 0.2 0.2; 1/3 0.2 0.6 0.2]; wg = [-27/96., 25. /96., 25. /96., 25. /96.];
    elseif ng == 6
        a1 = 0.816847572980459; b1 = 0.091576213509771; w1 = 0.5*0.109951743655322;
        a2 = 0.108103018168070; b2 = 0.445948490915965; w2 = 0.5*0.223381589678011;
        xg = [a1 b1 b1 a2 b2 b2; b1 a1 b1 b2 a2 b2]; wg = [w1, w1, w1, w2, w2, w2];
    else
        xg = []; wg = [];
    end
    return Quadrature2d(wg, xg)
end

function quadGL2d(ng::Int64, T::Array{Float64, 2})
    #= Returns quadrature on T =#
    ref = quadGL2d(ng)
    A  = T[:, 1]
    v1 = T[:, 2] - T[:, 1]
    v2 = T[:, 3] - T[:, 1]
    ac = abs(v1[1]*v2[2] - v1[2]*v2[1])
    xg = zeros(Float64, 2, ng)
    wg = zeros(Float64, ng)
    for i=1:ng
        xg[1, i] = A[1] + v1[1]*ref.xg[1, i] + v2[1]*ref.xg[2, i]
        xg[2, i] = A[2] + v1[2]*ref.xg[1, i] + v2[2]*ref.xg[2, i]
        wg[i]    = ac * ref.wg[i]
    end
    return Quadrature2d(wg, xg)
end

function quadGL2d(ref::Quadrature2d, T::Array{Float64, 2})
    #= Quadrature on triangle T from a reference quadrature ref =#
    ng = length(ref.wg)
    A  = T[:, 1]
    v1 = T[:, 2] - T[:, 1]
    v2 = T[:, 3] - T[:, 1]
    ac = abs(v1[1]*v2[2] - v1[2]*v2[1])
    xg = zeros(Float64, 2, ng)
    wg = zeros(Float64, ng)
    for i=1:ng
        xg[1, i] = A[1] + v1[1]*ref.xg[1, i] + v2[1]*ref.xg[2, i]
        xg[2, i] = A[2] + v1[2]*ref.xg[1, i] + v2[2]*ref.xg[2, i]
        wg[i]    = ac * ref.wg[i]
    end
    return Quadrature2d(wg, xg)
end

function integrate(quad::Quadrature1d, f::Function)
    #= Computes the integral of f using quad =#
    ng  = length(quad.wg)
    tot = zero(f(0.))
    for n=1:ng
        @inbounds tot += quad.wg[n] * f(quad.xg[n])
    end
    return tot
end

function integrate(quad::Quadrature2d, f::Function)
    #= same as before, but in 2d =#
    ng  = length(quad.wg)
    tot = zero(f([0., 0.]))
    for n=1:ng
        @inbounds tot += quad.wg[n] * f(quad.xg[:, n])
    end
    return tot
end