text:
### DCOPF without renewable energy sources
function DCOPF(time,
			   nodes::Nodes,
			   lines::Lines,
			   powerplants::PowerPlants,
			   )

	T = time
	G = powerplants.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	lines_pmax = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		lines_pmax[j, k] = lines.pmax[id]
		lines_pmax[k, j] = lines.pmax[id]
	end

	# Create the B_prime matrix in one step
	B_prime = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		X = lines.reactance[id]				# reactance

		B_prime[j, k] = 1/X					# Fill YBUS with susceptance
		B_prime[k, j] = 1/X					# Symmetric matrix
	end

	for j in 1:size(B_prime)[1]
		B_prime[j, j] = -sum(B_prime[j, :])
	end

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	K = nodes.id

	# Initialise JuMP model: DCOPF
	DCOPF_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(DCOPF_mod, P[T, G] >= 0)
	@variable(DCOPF_mod, Θ[T, J])
	@variable(DCOPF_mod, P_inj[T, J])
	@variable(DCOPF_mod, P_flow[T, K, J])

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	# Constraints
	# Market clearing/power balance
	@constraint(DCOPF_mod, MarketClearing[t = T],
	    sum(P[t, g] for g in G) == nodes.systemload[t]/Sbase);

	# Upper generation limit
	@constraint(DCOPF_mod, GenerationLimitUp[g = G, t = T],
	    P[t, g] <= powerplants.pmax[g]/Sbase);

	# Lower generation limit
	@constraint(DCOPF_mod, GenerationLimitDown[g = G, t = T],
	    P[t, g] >= powerplants.pmin[g]/Sbase);

	# Ramp-up limit
	@constraint(DCOPF_mod, RampUp[g = G, t = T; t > t[1]],
		(P[t, g]) - (P[t-1, g]) <= powerplants.rup[g])

	# Ramp-down limit
	@constraint(DCOPF_mod, RampDown[g = G, t = T; t > t[1]],
		-((P[t, g]) - (P[t-1, g])) <= powerplants.rdn[g])

	### DC power flow constraints
	# Power injection balance
	@constraint(DCOPF_mod, PowerInjectionBal[j = J, t = T],
		sum(P[t, g] for g in getKeyVector(powerplants.node, j)) -	# all generation units connected to a node
			nodes.load[j][t]/Sbase == P_inj[t, j]);

	# Power injection angle
	@constraint(DCOPF_mod, PowerInjectionAng[j = J, t = T],
		sum(B_prime[j, k] * (Θ[t, j] - Θ[t, k]) for k in K)
			== P_inj[t, j]);

	@constraint(DCOPF_mod, LinePowerFlow[k = K, j = J, t = T; j != k],
		B_prime[j, k] * (Θ[t, j] - Θ[t, k]) == P_flow[t, j, k]);

	@constraint(DCOPF_mod, LinePowerFlowMax[k = K, j = J, t = T],
		P_flow[t, j, k] <= lines_pmax[j, k]/Sbase);

	@constraint(DCOPF_mod, LinePowerFlowMin[k = K, j = J, t = T],
		P_flow[t, j, k] >= -lines_pmax[j, k]/Sbase);

	# Objective function
	@objective(DCOPF_mod, Min,
	    sum(powerplants.mc[g][1] * P[t, g] * Sbase for g in G, t in T));

	# Initiate optimisation process
	JuMP.optimize!(DCOPF_mod)

	# Export results
	P_opt = JuMP.value.(P)*Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	P_inj_opt = JuMP.value.(P_inj)*Sbase

	return(DCOPF_mod,
	       P_opt,
	       Θ_opt,
	       P_inj_opt,
	       P_flow_opt
	       )
end

### DCOPF with renewable energy sources
function DCOPF(time,
			   nodes::Nodes,
			   lines::Lines,
			   powerplants::PowerPlants,
			   renewables::Renewables
			   )

	T = time
	G = powerplants.unit
	R = renewables.unit

	slack = nodes.slack
	Sbase = nodes.bmva # MVA

	lines_pmax = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		lines_pmax[j, k] = lines.pmax[id]
		lines_pmax[k, j] = lines.pmax[id]
	end

	# Create the B_prime matrix in one step
	B_prime = zeros(length(nodes.id), length(nodes.id))

	for id in lines.id
		j = lines.from[id]
		k = lines.to[id]

		X = lines.reactance[id]				# reactance

		B_prime[j, k] = 1/X					# Fill YBUS with susceptance
		B_prime[k, j] = 1/X					# Symmetric matrix
	end

	for j in 1:size(B_prime)[1]
		B_prime[j, j] = -sum(B_prime[j, :])
	end

	###################
	# MODEL ###########
	###################

	# Create a subset of buses I and J
	J = nodes.id
	K = nodes.id

	# Initialise JuMP model: DCOPF
	DCOPF_mod = JuMP.Model(with_optimizer(Gurobi.Optimizer))

	# Variables
	@variable(DCOPF_mod, P[T, G] >= 0)
	@variable(DCOPF_mod, P_R[T, R] >= 0)
	@variable(DCOPF_mod, Θ[T, J])
	@variable(DCOPF_mod, P_inj[T, J])
	@variable(DCOPF_mod, P_flow[T, K, J])

	# Fix the voltage angle of the slack bus
	for t in T
		fix(Θ[t, slack], 0)
	end

	# Constraints
	# Market clearing/power balance
	@constraint(DCOPF_mod, MarketClearing[t = T],
	    sum(P[t, g] for g in G) +
			sum(P_R[t, r] for r in R) == nodes.systemload[t]/Sbase);

	# Upper generation limit
	@constraint(DCOPF_mod, GenerationLimitUp[g = G, t = T],
	    P[t, g] <= powerplants.pmax[g]/Sbase);

	# Lower generation limit
	@constraint(DCOPF_mod, GenerationLimitDown[g = G, t = T],
	    P[t, g] >= powerplants.pmin[g]/Sbase);

	# Upper generation limit for renewables
	@constraint(DCOPF_mod, ResGenerationLimitUp[r = R, t = T],
	    P_R[t, r] <= renewables.generation[r][t]/Sbase);

	# Ramp-up limit
	@constraint(DCOPF_mod, RampUp[g = G, t = T; t > t[1]],
		(P[t, g]) - (P[t-1, g]) <= powerplants.rup[g])

	# Ramp-down limit
	@constraint(DCOPF_mod, RampDown[g = G, t = T; t > t[1]],
		-((P[t, g]) - (P[t-1, g])) <= powerplants.rdn[g])

	### DC power flow constraints
	# Power injection balance
	@constraint(DCOPF_mod, PowerInjectionBal[j = J, t = T],
		sum(P[t, g] for g in getKeyVector(powerplants.node, j)) +
			sum(P_R[t, r] for r in getKeyVector(renewables.node, j)) -
				nodes.load[j][t]/Sbase == P_inj[t, j]);

	# Power injection angle
	@constraint(DCOPF_mod, PowerInjectionAng[j = J, t = T],
		sum(B_prime[j, k] * (Θ[t, j] - Θ[t, k]) for k in K)
			== P_inj[t, j]);

	@constraint(DCOPF_mod, LinePowerFlow[k = K, j = J, t = T; j != k],
		B_prime[j, k] * (Θ[t, j] - Θ[t, k]) == P_flow[t, j, k]);

	@constraint(DCOPF_mod, LinePowerFlowMax[k = K, j = J, t = T],
		P_flow[t, j, k] <= lines_pmax[j, k]/Sbase);

	@constraint(DCOPF_mod, LinePowerFlowMin[k = K, j = J, t = T],
		P_flow[t, j, k] >= -lines_pmax[j, k]/Sbase);

	# Objective function
	@objective(DCOPF_mod, Min,
	    sum(powerplants.mc[g][1] * P[t, g] * Sbase for g in G, t in T) +
			sum(renewables.mc[r] * P_R[t, r] * Sbase for r in R, t in T));

	# Initiate optimisation process
	JuMP.optimize!(DCOPF_mod)

	# Export results
	P_opt = JuMP.value.(P)*Sbase
	P_R_opt = JuMP.value.(P_R)*Sbase
	Θ_opt = (JuMP.value.(Θ)) * 360/(2π)
	P_flow_opt = JuMP.value.(P_flow)*Sbase
	P_inj_opt = JuMP.value.(P_inj)*Sbase

	return(DCOPF_mod,
	       P_opt,
		   P_R_opt,
	       Θ_opt,
	       P_inj_opt,
	       P_flow_opt
	       )
end


text:
#__precompile__(false) # see also warmup=false

module GPlot

using DocStringExtensions: SIGNATURES

using Parameters
using Colors
using DelimitedFiles: writedlm
using Statistics: quantile, mean

import Base: |>, take!, isempty

export Figure, gcf, gca, clf, cla, clo, cll, erase!,
    continuous_preview, preview,
    # General set function
    set, set_palette,
    # Layout
    layout, subplot,
    # Drawings
    plot!, scatter!, fill_between!, hist!, bar!,  # append mode
    plot, scatter, fill_between, hist, bar,       # new/overwrite mode
    # -- special drawings
    boxplot, heatmap, heatmap_ticks,
    # -- 3D drawings
    plot3!, scatter3!,
    plot3, scatter3,
    # Objects
    text, vline, hline, line, box, colorbar,
    # Axis / Axes
    # -- legend
    legend,
    # -- titles
    title, xtitle, x2title, ytitle, y2title,
    xlabel, x2label, ylabel, y2label,
    # -- lims
    xlim, x2lim, ylim, y2lim,
    # -- ticks
    xticks, x2ticks, yticks, y2ticks,
    # -- scale
    xscale, x2scale, yscale, y2scale,
    # -- misc
    axis, xaxis, x2axis, yaxis, y2axis, grid, math,
    # Preview / rendering / saving
    preview, savefig, isempty, debug_gle,
    # Simple macros for tex strings and color strings
    @t_str, @tex_str, @c_str

const PALETTE_1 = [ # imitated from tableau 10 - 2
    RGB(0.33, 0.47, 0.64),
    RGB(0.90, 0.57, 0.26),
    RGB(0.82, 0.37, 0.36),
    RGB(0.51, 0.70, 0.69),
    RGB(0.42, 0.62, 0.35),
    RGB(0.91, 0.79, 0.37),
    RGB(0.66, 0.49, 0.62),
    RGB(0.95, 0.63, 0.66),
    RGB(0.59, 0.46, 0.38),
    RGB(0.72, 0.69, 0.67) ]

const GP_ENV = Dict{String, Any}(
    "BACKEND"      => nothing,           # set in __init__
    "HAS_BACKEND"  => false,             # set in __init__ (hopefully to true)
    "VERBOSE"      => true,              # whether to write things in the REPL
    "TMP_PATH"     => mktempdir(),       # where intermedate files will be stored
    "DEL_INTERM"   => true,              # delete intermediate files
    "SHOW_GSERR"   => false,             # GLE related, ghostscript errors
    "PALETTE"      => PALETTE_1,         # default color palette
    "SIZE_PALETTE" => length(PALETTE_1),
    "CONT_PREVIEW" => true,              # whether to continuously preview or not
    "FORCE_PREVIEW" => false,            # even if cont prev is false (see render)
    )

# ain't nobody got time to write long type names
const ∅   = nothing
const AV  = AbstractVector
const AVM = AbstractVecOrMat
const AM  = AbstractMatrix
const AVR = AV{<:Real}
const T2F = NTuple{2,Float64}
const T3F = NTuple{3,Float64}
const T2R = NTuple{2,Real}
const T3R = NTuple{3,Real}

const PT_TO_CM   = 0.0352778         # 1pt in cm
const Option{T}  = Union{Nothing, T} # a useful type for optional values
const CanMiss{T} = Union{Missing, T}

include("utils.jl")

# Type of objects
include("types/style.jl")
include("types/drawing.jl")
include("types/drawing2.jl")
include("types/drawing3d.jl")
include("types/ax_elems.jl")
include("types/legend.jl")
include("types/object.jl")
include("types/ax.jl")
include("types/figure.jl")

# Naming of properties: specific and shared
include("set_prop/dicts_gle.jl")
include("set_prop/dicts_shared.jl")

# Set properties of objects
include("set_prop/style.jl")
include("set_prop/drawing.jl")
include("set_prop/object.jl")
include("set_prop/ax_elems.jl")
include("set_prop/legend.jl")
include("set_prop/ax.jl")
include("set_prop/figure.jl")
include("set_prop/properties.jl")

# Write objects to GLE buffer
include("apply_gle/style.jl")
include("apply_gle/drawing.jl")
include("apply_gle/drawing2.jl")
include("apply_gle/drawing3d.jl")
include("apply_gle/object.jl")
include("apply_gle/ax_elems.jl")
include("apply_gle/legend.jl")
include("apply_gle/ax.jl")
include("apply_gle/figure.jl")

# Main call for elements
include("drawing.jl")   # plot, scatter, fill, hist, bar
include("drawing2.jl")  # boxplot, heatmap
include("drawing3d.jl") # plot3, scatter3, surface
include("object.jl")    # text, hline, vline, box
include("ax.jl")
include("ax_elem.jl")
include("legend.jl")
include("figure.jl")
include("layout.jl")
include("gle_sub.jl")

# Rendering of things (preview/savefig)
include("render.jl")

# Extra few utils now that all types have been defined
include("utils2.jl")

include("build.jl")

# warmup script
include("init.jl")

end # module

ext: jl


text:
# ---------------------------------------------
# Royalty low level computations
# ---------------------------------------------

function η12(obs::ObservationRoyalty, theta::AbstractVector, zm::Real)
    l = _y(obs)
    # L = num_choices(obs)
    # η1 = l == 1 ? typemin(zm) : theta_royalty_κ(obs, theta, l-1) - zm
    # η2 = l == L ? typemax(zm) : theta_royalty_κ(obs, theta, l)   - zm
    η1, η2 = kappa_sums(obs, theta) .- zm
    η1 < η2 || throw(error("$η1 = η1 < η2=$η2 is false"))
    return η1, η2
end

@inline function kappa_sums(obs, theta::AbstractVector{T})::NTuple{2,T} where {T}
    l = _y(obs)
    L = num_choices(obs)

    kappas = theta_royalty_κ(obs, theta)
    k0 = first(kappas)
    ksum = 0.5*sumsq(kappas, 2, l-1)

    if l == 1
        return (typemin(T), k0)
    elseif 1 < l < L
        tmp = k0+ksum
        klsq = 0.5*kappas[l]^2
        return (tmp, tmp+klsq)
    elseif l == L
        return (k0+ksum, typemax(T))
    else
        throw(DomainError(l))
    end
end

# see
# https://github.com/cossio/TruncatedNormal.jl/blob/fc904152f2da11a257e3ccdd3e49ef118b81d437/notes/normal.pdf
# https://stats.stackexchange.com/questions/7200/evaluate-definite-interval-of-normal-distribution/7206#7206

function dlogcdf_trunc(a::Real, b::Real)
    # https://github.com/scipy/scipy/blob/a2ffe09aa751749f2372aa13c19c61b2dec5266f/scipy/stats/_continuous_distns.py
    # https://github.com/JuliaStats/Distributions.jl/blob/master/src/truncated/normal.jl
    # https://github.com/cossio/TruncatedNormal.jl/blob/master/notes/normal.pdf
    a == typemin(a) && return   normpdf(b) / normcdf(b)
    b == typemax(b) && return - normpdf(a) / normccdf(a)
    return -_tnmom1(a,b)
end

function lik_royalty(obs::ObservationRoyalty, η12::NTuple{2,Real})
    η1, η2 = η12
    l = _y(obs)
    if l == 1
        return normcdf(η2)
    elseif l == num_choices(obs)
        return normccdf(η1)
    else
        return normcdf(η2) - normcdf(η1)
    end
end

function loglik_royalty(obs::ObservationRoyalty, η12::NTuple{2,Real})
    η1, η2 = η12
    l = _y(obs)
    if l == 1
        return normlogcdf(η2)
    elseif l == num_choices(obs)
        return normlogccdf(η1)
    else
        return log(normcdf(η2) - normcdf(η1))
    end
    return F
end

function lik_loglik_royalty(obs::ObservationRoyalty, η12::NTuple{2,Real})
    l = _y(obs)
    F = lik_royalty(obs, η12)
    η1, η2 = η12
    if l == 1
        return  F, normlogcdf(η2)
    elseif l == num_choices(obs)
        return  F, normlogccdf(η1)
    else
        return F, log(F)
    end
end

# ---------------------------------------------
# Royalty likelihood
# ---------------------------------------------

"this is the main function. for each ψm, it computes LLm(roy|ψm) and ∇LLm(roy|ψm)"
function simloglik_royalty!(obs::ObservationRoyalty, theta::AbstractVector, sim::SimulationDrawsVector, dograd::Bool=true)

    am  = _am(sim)
    bm  = _bm(sim)
    cm  = _cm(sim)
    LLm = _llm(sim)
    M = _num_sim(sim)

    psi = _ψ1(sim)

    l = _y(obs)
    xbeta = _xbeta(obs)
    isfinite(xbeta) || throw(error("xbeta = $xbeta not finite"))

    check_finite(theta)

    βψ = theta_royalty_ψ(obs, theta)
    @inbounds for m in OneTo(M)
        zm  = xbeta + βψ * psi[m]
        eta12 = η12(obs, theta, zm)

        if dograd == false
            LL = loglik_royalty(obs, eta12)
            isfinite(LL) || @warn "LL = $LL not finite"
            LLm[m] += LL
        else
            η1, η2 = eta12
            F,LL  = lik_loglik_royalty(obs, eta12)
            isfinite(LL) || @warn "LL = $LL not finite. η1, η2 = $eta12 and choice =$l"
            isfinite(F) || F <= 0 || @warn "F = $F not finite or zero"
            LLm[m] += LL
            am[m] = normpdf(η1) / F
            bm[m] = normpdf(η2) / F
            cm[m] = dlogcdf_trunc(η1, η2)
        end
    end

    return nothing
end

# ---------------------------------------------
# Royalty gradient
# ---------------------------------------------

"gradient for simulated likelihood when integrating"
function grad_simloglik_royalty!(grad::AbstractVector, obs::ObservationRoyalty, theta::AbstractVector, sim::SimulationDrawsVector)

    length(grad) == length(theta) == _nparm(obs) || throw(DimensionMismatch())

    ψ1 = _ψ1(sim)
    dψ1dρ = _dψ1dθρ(sim)
    qm = _qm(sim) # Pr(um,vm | data) already computed!!!
    am = _am(sim) # computed in likelihood
    bm = _bm(sim) # computed in likelihood
    cm = _cm(sim) # computed in likelihood

    model = _model(obs)
    l = _y(obs)         # discrete choice info
    L = num_choices(obs)  # discrete choice info
    x = _x(obs)

    βψ = theta_royalty_ψ(obs, theta) # parameters

    aqm = dot(qm,am)
    bqm = dot(qm,bm)
    cqm = dot(qm,cm)

    # gradient
    grad[idx_royalty_ρ(obs)] -= sumprod3(dψ1dρ, qm, cm) * βψ
    grad[idx_royalty_ψ(obs)] -= sumprod3(ψ1,    qm, cm)
    grad[idx_royalty_β(obs)] -= cqm .* x

    if l == 1
        grad[idx_royalty_κ(obs,l)] += bqm
    elseif l < L
        idx = 2:l-1
        grad[idx_royalty_κ(obs,1)]    += cqm
        grad[idx_royalty_κ(obs,idx)] .+= cqm .* theta_royalty_κ(obs, theta, idx)
        grad[idx_royalty_κ(obs,l)]    += bqm  * theta_royalty_κ(obs, theta, l)
    elseif l == L
        idx = 2:l-1
        grad[idx_royalty_κ(obs,1)]   -= aqm
        grad[idx_royalty_κ(obs,idx)] -= aqm .* theta_royalty_κ(obs, theta, idx)
    else
        throw(DomainError(l))
    end
end

function simloglik!(grad::AbstractVector, grp::ObservationGroupRoyalty, theta, sim, dograd; kwargs...)
    simloglik_royalty!(first(grp), theta, sim, dograd)
end

function grad_simloglik!(grad, grp::ObservationGroupRoyalty, theta, sim)
    grad_simloglik_royalty!(grad, first(grp), theta, sim)
end

# ---------------------------------------------
# mini version of royalty computations for sake of testing
# ---------------------------------------------

"for testing only"
function llthreads!(grad, θ, data::DataRoyalty{<:RoyaltyModelNoHet}, dograd::Bool=true)

    model = _model(data)
    # theta_royalty_check(data, θ) || return -Inf

    k = _num_x(data)
    n = length(data)
    L = num_choices(data)
    ncoef = _nparm(data)

    ncoef == length(grad) || throw(DimensionMismatch())

    gradtmp = zeros(Float64, ncoef, n)
    LL      = Vector{Float64}(undef, n)

    beta = theta_royalty_β(data, θ)
    check_finite(beta)

    update_xbeta!(data, beta)

    for i in OneTo(n)
        gtmp = uview(gradtmp, :, i)
        LL[i] = ll_inner!(gtmp, data[i], dograd, θ)
    end
    dograd && sum!(reshape(grad, ncoef, 1), gradtmp)

    return sum(LL)
end

function ll_inner!(gradtmp::AbstractVector, grp::ObservationGroupRoyalty, dograd::Bool, θ::AbstractVector)

    RM = _model(_data(grp))
    obs = first(grp)
    eta12 = η12(obs, θ, _xbeta(obs))
    F, LL = lik_loglik_royalty(obs, eta12)
    l = _y(obs)
    L = num_choices(obs)
    theta = θ

    if dograd
        a, b = normpdf.(eta12) ./ F
        c = dlogcdf_trunc(eta12...)

        gradtmp[idx_royalty_β(obs)] .= - c .* _x(obs)
        # l > 1  && ( gradtmp[idx_royalty_κ(obs, l-1)] = -a )
        # l < L  && ( gradtmp[idx_royalty_κ(obs, l)  ] =  b )

        if l == 1
            gradtmp[idx_royalty_κ(obs,l)] += b
        elseif l < L
            idx = 2:l-1
            gradtmp[idx_royalty_κ(obs,1)]    += c
            gradtmp[idx_royalty_κ(obs,idx)] .+= c .* theta_royalty_κ(obs, theta, idx)
            gradtmp[idx_royalty_κ(obs,l)]    += b  * theta_royalty_κ(obs, theta, l)
        elseif l == L
            idx = 2:l-1
            gradtmp[idx_royalty_κ(obs,1)]   -= a
            gradtmp[idx_royalty_κ(obs,idx)] -= a .* theta_royalty_κ(obs, theta, idx)
        else
            throw(DomainError(l))
        end
    end
    return LL
end

function solve_model(data::DataRoyalty{RoyaltyModelNoHet}, theta)

    grad = similar(theta)
    f(parm) = - llthreads!(grad, parm, data, false)
    function fg!(g, parm)
        LL = llthreads!(g, parm, data, true)
        g .*= -1
        return -LL
    end

    odfg = OnceDifferentiable(f, fg!, fg!, theta)

    # check that it solves
    opts = Optim.Options(time_limit = 60, allow_f_increases=true)
    res = optimize(odfg, theta, BFGS(), opts)
    return res
end


text:
#############################################################################
#############################################################################
#
# This file contains defines unit tests for integer operations
#
#############################################################################
#############################################################################

"""
Tests Euclidean GCD algorithm for integers.
"""
function test_euclid_ints(;N::Int = 10^4)
    Random.seed!(0)
    for _ in 1:N
        n1 = rand(1:10^6)
        n2 = rand(1:10^6)
        g = euclid_alg(n1,n2)
        @assert mod(n1,g) == 0 &&  mod(n2,g) == 0
    end
    println("test_euclid_ints - PASSED")
end

"""
Tests the extended Euclidean GCD algorithm for integers.
"""
function test_ext_euclid_ints(;N::Int = 10^4)
    Random.seed!(0)
    for _ in 1:N
        n1 = rand(1:10^6)
        n2 = rand(1:10^6)
        g, s, t = ext_euclid_alg(n1,n2)
        @assert g == s*n1 + t*n2
        @assert mod(n1,g) == 0 &&  mod(n2,g) == 0
    end
    println("test_ext_euclid_ints - PASSED")
end

"""
Tests the computation of inverse mod for integers.
"""
function test_inverse_mod_ints(;prime::Int=101,N::Int=10^4)
    Random.seed!(0)
    for _ in 1:N
        n = rand(1:10^6)
        if mod(n, prime) == 0
            continue
        end
        im = int_inverse_mod(n,prime)
        @assert mod((n*im), prime) == 1
    end
    println("test_inverse_mod_ints - PASSED")
end


text:
# Autogenerated wrapper script for Sundials_jll for x86_64-linux-gnu-libgfortran5
export libsundials_arkode, libsundials_cvode, libsundials_cvodes, libsundials_ida, libsundials_idas, libsundials_kinsol, libsundials_nvecmanyvector, libsundials_nvecserial, libsundials_sunlinsolband, libsundials_sunlinsoldense, libsundials_sunlinsolklu, libsundials_sunlinsollapackband, libsundials_sunlinsollapackdense, libsundials_sunlinsolpcg, libsundials_sunlinsolspbcgs, libsundials_sunlinsolspfgmr, libsundials_sunlinsolspgmr, libsundials_sunlinsolsptfqmr, libsundials_sunmatrixband, libsundials_sunmatrixdense, libsundials_sunmatrixsparse, libsundials_sunnonlinsolfixedpoint, libsundials_sunnonlinsolnewton

using CompilerSupportLibraries_jll
using OpenBLAS_jll
using SuiteSparse_jll
JLLWrappers.@generate_wrapper_header("Sundials")
JLLWrappers.@declare_library_product(libsundials_arkode, "libsundials_arkode.so.4")
JLLWrappers.@declare_library_product(libsundials_cvode, "libsundials_cvode.so.5")
JLLWrappers.@declare_library_product(libsundials_cvodes, "libsundials_cvodes.so.5")
JLLWrappers.@declare_library_product(libsundials_ida, "libsundials_ida.so.5")
JLLWrappers.@declare_library_product(libsundials_idas, "libsundials_idas.so.4")
JLLWrappers.@declare_library_product(libsundials_kinsol, "libsundials_kinsol.so.5")
JLLWrappers.@declare_library_product(libsundials_nvecmanyvector, "libsundials_nvecmanyvector.so.5")
JLLWrappers.@declare_library_product(libsundials_nvecserial, "libsundials_nvecserial.so.5")
JLLWrappers.@declare_library_product(libsundials_sunlinsolband, "libsundials_sunlinsolband.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsoldense, "libsundials_sunlinsoldense.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsolklu, "libsundials_sunlinsolklu.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsollapackband, "libsundials_sunlinsollapackband.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsollapackdense, "libsundials_sunlinsollapackdense.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsolpcg, "libsundials_sunlinsolpcg.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspbcgs, "libsundials_sunlinsolspbcgs.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspfgmr, "libsundials_sunlinsolspfgmr.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspgmr, "libsundials_sunlinsolspgmr.so.3")
JLLWrappers.@declare_library_product(libsundials_sunlinsolsptfqmr, "libsundials_sunlinsolsptfqmr.so.3")
JLLWrappers.@declare_library_product(libsundials_sunmatrixband, "libsundials_sunmatrixband.so.3")
JLLWrappers.@declare_library_product(libsundials_sunmatrixdense, "libsundials_sunmatrixdense.so.3")
JLLWrappers.@declare_library_product(libsundials_sunmatrixsparse, "libsundials_sunmatrixsparse.so.3")
JLLWrappers.@declare_library_product(libsundials_sunnonlinsolfixedpoint, "libsundials_sunnonlinsolfixedpoint.so.2")
JLLWrappers.@declare_library_product(libsundials_sunnonlinsolnewton, "libsundials_sunnonlinsolnewton.so.2")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, OpenBLAS_jll, SuiteSparse_jll)
    JLLWrappers.@init_library_product(
        libsundials_arkode,
        "lib/libsundials_arkode.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_cvode,
        "lib/libsundials_cvode.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_cvodes,
        "lib/libsundials_cvodes.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_ida,
        "lib/libsundials_ida.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_idas,
        "lib/libsundials_idas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_kinsol,
        "lib/libsundials_kinsol.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_nvecmanyvector,
        "lib/libsundials_nvecmanyvector.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_nvecserial,
        "lib/libsundials_nvecserial.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolband,
        "lib/libsundials_sunlinsolband.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsoldense,
        "lib/libsundials_sunlinsoldense.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolklu,
        "lib/libsundials_sunlinsolklu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsollapackband,
        "lib/libsundials_sunlinsollapackband.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsollapackdense,
        "lib/libsundials_sunlinsollapackdense.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolpcg,
        "lib/libsundials_sunlinsolpcg.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspbcgs,
        "lib/libsundials_sunlinsolspbcgs.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspfgmr,
        "lib/libsundials_sunlinsolspfgmr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspgmr,
        "lib/libsundials_sunlinsolspgmr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolsptfqmr,
        "lib/libsundials_sunlinsolsptfqmr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixband,
        "lib/libsundials_sunmatrixband.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixdense,
        "lib/libsundials_sunmatrixdense.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixsparse,
        "lib/libsundials_sunmatrixsparse.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunnonlinsolfixedpoint,
        "lib/libsundials_sunnonlinsolfixedpoint.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunnonlinsolnewton,
        "lib/libsundials_sunnonlinsolnewton.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module NoiseRobustDifferentiation

using LinearAlgebra
using SparseArrays: spdiagm
using Preconditioners
using LinearMaps: LinearMap
using IterativeSolvers: cg

include("tvdiff.jl")

export tvdiff

end # module


text:
@inline _extract_val(::Val{T}) where {T} = T

"""
    Options

A general super type for all options.

# Fields
The following fields are assumed to be default. If you use different ones,
provide the access functions accordingly
* `x` a point with the current iterate
* `stop` a [`StoppingCriterion`](@ref).

"""
abstract type Options end

"""
    dispatch_options_decorator(o::Options)

Indicate internally, whether an [`Options`](@ref) `o` to be of decorating type, i.e.
it stores (encapsulates) options in itself, by default in the field `o. options`.

Decorators indicate this by returning `Val{true}` for further dispatch.

The default is `Val{false}`, i.e. by default an options is not decorated.
"""
dispatch_options_decorator(::Options) = Val(false)

"""
    is_options_decorator(o::Options)

Indicate, whether [`Options`](@ref) `o` are of decorator type.
"""
is_options_decorator(o::Options) = _extract_val(dispatch_options_decorator(o))

#
# StoppingCriterion meta
#
@doc raw"""
    StoppingCriterion

An abstract type for the functors representing stoping criteria, i.e. they are
callable structures. The naming Scheme follows functions, see for
example [`StopAfterIteration`](@ref).

Every StoppingCriterion has to provide a constructor and its function has to have
the interface `(p,o,i)` where a [`Problem`](@ref) as well as [`Options`](@ref)
and the current number of iterations are the arguments and returns a Bool whether
to stop or not.

By default each `StoppingCriterion` should provide a fiels `reason` to provide
details when a criteion is met (and that is empty otherwise).
"""
abstract type StoppingCriterion end

@doc raw"""
    StoppingCriterionGroup <: StoppingCriterion

An abstract type for a Stopping Criterion that itself consists of a set of
Stopping criteria. In total it acts as a stopping criterion itself. Examples
are [`StopWhenAny`](@ref) and [`StopWhenAll`](@ref) that can be used to
combine stopping criteria.
"""
abstract type StoppingCriterionSet <: StoppingCriterion end
#
#
# StepsizeOptions
#
#
"""
    Stepsize

An abstract type for the functors representing step sizes, i.e. they are callable
structurs. The naming scheme is `TypeOfStepSize`, e.g. `ConstantStepsize`.

Every Stepsize has to provide a constructor and its function has to have
the interface `(p,o,i)` where a [`Problem`](@ref) as well as [`Options`](@ref)
and the current number of iterations are the arguments
and returns a number, namely the stepsize to use.

# See also
[`Linesearch`](@ref)
"""
abstract type Stepsize end

@doc raw"""
    get_options(o::Options)

return the undecorated [`Options`](@ref) of the (possibly) decorated `o`.
As long as your decorated options store the options within `o.options` and
the [`dispatch_options_decorator`](@ref) is set to `Val{true}`,
the internal options are extracted.
"""
get_options(o::Options) = get_options(o, dispatch_options_decorator(o))
get_options(o::Options, ::Val{false}) = o
get_options(o::Options, ::Val{true}) = get_options(o.options)

@doc raw"""
    get_reason(o)

return the current reason stored within the [`StoppingCriterion`](@ref) from
within the [`Options`](@ref) This reason is empty if the criterion has never
been met.
"""
get_reason(o::Options) = get_reason(get_options(o).stop)

#
# Common Actions for decorated Options
#
@doc raw"""
    AbstractOptionsAction

a common `Type` for `AbstractOptionsActions` that might be triggered in decoraters,
for example [`DebugOptions`](@ref) or [`RecordOptions`](@ref).
"""
abstract type AbstractOptionsAction end

@doc raw"""
    StoreTupleAction <: AbstractOptionsAction

internal storage for [`AbstractOptionsAction`](@ref)s to store a tuple of fields from an
[`Options`](@ref)s

This functor posesses the usual interface of functions called during an
iteration, i.e. acts on `(p,o,i)`, where `p` is a [`Problem`](@ref),
`o` is an [`Options`](@ref) and `i` is the current iteration.

# Fields
* `values` – a dictionary to store interims values based on certain `Symbols`
* `keys` – an `NTuple` of `Symbols` to refer to fields of `Options`
* `once` – whether to update the internal values only once per iteration
* `lastStored` – last iterate, where this `AbstractOptionsAction` was called (to determine `once`

# Constructiors

    StoreOptionsAction([keys=(), once=true])

Initialize the Functor to an (empty) set of keys, where `once` determines
whether more that one update per iteration are effective

    StoreOptionsAction(keys, once=true])

Initialize the Functor to a set of keys, where the dictionary is initialized to
be empty. Further, `once` determines whether more that one update per iteration
are effective, otherwise only the first update is stored, all others are ignored.
"""
mutable struct StoreOptionsAction <: AbstractOptionsAction
    values::Dict{Symbol,<:Any}
    keys::NTuple{N,Symbol} where {N}
    once::Bool
    last_stored::Int
    function StoreOptionsAction(
        keys::NTuple{N,Symbol} where {N}=NTuple{0,Symbol}(), once=true
    )
        return new(Dict{Symbol,Any}(), keys, once, -1)
    end
end
function (a::StoreOptionsAction)(::P, o::O, i::Int) where {P<:Problem,O<:Options}
    #update values (maybe only once)
    if !a.once || a.last_stored != i
        merge!(
            a.values,
            Dict{Symbol,Any}(key => deepcopy(getproperty(o, key)) for key in a.keys),
        )
    end
    return a.last_stored = i
end

"""
    get_storage(a,key)

return the internal value of the [`StoreOptionsAction`](@ref) `a` at the
`Symbol` `key`.
"""
get_storage(a::StoreOptionsAction, key) = a.values[key]

"""
    get_storage(a,key)

return whether the [`StoreOptionsAction`](@ref) `a` has a value stored at the
`Symbol` `key`.
"""
has_storage(a::StoreOptionsAction, key) = haskey(a.values, key)

"""
    update_storage!(a,o)

update the [`StoreOptionsAction`](@ref) `a` internal values to the ones given on
the [`Options`](@ref) `o`.
"""
function update_storage!(a::StoreOptionsAction, o::O) where {O<:Options}
    return update_storage!(a, Dict(key => getproperty(o, key) for key in a.keys))
end

"""
    update_storage!(a,o)

update the [`StoreOptionsAction`](@ref) `a` internal values to the ones given in
the dictionary `d`. The values are merged, where the values from `d` are preferred.
"""
function update_storage!(a::StoreOptionsAction, d::Dict{Symbol,<:Any})
    merge!(a.values, d)
    # update keys
    return a.keys = Tuple(keys(a.values))
end


text:
import Base.|>

|>(x::Tuple, f) = f(x...)

|>(x, f) = f(x)

unzip(a) = map(x -> getfield.(a, x), fieldnames(eltype(a)))

max_stars_repo_name: manuelelucchi/Fibonacci

max_issues_repo_name: manuelelucchi/Fibonacci

max_forks_repo_name: manuelelucchi/Fibonacci


text:
module TimeSeries

using Dates
using DelimitedFiles
using Statistics

using RecipesBase
using Reexport

export TimeArray, AbstractTimeSeries,
       when, from, to, findwhen, find, timestamp, values, colnames, meta, head, tail,
       lag, lead, diff, percentchange, moving, upto,
       uniformspaced, uniformspace, dropnan,
       basecall,
       merge, collapse,
       readtimearray, writetimearray,
       update, rename, rename!

@reexport using Dates

###### include ##################

include(".timeseriesrc.jl")
include("timearray.jl")
include("utilities.jl")
include("split.jl")
include("apply.jl")
include("broadcast.jl")
include("combine.jl")
include("readwrite.jl")
include("modify.jl")
include("basemisc.jl")
include("deprecated.jl")
include("plotrecipes.jl")

end  # module TimeSeries


text:
"""
Holder of Shader context info, can also be used for dispatch
"""
abstract type AbstractContext end

"""
A Sampler, that supports interpolated access
"""
abstract type AbstractSampler{T, N} <: DenseArray{T, N} end

abstract type AbstractSamplerBuffer{T} <: DenseVector{T} end

"""
VertexArray, holds the vertex info a vertex shaders maps over.
"""
abstract type AbstractVertexArray{T} <: DenseVector{T} end

struct Sampler{T, N} <: AbstractSampler{T, N}
    data::T
    minfilter::Symbol
    magfilter::Symbol # magnification
    repeat::NTuple{N, Symbol}
    anisotropic::Float32
    color_swizzel::Vector{Symbol}
end

function Sampler(
        data::AbstractArray{T, N};
        minfilter = T <: Integer ? :nearest : :linear,
        magfilter = minfilter, # magnification
        x_repeat  = :clamp_to_edge, #wrap_s
        y_repeat  = x_repeat, #wrap_t
        z_repeat  = x_repeat, #wrap_r
        anisotropic = 1f0,
        color_swizzel = nothing
    ) where {T, N}

    swizzel = color_swizzel !== nothing ? color_swizzel : if T <: Gray
        Symbol[:RED, :RED, :RED, :ONE]
    elseif T <: GrayA
        Symbol[:RED,:_RED, :RED, :ALPHA]
    else
        Symbol[]
    end
    Sampler{T, N}(
        data, minfilter, magfilter,
        (x_reapt, y_repeat, z_repeat),
        anisotropic, swizzel
    )
end

struct BufferSampler{T} <: AbstractSamplerBuffer{T}
    data::T
end

struct VertexArray{ET, Data} <: AbstractVertexArray{ET}
    data::Data
end

Tables.schema(va::VertexArray) = Tables.schema(getfield(va, :data))

# GeometryBasics.column_names(va::VertexArray) = GeometryBasics.column_names(getfield(va, :data))
# GeometryBasics.column_types(va::VertexArray) = GeometryBasics.column_types(getfield(va, :data))

Base.size(x::VertexArray) = size(x.data)
Base.getindex(x::VertexArray, i) = getindex(x.data, i)

function VertexArray(data::AbstractArray{T}) where T
    return VertexArray{T, typeof(data)}(data)
end

function VertexArray(points, faces = nothing; kw_args...)
    vertices = if isempty(kw_args)
        points
    else
        GeometryBasics.meta(points; kw_args...)
    end
    data = if faces === nothing
        vertices
    else
        GeometryBasics.connect(vertices, faces)
    end
    VertexArray(data)
end

function VertexArray(mesh::GeometryTypes.AbstractMesh)
    data = GeometryTypes.attributes(mesh)
    vs = GeometryBasics.Point.(pop!(data, :vertices))
    fs = convert(Vector{GeometryBasics.TriangleFace{Cuint}}, pop!(data, :faces))
    m = GeometryBasics.Mesh(GeometryBasics.meta(vs; data...), fs)
    VertexArray{eltype(m), typeof(m)}(m)
end

function VertexArray(; meta...)
    m = StructArray(; meta...)
    VertexArray{eltype(m), typeof(m)}(m)
end


text:
# functions for I/O operations of data and results

doc"""
    export_data(problem::Problem, results::Result, jobname::String)

Exportes the resultes into three cvs files

# Files
* `output_res.csv`: contines the estimations of Theta, loglik and the accept vec
* `output_prior_dist.csv`: contines parameters for the prior distributions.
* `output_prior_dist_type.csv`: contines information reg. the prior distribution.
* `data_used.csv`: contins the data set.
"""
function export_data(problem::Problem, results::Result, jobname::String)

  # set parameters
  theta_true = problem.model_param.theta_true
  burn_in = problem.alg_param.burn_in
  prior_parameters = problem.prior_dist.prior_parameters
  dist_type = problem.prior_dist.dist
  Theta_est = results.Theta_est'
  loglik_est = results.loglik_est
  accept_vec_est = results.accept_vec

  # export results
  data_param = zeros(1, length(theta_true)+1)
  data_param[1:length(theta_true)] = theta_true
  data_param[length(theta_true)+1:length(theta_true)+1] = burn_in
  writetable("./Results/output_param"*jobname*".csv", convert(DataFrame, data_param))

  data_output = zeros(length(loglik_est), length(theta_true) + 2)
  data_output[:,1:length(theta_true)] = Theta_est
  data_output[:,length(theta_true)+1:length(theta_true)+1] = loglik_est
  data_output[:,end] = accept_vec_est

  #
  data_used = zeros(length(problem.data.Z),2)
  data_used[:,1] = problem.data.Z
  writetable("./Results/output_res"*jobname*".csv", convert(DataFrame, data_output))

  writetable("./Results/output_prior_dist"*jobname*".csv", convert(DataFrame, prior_parameters))

  writetable("./Results/output_prior_dist_type"*jobname*".csv", convert(DataFrame, [1 dist_type]))

  writetable("./Results/data_used"*jobname*".csv", convert(DataFrame,data_used))

end

doc"""
    export_data(problem::gpProblem, results::Result, jobname::String)

Exportes the resultes into three cvs files

# Files
* `output_res.csv`: contines the estimations of Theta, loglik and the accept vec
* `output_prior_dist.csv`: contines parameters for the prior distributions.
* `output_prior_dist_type.csv`: contines information reg. the prior distribution.
* `data_used.csv`: contins the data set.
"""
function export_data(problem::gpProblem, results::Result, jobname::String)

  # set parameters
  theta_true = problem.model_param.theta_true
  burn_in = problem.alg_param.burn_in
  length_training_data = problem.alg_param.length_training_data
  prior_parameters = problem.prior_dist.prior_parameters
  dist_type = problem.prior_dist.dist
  Theta_est = results.Theta_est'
  loglik_est = results.loglik_est
  accept_vec_est = results.accept_vec

  # export results
  data_param = zeros(1, length(theta_true)+1)
  data_param[1:length(theta_true)] = theta_true
  data_param[length(theta_true)+1:length(theta_true)+1] = burn_in + length_training_data
  writetable("./Results/output_param_dagp"*jobname*".csv", convert(DataFrame, data_param))

  data_output = zeros(length(loglik_est), length(theta_true) + 2)
  data_output[:,1:length(theta_true)] = Theta_est
  data_output[:,length(theta_true)+1:length(theta_true)+1] = loglik_est
  data_output[:,end] = accept_vec_est

  #
  data_used = zeros(length(problem.data.Z),2)
  data_used[:,1] = problem.data.Z
  writetable("./Results/output_res_dagp"*jobname*".csv", convert(DataFrame, data_output))

  writetable("./Results/output_prior_dist_dagp"*jobname*".csv", convert(DataFrame, prior_parameters))

  writetable("./Results/output_prior_dist_type_dagp"*jobname*".csv", convert(DataFrame, [1 dist_type]))

  writetable("./Results/data_used_dagp"*jobname*".csv", convert(DataFrame,data_used))

end

doc"""
    export_parameters(lamba_t::Array)

Exportes the parmeters for the AM gen/AM gen cw algorithm into a cvs files

# Files
* `output_time_dep_scaling.csv`: contines the parameters values .

"""
function export_parameters(lamba_t::Array)

  writetable("output_time_dep_scaling.csv", convert(DataFrame, lamba_t))

end

doc"""
    load_data()

Loads the protien folding data into an array

# Outputs
* `Z_data`: array containing the data.
"""
function load_data()
  file = open("1LE1_L.dat")
  data = readlines(file)
  close(file)

  Z_data = zeros(1,length(data))

  for i = 1:length(data)
    Z_data[i] = readdlm(IOBuffer(data[i]), Float64)[2]
  end
  return Z_data
end

#=

doc"""
    export_pf_diagnostics(loglik::Array, weigths::Array, particels::Array, Z::Array, X::Array)

Exports data realated to the diagnostics res. for the pf to the files; `output_pf_diagnistics_loglik.csv`,
`output_pf_diagnistics_particles.csv`,`output_pf_diagnistics_particles.csv`, and `output_pf_diagnistics_Z_X.csv`.
"""
function export_pf_diagnostics(loglik::Array, weigths::Array, particels::Array, Z::Array, X::Array)

  writetable("output_pf_diagnistics_loglik.csv", convert(DataFrame, loglik))
  writetable("output_pf_diagnistics_weigths.csv", convert(DataFrame, weigths))
  writetable("output_pf_diagnistics_particles.csv", convert(DataFrame, particels))

  data_Z_X = zeros(length(Z), 2)
  data_Z_X[:,1] = Z'
  data_Z_X[:,2] = X'

  writetable("output_pf_diagnistics_Z_X.csv", convert(DataFrame, data_Z_X))

end

doc"""
    export_pf_diagnostics(loglik::Array, weigths::Array, particels::Array, Z::Array, X::Array)

Exports data realated to the diagnostics res. for the pf to the files; `output_pf_diagnistics_loglik.csv`,
`output_pf_diagnistics_particles.csv`,`output_pf_diagnistics_particles.csv`, and `output_pf_diagnistics_Z_X.csv`.
"""
function export_pf_diagnostics(loglik::RowVector, weigths::Array, particels::Array, Z::Array, X::Array)

  writetable("output_pf_diagnistics_loglik.csv", convert(DataFrame, loglik))
  writetable("output_pf_diagnistics_weigths.csv", convert(DataFrame, weigths))
  writetable("output_pf_diagnistics_particles.csv", convert(DataFrame, particels))

  data_Z_X = zeros(length(Z), 2)
  data_Z_X[:,1] = Z'
  data_Z_X[:,2] = X'

  writetable("output_pf_diagnistics_Z_X.csv", convert(DataFrame, data_Z_X))

end

=#

max_stars_repo_path: DWPSDE model/io_functions.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: DWPSDE model/io_functions.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: DWPSDE model/io_functions.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
# -------------------------------------------------------------------
# JULIA TOOL FOR EXPORTING INTERACTION MATRICES AS C++-HEADER FILES
# AS USED IN KITAEV_QMC AND KITAEV_QMC_KPM
# Author: Tim Eschmann
# Date: 2017-01-26
# Modified for new version of LatticePhysics: 2019-01-31
#
# This tool requires the Julia package "LatticePhysics.jl"
# (Copyright (c) 2017: Jan Attig) which is available under
# https://github.com/janattig/LatticePhysics.jl
# -------------------------------------------------------------------

include("matrix_creator.jl")

# use the modules
using LatticePhysics
using MatrixCreator

#####################################################################
# 1. Define System:

# Lattice:
lat = "(10,3)a"

# Size:
Lx = 8
Ly = 8
Lz = 8

# Periodicity:
per_dir = "periodic"
whichone = "Any"

# Gauge field:
gauge = "All"

# Loop length:
l = 10

# Coupling Constants
Jx = 1/3.
Jy = 1/3.
Jz = 1/3.

#####################################################################
# 2. Choose unit cell:

if (lat == "(6,3)a")
	uc = getUnitcellHoneycomb(4)
	filename = "honeycomb_$(Lx)_$(Ly)_Jz_$(Jz)_disorder.hpp"
elseif (lat == "(8,3)a")
	uc = getUnitcell_8_3_a(4)
	filename = "8_3_a_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
elseif (lat == "(8,3)b")
	uc = getUnitcell_8_3_b(4)
	filename = "8_3_b_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
elseif (lat == "(8,3)c")
	uc = getUnitcell_8_3_c(4)
	filename = "8_3_c_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
elseif (lat == "(8,3)n")
	uc = getUnitcell_8_3_n(4)
	filename = "8_3_n_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
elseif (lat == "(9,3)a")
	uc = getUnitcell_9_3_a(4)
	filename = "9_3_a_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
elseif (lat == "(10,3)a")
	uc = getUnitcell_10_3_a(4)
	filename = "10_3_a_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
elseif (lat == "(10,3)b")
	uc = getUnitcell_10_3_b(4)
	filename = "10_3_b_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
elseif (lat == "(10,3)c")
	uc = getUnitcell_10_3_c(4)
	filename = "10_3_c_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
elseif (lat == "(10,3)d")
	uc = getUnitCell_10_3_d(4)
	filename = "10_3_d_$(Lx)x$(Ly)x$(Lz)_$(per_dir)_$(l)flux_Jz_$(Jz).hpp"
else
	println("ERROR: Choose unitcell")
end

#####################################################################
# 3. Generate Lattice:

# Boundary conditions:
if (lat == "(6,3)a" && per_dir == "periodic")
	lt = getLatticePeriodic2D(uc, [Lx, Ly])
elseif (lat == "(6,3)a" && per_dir == "open")
	lt = getLatticeOpen2D(uc, [Lx, Ly])
elseif (per_dir == "open")
	lt = getLatticeOpen3D(uc, [Lx, Ly, Lz])
#elseif (per_dir == "poo")
#	lt = getStripLattice3D(uc, -Lx, Ly, Lz, whichone)
#elseif (per_dir == "opo")
#	lt = getStripLattice3D(uc, Lx, -Ly, Lz, whichone)
#elseif (per_dir == "oop")
#	lt = getStripLattice3D(uc, Lx, Ly, -Lz, whichone)
#elseif (per_dir == "ppo")
#	lt = getStripLattice3D(uc, -Lx, -Ly, Lz, whichone)
#elseif (per_dir == "pop")
#	lt = getStripLattice3D(uc, -Lx, Ly, -Lz, whichone)
#elseif (per_dir == "opp")
#	lt = getStripLattice3D(uc, Lx, -Ly, -Lz, whichone)
elseif (per_dir == "periodic")
	lt = getLatticePeriodic3D(uc, [Lx,Ly,Lz])
end

#####################################################################
# 4. Create Kitaev interaction matrix A from lattice:

matrix = getInteractionMatrixCouplings(lt, Jx, Jy, Jz, gauge)

# Disordered version:
#p = 0.25
#dJ = 0.8
#matrix = getInteractionMatrixDisorder(lt, Jx, Jy, Jz, dJ, p)

#####################################################################
# 5. Find closed (even-length) loops inside the lattice:
netto_plaquettes = getPlaquettesOfLattice(lt,l)

#####################################################################
# 6. Write Header-File for C++ Monte Carlo Program:

writehpp(lt, netto_plaquettes, l, matrix, filename)

# Include only zz-bonds
writehpp_zz(lt, netto_plaquettes, l, matrix, filename)
#####################################################################

exit()

max_stars_repo_name: timeschmann/LatticeFilesForKitaevQMC

max_issues_repo_name: timeschmann/LatticeFilesForKitaevQMC

max_forks_repo_name: timeschmann/LatticeFilesForKitaevQMC


text:
module SmolyakApprox

using LinearAlgebra
using ThreadPools

include("chebyshev_gauss_lobatto.jl")
include("clenshaw_curtis_equidistant.jl")
include("chebyshev_polynomial.jl")
include("smolyak_weights.jl")
include("smolyak_evaluate.jl")
include("smolyak_grid.jl")
include("generate_multi_index.jl")
include("combine_nodes.jl")
include("m_i.jl")
include("chebyshev_polynomial_derivative.jl")
include("smolyak_derivative_finite_difference.jl")
include("smolyak_derivative.jl")
include("scale_nodes.jl")
include("smolyak_full_grid.jl")
include("smolyak_piecewise_linear.jl")

export chebyshev_gauss_lobatto,
       clenshaw_curtis_equidistant,
       smolyak_grid,
       smolyak_weights,
       smolyak_weights_threaded,
       smolyak_inverse_interpolation_matrix,
       smolyak_inverse_interpolation_matrix_threaded,
       smolyak_evaluate,
       smolyak_derivative_finite_difference,
       smolyak_derivative,
       smolyak_gradient,
       smolyak_grid_full,
       smolyak_weights_full,
       smolyak_evaluate_full,
       smolyak_derivative_full,
       smolyak_gradient_full,
       smolyak_pl_weights,
       smolyak_pl_weights_threaded,
       smolyak_pl_evaluate

end


text:
############################
## read and write to disk ##
############################

function readTimedata(io::IO)

    df = readtable(io)
    td = convert(AbstractTimedata, df)
end

function readTimedata(filename::String)

    df = readtable(filename, nastrings=[".", "", "NA"])
    td = convert(AbstractTimedata, df)
end

function writeTimedata(filename::String, td::AbstractTimedata)
    ## create large dataframe
    idxDf = DataFrame(idx = idx(td));
    df = [idxDf td.vals];
    writetable(filename, df)
end

######################
## org-babel output ##
######################

function Base.writedlm(io::IO, td::AbstractTimedata, dlm; opts...)
    nams = [:idx, names(td)]
    vals = [idx(td) arrAny = get(td)]
    arrAny = [nams', vals]
    writedlm(io, arrAny, dlm; opts...)
end


text:
################################################################################
# Data Jacobians
################################################################################
function joint_constraint_jacobian_body_data(mechanism::Mechanism, joint::JointConstraint{T,N}, body::Body{T}) where {T,N}
    Nd = data_dim(body)
    ∇m = szeros(T,N,1)
    ∇J = szeros(T,N,6)
    ∇v15 = szeros(T,N,3)
    ∇ϕ15 = szeros(T,N,3)
    ∇z2 = -constraint_jacobian_configuration(mechanism, joint, body) *
        integrator_jacobian_configuration(body, mechanism.timestep, attjac=true)
    ∇g = [∇m ∇J ∇v15 ∇ϕ15 ∇z2]
    return ∇g
end

function body_constraint_jacobian_body_data(mechanism::Mechanism, body::Body{T}) where T
    Δt = mechanism.timestep
    N = 6
    x1, v15, q1, ϕ15 = previous_configuration_velocity(body.state)
    x2, v25, q2, ϕ25 = current_configuration_velocity(body.state)
    x3, q3 = next_configuration(body.state, Δt)
    # Mass
    gravity=mechanism.gravity
    ∇m = [1 / Δt * (x2 - x1) + Δt/2 * gravity - 1 / Δt * (x3 - x2) + Δt/2 * gravity;
          szeros(T,3,1)]
    # Inertia
    ∇J = 2 / Δt * LVᵀmat(q2)' * LVᵀmat(q1) * ∂Jp∂J(VLᵀmat(q1) * vector(q2))
    ∇J += 2 / Δt * LVᵀmat(q2)' * Tmat() * RᵀVᵀmat(q3) * ∂Jp∂J(VLᵀmat(q2) * vector(q3))
    ∇J = [szeros(T,3,6); ∇J]

    # initial conditions: v15, ϕ15
    ∇v15 = body.mass * SMatrix{3,3,T,9}(Diagonal(sones(T,3)))
    ∇q1 = -2 / Δt * LVᵀmat(q2)' * ∂LVᵀmat∂q(body.inertia * VLᵀmat(q1) * vector(q2))
    ∇q1 += -2 / Δt * LVᵀmat(q2)' * LVᵀmat(q1) * body.inertia * ∂VLᵀmat∂q(vector(q2))
    ∇ϕ15 = ∇q1 * rotational_integrator_jacobian_velocity(q2, -ϕ15, Δt)
    ∇15 = [∇v15 szeros(T,3,3);
           szeros(T,3,3) ∇ϕ15]

    # current configuration: z2 = x2, q2
    # manipulator's equation contribution
    ∇tra_x2 = - 2 / Δt * body.mass * SMatrix{3,3,T,9}(Diagonal(sones(T,3)))
    ∇tra_q2 = szeros(T,3,3)
    ∇rot_x2 = szeros(T,3,3)
    ∇rot_q2 = -2 / Δt * VLᵀmat(q2) * LVᵀmat(q1) * body.inertia * VLᵀmat(q1)
    ∇rot_q2 += -2 / Δt * VLᵀmat(q2) * Tmat() * RᵀVᵀmat(q3) * body.inertia * ∂VLᵀmat∂q(vector(q3))
    ∇rot_q2 += -2 / Δt * ∂VLᵀmat∂q(LVᵀmat(q1) * body.inertia * VLᵀmat(q1) * vector(q2) + Tmat() * RᵀVᵀmat(q3) * body.inertia * VLᵀmat(q2) * vector(q3))
    ∇rot_q2 *= LVᵀmat(q2)
    ∇z2 = [∇tra_x2 ∇tra_q2;
           ∇rot_x2 ∇rot_q2]
    # @show ∇z2
    # TODO
    # # constact constraints impulses contribution
    # @warn "000"
    return [∇m ∇J ∇15 0.0000000000*∇z2] #TODO not sure why we need to zero out this block, maybe finite diff is not correct and we try to match finite diff
end

function body_constraint_jacobian_body_data(mechanism::Mechanism, pbody::Node{T},
        cbody::Node{T}, joint::JointConstraint{T,N,Nc}) where {T,N,Nc}
    # Jacobian of pbody's dynamics constraints wrt cbody's data (x2b, q2b)
    # this comes from the fact that the Joint Constraint force mapping of pbody
    # depends on cbody's data (x2b, q2b)
    # This is the same for spring and damper forces.
    timestep= mechanism.timestep

    ∇z2_aa = szeros(T,6,6)
    ∇z2_ab = szeros(T,6,6)
    # joint constraints impulses contribution
    for i = 1:Nc
        λ = get_joint_impulses(joint, i)
        if cbody.id == joint.child_id
            ∇z2_aa += impulse_map_jacobian(:parent, :parent, [joint.translational, joint.rotational][i],
                pbody, cbody, λ)
            ∇z2_ab += impulse_map_jacobian(:parent, :child, [joint.translational, joint.rotational][i],
                pbody, cbody, λ)
        elseif pbody.id == joint.child_id
            ∇z2_aa += impulse_map_jacobian(:child, :child, [joint.translational, joint.rotational][i],
                cbody, pbody, λ)
            ∇z2_ab += impulse_map_jacobian(:child, :parent, [joint.translational, joint.rotational][i],
                cbody, pbody, λ)
        end
    end
    # spring and damper impulses contribution
    if joint.spring
        for i = 1:Nc
            λ = get_joint_impulses(joint, i)
            if cbody.id == joint.child_id
                ∇z2_aa += spring_jacobian_configuration(
                    :parent, :parent,
                    [joint.translational, joint.rotational][i], pbody, cbody, timestep)
                ∇z2_ab += spring_jacobian_configuration(
                    :parent, :child,
                    [joint.translational, joint.rotational][i], pbody, cbody, timestep)
            elseif pbody.id == joint.child_id
                ∇z2_aa += spring_jacobian_configuration(
                    :child, :child,
                    [joint.translational, joint.rotational][i], cbody, pbody, timestep)
                ∇z2_ab += spring_jacobian_configuration(
                    :child, :parent,
                    [joint.translational, joint.rotational][i], cbody, pbody, timestep)
            end
        end
    end
    if joint.damper
        for i = 1:Nc
            λ = get_joint_impulses(joint, i)
            if cbody.id == joint.child_id
                ∇z2_aa += damper_jacobian_configuration(
                    :parent, :parent,
                    [joint.translational, joint.rotational][i], pbody, cbody, timestep)
                ∇z2_ab += damper_jacobian_configuration(
                    :parent, :child,
                    [joint.translational, joint.rotational][i], pbody, cbody, timestep)
            elseif pbody.id == joint.child_id
                ∇z2_aa += damper_jacobian_configuration(
                    :child, :child,
                    [joint.translational, joint.rotational][i], cbody, pbody, timestep)
                ∇z2_ab += damper_jacobian_configuration(
                    :child, :parent,
                    [joint.translational, joint.rotational][i], cbody, pbody, timestep)
            end
        end
    end
    return [szeros(T,6,13) ∇z2_aa], [szeros(T,6,13) ∇z2_ab]
end

function body_constraint_jacobian_body_data(mechanism::Mechanism, body::Node{T},
        contact::ContactConstraint{T,N,Nc}) where {T,N,Nc}
    # Jacobian of the Body's dynamics constraints wrt the Body's data (x2, q2)
    # this comes from the fact that the Contact Constraint force mapping depends
    # on the Body's data (x2, q2)
    # contact constraints impulses contribution
    ∇z3 = impulse_map_jacobian_configuration(mechanism, body, contact)
    ∇z2 = ∇z3 * integrator_jacobian_configuration(body, mechanism.timestep)
    return [szeros(T,6,13) ∇z2]
end

function body_constraint_jacobian_joint_data(mechanism::Mechanism{T}, body::Body{T},
        joint::JointConstraint{T}) where {T}
    Δt = mechanism.timestep
    Nd = data_dim(joint)
    N = 6
    x1, v15, q1, ϕ15 = previous_configuration_velocity(body.state)
    x2, v25, q2, ϕ25 = current_configuration_velocity(body.state)
    x3, q3 = next_configuration(body.state, Δt)
    # ∇u = Diagonal(SVector{6,T}(1,1,1,2,2,2)) * input_jacobian_control(mechanism, joint, body)
    ∇u = Diagonal(SVector{6,T}(1,1,1,1,1,1)) * input_jacobian_control(mechanism, joint, body)
    ∇spring = joint.spring ? spring_impulses(mechanism, joint, body, unitary=true) : szeros(T,6,1)
    ∇damper = joint.damper ? damper_impulses(mechanism, joint, body, unitary=true) : szeros(T,6,1)
    return [∇u ∇spring ∇damper]
end

function body_constraint_jacobian_contact_data(mechanism::Mechanism, body::Body{T},
        contact::ContactConstraint{T,N,Nc,Cs,N½}) where {T,N,Nc,Cs<:NonlinearContact{T,N},N½}
    Nd = data_dim(contact)
    model = contact.model
    offset = model.offset
    x3, q3 = next_configuration(body.state, mechanism.timestep)
    γ = contact.impulses[2]

    ∇friction_coefficient = szeros(T,3,1)

    X = force_mapping(model, x3, q3)
    # this what we differentiate: Qᵀγ = - skew(p - vector_rotate(offset, inv(q3))) * VRmat(q3) * LᵀVᵀmat(q3) * X' * γ
    ∇p = - ∂skew∂p(VRmat(q3) * LᵀVᵀmat(q3) * X' * γ)
    ∇off = - ∂skew∂p(VRmat(q3) * LᵀVᵀmat(q3) * X' * γ) * -∂vector_rotate∂p(offset, inv(q3))

    ∇X = szeros(T,3,Nd)
    ∇Q = -[∇friction_coefficient ∇off ∇p]
    return [∇X; ∇Q]
end

function contact_constraint_jacobian_contact_data(mechanism::Mechanism, contact::ContactConstraint{T,N,Nc,Cs,N½}, body::Body{T}) where {T,N,Nc,Cs<:NonlinearContact{T,N},N½}
    Nd = data_dim(contact)
    model = contact.model
    p = model.contact_point
    offset = model.offset
    x2, v25, q2, ϕ25 = current_configuration_velocity(body.state)
    x3, q3 = next_configuration(body.state, mechanism.timestep)
    s = contact.impulses_dual[2]
    γ = contact.impulses[2]

    ∇friction_coefficient = SA[0,γ[1],0,0]
    ∇off = [-model.surface_normal_projector; szeros(T,1,3); -model.surface_projector * skew(vector_rotate(ϕ25, q3))]
    ∇p = [model.surface_normal_projector * ∂vector_rotate∂p(model.contact_point, q3); szeros(T,1,3); model.surface_projector * skew(vector_rotate(ϕ25, q3)) * ∂vector_rotate∂p(model.contact_point, q3)]

    ∇compμ = szeros(T,N½,Nd)
    ∇g = -[∇friction_coefficient ∇off ∇p]
    return [∇compμ; ∇g]
end

function contact_constraint_jacobian_body_data(mechanism::Mechanism, contact::ContactConstraint{T,N,Nc,Cs,N½}, body::Body{T}) where {T,N,Nc,Cs,N½}
    Nd = data_dim(body)
    ∇compμ = szeros(T,N½,Nd)
    ∇m = szeros(T,N½,1)
    ∇J = szeros(T,N½,6)
    ∇v15 = szeros(T,N½,3)
    ∇ϕ15 = szeros(T,N½,3)
    ∇z3 = - constraint_jacobian_configuration(mechanism, contact, body) # minus sign coming from res = [-compμ; -constraint]
    ∇z2 = ∇z3 * integrator_jacobian_configuration(body, mechanism.timestep, attjac=true) # 4x7 * 7x6 = 4x6
    ∇g = [∇m ∇J ∇v15 ∇ϕ15 ∇z2]
    return [∇compμ; ∇g]
end

################################################################################
# System Data Jacobians
################################################################################
function data_adjacency_matrix(joints::Vector{<:JointConstraint}, bodies::Vector{<:Body},
        contacts::Vector{<:ContactConstraint})
    # mode can be impulses or data depending on whi
    nodes = [joints; bodies; contacts]
    n = length(nodes)
    A = zeros(Bool, n, n)

    for node1 in nodes
        for node2 in nodes
            T1 = typeof(node1)
            T2 = typeof(node2)
            if T1 <: Body
                if T2 <: Body
                    (node1.id == node2.id) && (A[node1.id, node2.id] = 1) # self loop
                    linked = length(indirect_link(node1.id, node2.id, [joints; contacts])) > 0
                    linked && (A[node1.id, node2.id] = 1) # linked through a common joint
                elseif T2 <: JointConstraint
                    (node1.id == node2.parent_id || node1.id == node2.child_id) && (A[node1.id, node2.id] = 1) # linked
                elseif T2 <: ContactConstraint
                    (node1.id == node2.parent_id || node1.id == node2.child_id) && (A[node1.id, node2.id] = 1) # linked
                end
            elseif T1 <: JointConstraint
                if T2 <: Body
                    (node2.id == node1.parent_id || node2.id == node1.child_id) && (A[node1.id, node2.id] = 1) # linked
                end
            elseif T1 <: ContactConstraint
                if T2 <: Body
                    (node2.id == node1.parent_id || node2.id == node1.child_id) && (A[node1.id, node2.id] = 1) # linked
                elseif T2 <: ContactConstraint
                    (node1.id == node2.id) && (A[node1.id, node2.id] = 1) # self loop
                end
            end
        end
    end
    A = convert(Matrix{Int64}, A)
    return A
end

function indirect_link(id1, id2, nodes::Vector{S}) where {S<:Node}
    ids = zeros(Int, 0)
    for node in nodes
        parent_id = node.parent_id
        (parent_id == nothing) && (parent_id = 0) #handle the origin's corner case
        linked = (id1 == node.child_id) && (id2 == parent_id)
        linked |= (id2 == node.child_id) && (id1 == parent_id)
        linked && push!(ids, node.id)
    end
    return ids
end

function create_data_matrix(joints::Vector{<:JointConstraint}, bodies::Vector{B}, contacts::Vector{<:ContactConstraint};
        force_static::Bool=false) where {T,B<:Body{T}}
    nodes = [joints; bodies; contacts]
    A = data_adjacency_matrix(joints, bodies, contacts)
    dimrow = length.(nodes)
    dimcol = data_dim.(nodes)

    N = length(dimrow)
    static = force_static || (all(dimrow.<=10) && all(dimcol.<=10))
    data_matrix = spzeros(Entry,N,N)

    for i = 1:N
        for j = 1:N
            if A[i,j] == 1
                data_matrix[i,j] = Entry{T}(dimrow[i], dimcol[j], static = static)
            end
        end
    end
    return data_matrix
end

function jacobian_data!(data_matrix::SparseMatrixCSC, mechanism::Mechanism)
    jacobian_contact_data!(data_matrix, mechanism)
    jacobian_body_data!(data_matrix, mechanism)
    jacobian_joint_data!(data_matrix, mechanism)
    return nothing
end

function jacobian_contact_data!(data_matrix::SparseMatrixCSC, mechanism::Mechanism{T}) where {T}
    # ∂body∂ineqcdata
    for contact in mechanism.contacts
        pbody = get_body(mechanism, contact.parent_id)
        data_matrix[pbody.id, contact.id].value += body_constraint_jacobian_contact_data(mechanism, pbody, contact)
    end
    # ∂contact∂contactdata
    for contact in mechanism.contacts
        pbody = get_body(mechanism, contact.parent_id)
        data_matrix[contact.id, contact.id].value += contact_constraint_jacobian_contact_data(mechanism, contact, pbody)
    end
    return nothing
end

function jacobian_joint_data!(data_matrix::SparseMatrixCSC, mechanism::Mechanism{T}) where T
    # ∂body∂jointdata
    # TODO adapt this to handle cycles
    for body in mechanism.bodies
        for joint in mechanism.joints
            if (body.id == joint.parent_id) || (body.id == joint.child_id)
                data_matrix[body.id, joint.id].value += body_constraint_jacobian_joint_data(mechanism, body, joint)
            end
        end
    end
    return nothing
end

function jacobian_body_data!(data_matrix::SparseMatrixCSC, mechanism::Mechanism{T}) where T
    # ∂joint∂bodydata
    # TODO adapt this to handle cycles
    for body in mechanism.bodies
        for joint in mechanism.joints
            if (body.id == joint.parent_id) || (body.id == joint.child_id)
                data_matrix[joint.id, body.id].value += joint_constraint_jacobian_body_data(mechanism, joint, body)
            end
        end
    end
    # ∂body∂bodydata
    for pbody in mechanism.bodies
        data_matrix[pbody.id, pbody.id].value += body_constraint_jacobian_body_data(mechanism, pbody)

        for cbody in [mechanism.bodies; mechanism.origin]
            joint_links = indirect_link(pbody.id, cbody.id, mechanism.joints)
            joints = [get_joint(mechanism, id) for id in joint_links]
            for joint in joints
                ∇11, ∇12 = body_constraint_jacobian_body_data(mechanism, pbody, cbody, joint)
                (typeof(pbody) <: Body) && (data_matrix[pbody.id, pbody.id].value += ∇11)
                (typeof(pbody) <: Body && typeof(cbody) <: Body) && (data_matrix[pbody.id, cbody.id].value += ∇12)
            end
            # pretty sure this is useless because contact is never linked to two bodies
            # contact_links = indirect_link(pbody.id, cbody.id, mechanism.contacts)
            # contacts = [get_joint(mechanism, id) for id in contact_links]
            # for contact in contacts
                # data_matrix[pbody.id, cbody.id].value += body_constraint_jacobian_body_data(mechanism, pbody, cbody, contact)
            # end
        end
    end
    for contact in mechanism.contacts
        body = get_body(mechanism, contact.parent_id)
        data_matrix[body.id, body.id].value += body_constraint_jacobian_body_data(mechanism, body, contact)
    end
    # ∂contact∂bodydata
    for contact in mechanism.contacts
        pbody = get_body(mechanism, contact.parent_id)
        data_matrix[contact.id, pbody.id].value += contact_constraint_jacobian_body_data(mechanism, contact, pbody)
    end
    return nothing
end


text:
# binomial.jl
# Binomial tests
#
# Copyright (C) 2013   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export BinomialTest, SignTest

## BINOMIAL TEST

struct BinomialTest <: HypothesisTest
    p::Float64
    x::Int
    n::Int

    BinomialTest(x::Real, n::Real, p::Real=0.5) = new(p, x, n)
end

"""
    BinomialTest(x::Integer, n::Integer, p::Real = 0.5)
    BinomialTest(x::AbstractVector{Bool}, p::Real = 0.5)

Perform a binomial test of the null hypothesis that the distribution from which `x`
successes were encountered in `n` draws (or alternatively from which the vector `x` was
drawn) has success probability `p` against the alternative hypothesis that the success
probability is not equal to `p`.

Computed confidence intervals ([`confint`](@ref)) by default are Clopper-Pearson intervals.

Implements: [`pvalue`](@ref), [`confint(::BinomialTest)`](@ref)
"""
BinomialTest(x::AbstractVector{Bool}, p=0.5) =
    BinomialTest(sum(x), length(x), p)

"""
    testname(::HypothesisTest)

Returns the string value, e.g. "Binomial test" or "Sign Test".
"""
testname(::BinomialTest) = "Binomial test"
population_param_of_interest(x::BinomialTest) = ("Probability of success", x.p, x.x/x.n) # parameter of interest: name, value under h0, point estimate
default_tail(test::BinomialTest) = :both

function show_params(io::IO, x::BinomialTest, ident="")
    println(io, ident, "number of observations: $(x.n)")
    println(io, ident, "number of successes:    $(x.x)")
end

pvalue(x::BinomialTest; tail=:both) = pvalue(Binomial(x.n, x.p), x.x; tail=tail)

# Confidence interval
"""
    confint(test::BinomialTest, alpha = 0.05; tail = :both, method = :clopper_pearson)

Compute a confidence interval with coverage 1-`alpha` for a binomial proportion using one
of the following methods. Possible values for `method` are:

  - `:clopper_pearson` (default): Clopper-Pearson interval is based on the binomial
    distribution. The empirical coverage is never less than the nominal coverage of
    1-`alpha`; it is usually too conservative.
  - `:wald`: Wald (or normal approximation) interval relies on the standard approximation of
    the actual binomial distribution by a normal distribution. Coverage can be erratically
    poor for success probabilities close to zero or one.
  - `:wilson`: Wilson score interval relies on a normal approximation. In contrast to `:wald`,
    the standard deviation is not approximated by an empirical estimate, resulting in good
    empirical coverages even for small numbers of draws and extreme success probabilities.
  - `:jeffrey`: Jeffreys interval is a Bayesian credible interval obtained by using a
    non-informative Jeffreys prior. The interval is very similar to the Wilson interval.
  - `:agresti_coull`: Agresti-Coull interval is a simplified version of the Wilson interval;
    both are centered around the same value. The Agresti Coull interval has higher or equal
    coverage.
  - `:arcsine`: Confidence interval computed using the arcsine transformation to make
    ``var(p)`` independent of the probability ``p``.

# References

  * Brown, L.D., Cai, T.T., and DasGupta, A. Interval estimation for a binomial proportion.
    Statistical Science, 16(2):101–117, 2001.

# External links

  * [Binomial confidence interval on Wikipedia](https://en.wikipedia.org/wiki/
    Binomial_proportion_confidence_interval)
"""
function StatsBase.confint(x::BinomialTest, alpha::Float64=0.05; tail=:both, method=:clopper_pearson)
    check_alpha(alpha)

    if tail == :left
        (0.0, StatsBase.confint(x, alpha*2, method=method)[2])
    elseif tail == :right
        (StatsBase.confint(x, alpha*2, method=method)[1], 1.0)
    elseif tail == :both
        if method == :clopper_pearson
            ci_clopper_pearson(x, alpha)
        elseif method == :wald
            ci_wald(x, alpha)
        elseif method == :wilson
            ci_wilson(x, alpha)
        elseif method == :jeffrey
            ci_jeffrey(x, alpha)
        elseif method == :agresti_coull
            ci_agresti_coull(x, alpha)
        elseif method == :arcsine
            ci_arcsine(x, alpha)
        else
            throw(ArgumentError("method=$(method) is not implemented yet"))
        end
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

# Clopper-Pearson interval (confidence interval by inversion)
function ci_clopper_pearson(x::BinomialTest, alpha::Float64=0.05)
    (x.x == 0 ? 0.0 : quantile(Beta(x.x, x.n - x.x + 1), alpha/2),
     x.x == x.n ? 1.0 : quantile(Beta(x.x + 1, x.n - x.x), 1-alpha/2))
end

# Wald interval / normal approximation interval
function ci_wald(x::BinomialTest, alpha::Float64=0.05)
    μ = x.x / x.n
    σ = sqrt(μ*(1-μ)/x.n)
    (quantile(Normal(μ, σ), alpha/2), quantile(Normal(μ, σ), 1-alpha/2))
end

# Jeffreys interval
function ci_jeffrey(x::BinomialTest, alpha::Float64=0.05)
    (quantile(Beta(x.x + 1/2, x.n - x.x + 1/2), alpha/2), quantile(Beta(x.x + 1/2, x.n - x.x + 1/2), 1-alpha/2))
end

# Agresti-Coull interval
function ci_agresti_coull(x::BinomialTest, alpha::Float64=0.05)
    q = quantile(Normal(), 1-alpha/2)
    n = x.n + q^2
    μ = (x.x + q^2/2)/n
    σ = sqrt(μ*(1-μ)/n)
    (μ-q*σ, μ+q*σ)
end

# Wilson score interval
function ci_wilson(x::BinomialTest, alpha::Float64=0.05)
    q = quantile(Normal(), 1-alpha/2)
    p = x.x / x.n
    denominator = 1 + q^2/x.n
    μ = p + q^2/(2*x.n)
    μ /= denominator
    σ = sqrt(p*(1-p)/x.n + q^2/(4x.n^2))
    σ /= denominator
    (μ-q*σ, μ+q*σ)
end

# Arcsine transformation interval as based on Cohen's H: https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Arcsine_transformation
function ci_arcsine(x::BinomialTest, alpha::Float64=0.05)
    q = quantile(Normal(), 1-alpha/2)
    p = x.x / x.n
    z = q/(2*sqrt(x.n))
    (sin(asin(sqrt(p))-z)^2, sin(asin(sqrt(p))+z)^2)
end

## SIGN TEST

struct SignTest <: HypothesisTest
    median::Float64
    x::Int
    n::Int
    data
end

"""
    SignTest(x::AbstractVector{T<:Real}, median::Real = 0)
    SignTest(x::AbstractVector{T<:Real}, y::AbstractVector{T<:Real}, median::Real = 0)

Perform a sign test of the null hypothesis that the distribution from which `x`
(or `x - y` if `y` is provided) was drawn has median `median` against the alternative
hypothesis that the median is not equal to `median`.

Implements: [`pvalue`](@ref), [`confint`](@ref)
"""
SignTest(x::AbstractVector{T}, median::Real=0) where {T<:Real} =
    SignTest(median, sum(x .> median), sum(x .!= median), sort(x))
SignTest(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real} =
    SignTest(x - y, 0.0)

testname(::SignTest) = "Sign Test"
population_param_of_interest(x::SignTest) = ("Median", x.median, median(x.data)) # parameter of interest: name, value under h0, point estimate
default_tail(test::SignTest) = :both

function show_params(io::IO, x::SignTest, ident="")
    text1 = "number of observations:"
    text2 = "observations larger than $(x.median): "
    maxlen = length(text2)

    println(io, ident, text1, repeat(" ", maxlen-length(text1)), x.n)
    println(io, ident, text2, x.x)
end

pvalue(x::SignTest; tail=:both) = pvalue(Binomial(x.n, 0.5), x.x; tail=tail)

function StatsBase.confint(x::SignTest, alpha::Float64=0.05; tail=:both)
    check_alpha(alpha)

    if tail == :left
        q = Int(quantile(Binomial(x.n, 0.5), alpha))
        (x.data[q+1], x.median)
    elseif tail == :right
        q = Int(quantile(Binomial(x.n, 0.5), alpha))
        (x.median, x.data[end-q])
    elseif tail == :both
        q = Int(quantile(Binomial(x.n, 0.5), alpha/2))
        (x.data[q+1], x.data[end-q])
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end


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

"""
Unit tests for the dot(x, y) function
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using LinearAlgebra: norm, ⋅
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- Tests

# ------ M::Multivector

@testset "dot(M::Multivector, N::Multivector)" begin
    @test_skip 1
end

@testset "dot(M::Multivector, N::Blade)" begin
    @test_skip 1
end

@testset "dot(M::Multivector, N::Pseudoscalar)" begin
    @test_skip 1
end

@testset "dot(M::Multivector, N::Scalar)" begin
    @test_skip 1
end

@testset "dot(M::Multivector, N::One)" begin
    @test_skip 1
end

@testset "dot(M::Multivector, N::Zero)" begin
    @test_skip 1
end

@testset "dot(M::Multivector, N::Real)" begin
    @test_skip 1
end

@testset "dot(M::Multivector, N::Vector)" begin
    @test_skip 1
end

# ------ B::Blade

@testset "dot(B::Blade, M::Multivector)" begin
    @test_skip 1
end

@testset "dot(B::Blade, C::Blade)" begin
    @test_skip 1
end

@testset "dot(B::Blade, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "dot(B::Blade, C::Scalar)" begin
    @test_skip 1
end

@testset "dot(B::Blade, C::One)" begin
    @test_skip 1
end

@testset "dot(B::Blade, C::Zero)" begin
    @test_skip 1
end

@testset "dot(B::Blade, C::Real)" begin
    @test_skip 1
end

@testset "dot(B::Blade, C::Vector)" begin
    @test_skip 1
end

# ------ B::Pseudoscalar

@testset "dot(B::Pseudoscalar, M::Multivector)" begin
    @test_skip 1
end

@testset "dot(B::Pseudoscalar, C::Blade)" begin
    @test_skip 1
end

@testset "dot(B::Pseudoscalar, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "dot(B::Pseudoscalar, C::Scalar)" begin
    @test_skip 1
end

@testset "dot(B::Pseudoscalar, C::One)" begin
    @test_skip 1
end

@testset "dot(B::Pseudoscalar, C::Zero)" begin
    @test_skip 1
end

@testset "dot(B::Pseudoscalar, C::Real)" begin
    @test_skip 1
end

@testset "dot(B::Pseudoscalar, C::Vector)" begin
    @test_skip 1
end

# ------ B::Scalar

@testset "dot(B::Scalar, M::Multivector)" begin
    @test_skip 1
end

@testset "dot(B::Scalar, C::Blade)" begin
    @test_skip 1
end

@testset "dot(B::Scalar, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "dot(B::Scalar, C::Scalar)" begin
    @test_skip 1
end

@testset "dot(B::Scalar, C::One)" begin
    @test_skip 1
end

@testset "dot(B::Scalar, C::Zero)" begin
    @test_skip 1
end

@testset "dot(B::Scalar, C::Real)" begin
    @test_skip 1
end

@testset "dot(B::Scalar, C::Vector)" begin
    @test_skip 1
end

# ------ B::One

@testset "dot(B::One, M::Multivector)" begin
    @test_skip 1
end

@testset "dot(B::One, C::Blade)" begin
    @test_skip 1
end

@testset "dot(B::One, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "dot(B::One, C::Scalar)" begin
    @test_skip 1
end

@testset "dot(B::One, C::One)" begin
    @test_skip 1
end

@testset "dot(B::One, C::Zero)" begin
    @test_skip 1
end

@testset "dot(B::One, C::Real)" begin
    @test_skip 1
end

@testset "dot(B::One, C::Vector)" begin
    @test_skip 1
end

# ------ B::Zero

@testset "dot(B::Zero, M::Multivector)" begin
    @test_skip 1
end

@testset "dot(B::Zero, C::Blade)" begin
    @test_skip 1
end

@testset "dot(B::Zero, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "dot(B::Zero, C::Scalar)" begin
    @test_skip 1
end

@testset "dot(B::Zero, C::One)" begin
    @test_skip 1
end

@testset "dot(B::Zero, C::Zero)" begin
    @test_skip 1
end

@testset "dot(B::Zero, C::Real)" begin
    @test_skip 1
end

@testset "dot(B::Zero, C::Vector)" begin
    @test_skip 1
end

# ------ B::Real

@testset "dot(B::Real, M::Multivector)" begin
    @test_skip 1
end

@testset "dot(B::Real, C::Blade)" begin
    @test_skip 1
end

@testset "dot(B::Real, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "dot(B::Real, C::Scalar)" begin
    @test_skip 1
end

@testset "dot(B::Real, C::One)" begin
    @test_skip 1
end

@testset "dot(B::Real, C::Zero)" begin
    @test_skip 1
end

# ------ B::Vector

@testset "dot(B::Vector, M::Multivector)" begin
    @test_skip 1
end

@testset "dot(B::Vector, C::Blade)" begin
    @test_skip 1
end

@testset "dot(B::Vector, C::Pseudoscalar)" begin
    @test_skip 1
end

@testset "dot(B::Vector, C::Scalar)" begin
    @test_skip 1
end

@testset "dot(B::Vector, C::One)" begin
    @test_skip 1
end

@testset "dot(B::Vector, C::Zero)" begin
    @test_skip 1
end


text:
export solve!, get_eigen, get_Jacobian, get_elastic_preconditioner, get_elastic_matrix

"""
```
	function solve!(
		acft::Aircraft{Fg},
		x::AbstractVector,
		U∞::Fg;
		n_iter::Int64 = 1,
		verbose::Bool = false,
		kwargs...
	) where Fg
```

Solve aircraft aerostructural problem using Newton method

* `acft`: aircraft object
* `q`: variables of interest
* `U∞`: freestream velocity
* `n_iter`: number of iterations. Defaults to 1
* `verbose`: display residual at each iteration
* `kwargs`: keyword arguments for `state_space`
"""
function solve!(
	acft::Aircraft{Fg},
	x::AbstractVector,
	U∞::Fg;
	n_iter::Int64 = 1,
	verbose::Bool = false,
	kwargs...
) where Fg

	f = x -> begin
		qd, _ = state_space(
			acft,
			x,
			U∞;
			kwargs...
		)

		qd
	end

	for nit = 1:n_iter
		r = f(x)

		if verbose
			println("Iteration $nit - residual norm $(norm(r))")
		end

		J = ForwardDiff.jacobian(
			f,
			x
		)

		x .-= J \ r
	end

end

"""
Get eigenvalues and eigenvectors of the state-space description
of the aircraft

* `acft`: aircraft object
* `q`: variables of interest
* `U∞`: freestream velocity
* `kwargs`: keyword arguments for `state_space`
"""
function get_eigen(
	acft::Aircraft{Fg},
	x::AbstractVector,
	U∞::Fg;
	kwargs...
) where Fg

	f = x -> begin
		qd, _ = state_space(
			acft,
			x,
			U∞;
			kwargs...
		)

		qd
	end

	J = ForwardDiff.jacobian(
		f,
		x
	)

	eigen(J)

end

"""
Get Jacobian of the state-space description
of the aircraft

* `acft`: aircraft object
* `q`: variables of interest
* `U∞`: freestream velocity
* `kwargs`: keyword arguments for `state_space`
"""
function get_Jacobian(
	acft::Aircraft{Fg},
	x::AbstractVector,
	U∞::Fg;
	kwargs...
) where Fg

	f = x -> begin
		qd, _ = state_space(
			acft,
			x,
			U∞;
			kwargs...
		)

		qd
	end

	J = q -> ForwardDiff.jacobian(
		f,
		x .+ q
	)

	J

end

"""
	```
	function get_elastic_matrix(
		acft::Aircraft{Fg};
		fixed_points::Vector{Int64} = Int64[],
		amplification_factor::Real = Fmax
	) where Fg
	```

Get matrix `inv(M) * K` in sparse format
"""
function get_elastic_matrix(
	acft::Aircraft{Fg};
	fixed_points::Vector{Int64} = Int64[],
	amplification_factor::Real = Fmax
) where Fg

	q = get_state(acft, 1.0)

	f = (y, x) -> begin
		qd, _ = state_space(
			acft,
			x,
			1.0;
			aerodynamic_forces = false,
			structural_damping = false,
			fixed_points = fixed_points,
			amplification_factor = amplification_factor
		)

		y .= qd

		y
	end

	HD = ndofs(acft) ÷ 2

	n_elements = 0

	for b in acft.beams
		for j = (6 * (b.ipt1 - 1) + 1):(6 * b.ipt1)
			for i = (6 * (b.ipt1 - 1) + 1 + HD):(6 * b.ipt1 + HD)
				n_elements += 1
			end

			for i = (6 * (b.ipt2 - 1) + 1 + HD):(6 * b.ipt2 + HD)
				n_elements += 1
			end
		end

		for j = (6 * (b.ipt2 - 1) + 1):(6 * b.ipt2)
			for i = (6 * (b.ipt1 - 1) + 1 + HD):(6 * b.ipt1 + HD)
				n_elements += 1
			end

			for i = (6 * (b.ipt2 - 1) + 1 + HD):(6 * b.ipt2 + HD)
				n_elements += 1
			end
		end
	end

	for i = 1:HD
		n_elements += 1
	end

	for f in fixed_points
		n_elements += 12
	end

	for lnk in acft.links
		for i = (6 * (lnk.i1 - 1) + 1 + HD):(6 * lnk.i1 + HD)
			for j = (6 * (lnk.i1 - 1) + 1):(6 * lnk.i1)
				n_elements += 1
			end

			if lnk.i2 != 0
				for j = (6 * (lnk.i2 - 1) + 1):(6 * lnk.i2)
					n_elements += 1
				end
			end

			for j = (6 * (lnk.i1 - 1) + 1 + HD):(6 * lnk.i1 + HD)
				n_elements += 1
			end

			if lnk.i2 != 0
				for j = (6 * (lnk.i2 - 1) + 1 + HD):(6 * lnk.i2 + HD)
					n_elements += 1
				end
			end
		end

		if lnk.i2 != 0
			for i = (6 * (lnk.i2 - 1) + 1 + HD):(6 * lnk.i2 + HD)
				for j = (6 * (lnk.i2 - 1) + 1):(6 * lnk.i2)
					n_elements += 1
				end

				for j = (6 * (lnk.i2 - 1) + 1 + HD):(6 * lnk.i2 + HD)
					n_elements += 1
				end

				for j = (6 * (lnk.i1 - 1) + 1):(6 * lnk.i1)
					n_elements += 1
				end

				for j = (6 * (lnk.i1 - 1) + 1 + HD):(6 * lnk.i1 + HD)
					n_elements += 1
				end
			end
		end
	end

	data = ones(Fg, n_elements)
	n_rows = Vector{Int64}(undef, n_elements)
	n_cols = Vector{Int64}(undef, n_elements)

	n_elements = 0

	for b in acft.beams
		for j = (6 * (b.ipt1 - 1) + 1):(6 * b.ipt1)
			for i = (6 * (b.ipt1 - 1) + 1 + HD):(6 * b.ipt1 + HD)
				n_elements += 1

				n_rows[n_elements] = i
				n_cols[n_elements] = j
			end

			for i = (6 * (b.ipt2 - 1) + 1 + HD):(6 * b.ipt2 + HD)
				n_elements += 1

				n_rows[n_elements] = i
				n_cols[n_elements] = j
			end
		end

		for j = (6 * (b.ipt2 - 1) + 1):(6 * b.ipt2)
			for i = (6 * (b.ipt1 - 1) + 1 + HD):(6 * b.ipt1 + HD)
				n_elements += 1

				n_rows[n_elements] = i
				n_cols[n_elements] = j
			end

			for i = (6 * (b.ipt2 - 1) + 1 + HD):(6 * b.ipt2 + HD)
				n_elements += 1

				n_rows[n_elements] = i
				n_cols[n_elements] = j
			end
		end
	end

	for i = 1:HD
		n_elements += 1

		n_rows[n_elements] = i
		n_cols[n_elements] = i + HD
	end

	for f in fixed_points
		for i = (6 * (f - 1) + 1):(6 * f)
			n_elements += 1

			n_rows[n_elements] = i
			n_cols[n_elements] = i

			n_elements += 1

			n_rows[n_elements] = i + HD
			n_cols[n_elements] = i + HD
		end
	end

	for lnk in acft.links
		for i = (6 * (lnk.i1 - 1) + 1 + HD):(6 * lnk.i1 + HD)
			for j = (6 * (lnk.i1 - 1) + 1):(6 * lnk.i1)
				n_elements += 1

				n_rows[n_elements] = i
				n_cols[n_elements] = j
			end

			if lnk.i2 != 0
				for j = (6 * (lnk.i2 - 1) + 1):(6 * lnk.i2)
					n_elements += 1

					n_rows[n_elements] = i
					n_cols[n_elements] = j
				end
			end

			for j = (6 * (lnk.i1 - 1) + 1 + HD):(6 * lnk.i1 + HD)
				n_elements += 1

				n_rows[n_elements] = i
				n_cols[n_elements] = j
			end

			if lnk.i2 != 0
				for j = (6 * (lnk.i2 - 1) + 1 + HD):(6 * lnk.i2 + HD)
					n_elements += 1

					n_rows[n_elements] = i
					n_cols[n_elements] = j
				end
			end
		end

		if lnk.i2 != 0
			for i = (6 * (lnk.i2 - 1) + 1 + HD):(6 * lnk.i2 + HD)
				for j = (6 * (lnk.i2 - 1) + 1):(6 * lnk.i2)
					n_elements += 1

					n_rows[n_elements] = i
					n_cols[n_elements] = j
				end

				for j = (6 * (lnk.i2 - 1) + 1 + HD):(6 * lnk.i2 + HD)
					n_elements += 1

					n_rows[n_elements] = i
					n_cols[n_elements] = j
				end

				for j = (6 * (lnk.i1 - 1) + 1):(6 * lnk.i1)
					n_elements += 1

					n_rows[n_elements] = i
					n_cols[n_elements] = j
				end

				for j = (6 * (lnk.i1 - 1) + 1 + HD):(6 * lnk.i1 + HD)
					n_elements += 1

					n_rows[n_elements] = i
					n_cols[n_elements] = j
				end
			end
		end
	end

	jac = sparse(n_rows, n_cols, data)

	forwarddiff_color_jacobian!(
		jac,
		f,
		q;
		colorvec = matrix_colors(jac)
	)

	jac

end

"""
	```
	function get_elastic_preconditioner(
		acft::Aircraft{Fg};
		fixed_points::Vector{Int64} = Int64[]
	) where Fg
	```

Get elastic matrix `inv(K) * M` preconditioner (LU decomposition
of Jacobian when disregarded elastic forces)
"""
function get_elastic_preconditioner(
	acft::Aircraft{Fg};
	fixed_points::Vector{Int64} = Int64[]
) where Fg

	jac = get_elastic_matrix(acft; fixed_points = fixed_points)

	lu(jac)

end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
module UnregisteredWithProject

using Example

export f
f(x) = x

end


text:
# Adapted from:
# SOSDEMO10 --- Set containment
# Section 3.10 of SOSTOOLS User's Manual

function sosdemo10_test(optimizer, config::MOIT.TestConfig)
    @polyvar x[1:2]

    eps = 1e-6
    p = x[1]^2+x[2]^2
    gamma = 1
    g0 = 2*x[1]
    theta = 1

    model = _model(optimizer)
    PolyJuMP.setpolymodule!(model, SumOfSquares)

    # FIXME s should be sos ?
    # in SOSTools doc it is said to be SOS
    # but in the demo it is not constrained so
    Z = monomials(x, 0:4)
    @variable(model, s, Poly(Z))

    Z = monomials(x, 2:3)
    @variable(model, g1, Poly(Z))

    Sc = [theta^2-s*(gamma-p) g0+g1; g0+g1 1]

    @SDconstraint(model, Matrix(eps * I, 2, 2) ⪯ Sc)

    JuMP.optimize!(model)

    # Program is feasible, that is, the set
    # { x |((g0+g1) + theta)(theta - (g0+g1)) >=0 }
    # contains the set
    # { x | p <= gamma }

    # ALMOST_OPTIMAL and NEARLY_FEASIBLE_POINT for CSDP
    @test JuMP.termination_status(model) in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
    @test JuMP.primal_status(model) in [MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT]

    #@show JuMP.value(s)
    #@show JuMP.value(g1)
end
sd_tests["sosdemo10"] = sosdemo10_test


text:
module FstFileFormat

using RCall
# using IterableTables # to enable conversion between table types
import DataFrames.DataFrame

export read, write, readmeta, install_fst, fst_installed

const FST_NOT_INSTALLED_ERR_MSG = "fst package not installed\n run 'FstFileFormat.install_fst()' or install fst manually in R using 'install.packages('fst')'"

"""
Install fst package if not already installed
"""
function install_fst()
  R"""
    if(!require(fst)) {
      install.packages("fst")
    }
  """
end

"""
Returns `true` if fst is installed and `false` otherwise
"""
function fst_installed()
  R"""
    memory.limit(1e10)
    fst_installed <- require(fst)
  """
  @rget fst_installed
  fst_installed
end

"""
Reads a fst file and return a DataFrame

**Arguments**
* `path` : path to the fst file
* `columns` : an array of strings return only these columns in the DataFrame
* `from`: read from which row; defaults from `1`
* `to` : read to which row; defaults to `[]` which reads all rows
"""
function read(path; columns = [], from = 1, to = [])
  if !fst_installed()
    throw(ErrorException(FST_NOT_INSTALLED_ERR_MSG))
  end

  @rput columns
  @rput to
  R"""
    if(length(columns) == 0) {
      columns = NULL
    }
    if(length(to) == 0) {
      to = NULL
    }
    dt = fst::read_fst($path, columns = columns, from = $from, to = to)
  """
  @rget dt
  return dt
end

"""
Write a fst file and return

**Arguments**
* `x` : a DataFrame or an object that can be turned into a DataFrame via DataFrame(x)
* `path` : path to the fst file
* `compress` : an integer from 0 to 100 indicating the compression level; higher means more compressed
"""
function write(x, path)
  write(x, path, 30)
end

function write(x, path, compress)
  if !fst_installed()
    throw(ErrorException(FST_NOT_INSTALLED_ERR_MSG))
  end

  xdf = DataFrame(x)
  @rput xdf
  R"""
    dt <- fst::write_fst(xdf, $path, compress = $compress)
  """
  @rget dt
  return dt
end

"""
Reads a fst file's metadata and returns a Dict of the attributes

**Arguments**
* `path` : path to the fst file
"""
function readmeta(path)
  if !fst_installed()
    throw(ErrorException(FST_NOT_INSTALLED_ERR_MSG))
  end

  R"""
    meta <- fst::fst.metadata($path)
  """
  @rget meta
  return meta
end

end # module


text:
# Lebesgue measure

export Lebesgue
struct Lebesgue{X} <: AbstractMeasure{X} end

Lebesgue(X) = Lebesgue{X}()

baseMeasure(μ::Lebesgue{X}) where {X} = μ


text:
# Test disjoint set

using DataStructures
using Base.Test

s = IntDisjointSets(10)

@test length(s) == 10
@test num_groups(s) == 10

for i = 1:10
    @test find_root(s, i) == i
end
@test_throws BoundsError find_root(s, 11)

@test !in_same_set(s, 2, 3)

union!(s, 2, 3)
@test num_groups(s) == 9
@test in_same_set(s, 2, 3)
@test find_root(s, 3) == 2
union!(s, 3, 2)
@test num_groups(s) == 9
@test in_same_set(s, 2, 3)
@test find_root(s, 3) == 2

# We cannot support arbitrary indexing and still use @inbounds with IntDisjointSets
# (and it's not useful anyway)
@test_throws MethodError push!(s, 22)

@test push!(s) == 11
@test num_groups(s) == 10

union!(s, 8, 7)
union!(s, 5, 6)
union!(s, 8, 5)
@test num_groups(s) == 7
@test find_root(s, 6) == 8
union!(s, 2, 6)
@test find_root(s, 2) == 8

s = DisjointSets{Int}(1:10)

@test length(s) == 10
@test num_groups(s) == 10

r = [find_root(s, i) for i in 1 : 10]
@test isa(r, Vector{Int})
@test isequal(r, [1:10])

for i = 1 : 5
    x = 2 * i - 1
    y = 2 * i
    union!(s, x, y)
end

@test length(s) == 10
@test num_groups(s) == 5

r0 = [1, 1, 3, 3, 5, 5, 7, 7, 9, 9]
r = [find_root(s, i) for i in 1 : 10]
@test isa(r, Vector{Int})
@test isequal(r, r0)

union!(s, 1, 4)
union!(s, 3, 5)
union!(s, 7, 9)

@test length(s) == 10
@test num_groups(s) == 2

r0 = [1, 1, 1, 1, 1, 1, 7, 7, 7, 7]
r = [find_root(s, i) for i in 1 : 10]
@test isa(r, Vector{Int})
@test isequal(r, r0)

push!(s, 17)

@test length(s) == 11
@test num_groups(s) == 3

r0 = [1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 11]
r = [find_root(s, i) for i in [1 : 10, 17] ]
@test isa(r, Vector{Int})
@test isequal(r, r0)


text:
immutable DifferentiableFunction{F <: Function, G <: Function} <: MathProgBase.AbstractNLPEvaluator
    func::F
    gradient!::G
end

DifferentiableFunction(f::Function) = DifferentiableFunction(f, (result, x) -> ForwardDiff.gradient!(result, f, x))

function MathProgBase.initialize(gd::DifferentiableFunction, requested_features::Vector{Symbol})
end

MathProgBase.features_available(gd::DifferentiableFunction) = [:Grad]

MathProgBase.jac_structure(gd::DifferentiableFunction) = (tuple(), tuple())

MathProgBase.eval_grad_f(gd::DifferentiableFunction, g, x) = gd.gradient!(g, x)

MathProgBase.eval_f(gd::DifferentiableFunction, x) = gd.func(x)

function getmodel(problem::DifferentiableFunction,
                num_vars::Integer,
                solver=MathProgBase.defaultNLPsolver,
                lower_bound=[-Inf for i in 1:num_vars],
                upper_bound=[Inf for i in 1:num_vars])
    model = MathProgBase.NonlinearModel(solver)
    MathProgBase.loadproblem!(model, num_vars, 0, lower_bound, upper_bound, Float64[], Float64[], :Min, problem)
    model
end

getmodel(f::Function, args...) = getmodel(DifferentiableFunction(f), args...)