text:
"""
    max_min_driving_force(
        model::MetabolicModel,
        reaction_standard_gibbs_free_energies::Dict{String,Float64},
        optimizer;
        flux_solution::Dict{String,Float64} = Dict{String,Float64}(),
        proton_ids::Vector{String} = ["h_c", "h_e"],
        water_ids::Vector{String} = ["h2o_c", "h2o_e"],
        constant_concentrations::Dict{String,Float64} = Dict{String,Float64}(),
        concentration_ratios::Dict{Tuple{String,String},Float64} = Dict{
            Tuple{String,String},
            Float64,
        }(),
        concentration_lb = 1e-9,
        concentration_ub = 100e-3,
        T = _constants.T,
        R = _constants.R,
        small_flux_tol = 1e-6,
        modifications = [],
        ignore_reaction_ids = [],
    )

Perform a max-min driving force analysis on the `model`, as defined by Noor, et al.,
"Pathway thermodynamics highlights kinetic obstacles in central metabolism.", PLoS
computational biology, 2014.

The function uses the supplied `optimizer` and `reaction_standard_gibbs_free_energies`.
Optionally, `flux_solution` can be used to set the directions of each reaction in `model`
(all reactions are assumed to proceed forward and are active by default). The supplied
`flux_solution` should be free of internal cycles i.e. thermodynamically consistent. This
optional input is important if a reaction in `model` normally runs in reverse (negative
flux). Note, reactions in `flux_solution` that are smaller than `small_flux_tol` are also
ignored in the analysis function (for numerical stability).

The max-min driving force algorithm returns the Gibbs free energy of the reactions, the
concentrations of metabolites and the actual maximum minimum driving force. The optimization
problem solved is:
```
max min -ΔᵣG
s.t. ΔᵣG = ΔrG⁰ + R T S' ln(C)
     ΔᵣG ≤ 0
     ln(Cₗ) ≤ ln(C) ≤ ln(Cᵤ)
```
where `ΔrG` are the Gibbs energies dissipated by the reactions, R is the gas constant, T is
the temperature, S is the stoichiometry of the model, and C is the vector of metabolite
concentrations (and their respective lower and upper bounds).

In case no feasible solution exists, `nothing` is returned.

Reactions specified in `ignore_reaction_ids` are internally ignored when calculating the
max-min driving force. This should include water and proton importers.

Since biochemical thermodynamics are assumed, the `proton_ids` and `water_ids` need to be
specified so that they can be ignored in the calculations. Effectively this assumes an
aqueous environment at constant pH is used.

`constant_concentrations` is used to fix the concentrations of certain metabolites (such as
CO₂). `concentration_ratios` is used to specify additional constraints on metabolite pair
concentrations (typically, this is done with various cofactors such as the ATP/ADP ratio.
For example, you can fix the concentration of ATP to be always 5× higher than of ADP by
specifying `Dict(("ATP","ADP") => 5.0)`

`concentration_lb` and `concentration_ub` set the `Cₗ` and `Cᵤ` in the
optimization problems.

`T` and `R` can be specified in the corresponding units; defaults are K and kJ/K/mol.
"""
function max_min_driving_force(
    model::MetabolicModel,
    reaction_standard_gibbs_free_energies::Dict{String,Float64},
    optimizer;
    flux_solution::Dict{String,Float64} = Dict{String,Float64}(),
    proton_ids::Vector{String} = ["h_c", "h_e"],
    water_ids::Vector{String} = ["h2o_c", "h2o_e"],
    constant_concentrations::Dict{String,Float64} = Dict{String,Float64}(),
    concentration_ratios::Dict{Tuple{String,String},Float64} = Dict{
        Tuple{String,String},
        Float64,
    }(),
    concentration_lb = 1e-9,
    concentration_ub = 100e-3,
    T = _constants.T,
    R = _constants.R,
    small_flux_tol = 1e-6,
    modifications = [],
    ignore_reaction_ids = [],
)
    opt_model = Model(optimizer)

    @variables opt_model begin
        mmdf
        logcs[1:n_metabolites(model)]
        dgrs[1:n_reactions(model)]
    end

    # set proton log concentration to zero so that it won't impact any calculations (biothermodynamics assumption)
    proton_idxs = Int.(indexin(proton_ids, metabolites(model)))
    for idx in proton_idxs
        JuMP.fix(logcs[idx], 0.0)
    end

    # set water concentrations to zero (aqueous condition assumptions)
    water_idxs = Int.(indexin(water_ids, metabolites(model)))
    for idx in water_idxs
        JuMP.fix(logcs[idx], 0.0)
    end

    # only consider reactions with supplied thermodynamic data AND have a flux bigger than
    # small_flux_tol => finds a thermodynamic profile that explains flux_solution
    active_rids = filter(
        rid ->
            haskey(reaction_standard_gibbs_free_energies, rid) &&
                abs(get(flux_solution, rid, small_flux_tol / 2)) > small_flux_tol &&
                !(rid in ignore_reaction_ids),
        reactions(model),
    )
    active_ridxs = Int.(indexin(active_rids, reactions(model)))

    # give dummy dG0 for reactions that don't have data
    dg0s =
        [get(reaction_standard_gibbs_free_energies, rid, 0.0) for rid in reactions(model)]

    S = stoichiometry(model)

    @constraint(opt_model, dgrs .== dg0s .+ (R * T) * S' * logcs)

    # thermodynamics should correspond to the fluxes
    flux_signs = [sign(get(flux_solution, rid, 1.0)) for rid in reactions(model)]

    # only constrain reactions that have thermo data
    @constraint(opt_model, dgrs[active_ridxs] .* flux_signs[active_ridxs] .<= 0)

    # add the absolute bounds
    missing_mets =
        [mid for mid in keys(constant_concentrations) if !(mid in metabolites(model))]
    !isempty(missing_mets) &&
        throw(DomainError(missing_mets, "metabolite(s) not found in model."))
    for (midx, mid) in enumerate(metabolites(model)) # idx in opt_model (missing ignore_metabolites)
        midx in water_idxs && continue
        midx in proton_idxs && continue
        if haskey(constant_concentrations, mid)
            JuMP.fix(logcs[midx], log(constant_concentrations[mid]))
        else
            # this metabolite needs bounds
            @constraint(
                opt_model,
                log(concentration_lb) <= logcs[midx] <= log(concentration_ub)
            )
        end
    end

    # add the relative bounds
    for ((mid1, mid2), val) in concentration_ratios
        idxs = indexin([mid1, mid2], metabolites(model)) # TODO: this is not performant
        any(isnothing.(idxs)) &&
            throw(DomainError((mid1, mid2), "metabolite pair not found in model."))
        @constraint(opt_model, logcs[idxs[1]] == log(val) + logcs[idxs[2]])
    end

    @constraint(opt_model, mmdf .<= -dgrs[active_ridxs] .* flux_signs[active_ridxs])

    @objective(opt_model, Max, mmdf)

    # apply the modifications, if any
    for mod in modifications
        mod(model, opt_model)
    end

    optimize!(opt_model)

    is_solved(opt_model) || return nothing

    return (
        mmdf = value(opt_model[:mmdf]),
        dg_reactions = Dict(
            rid => value(opt_model[:dgrs][i]) for (i, rid) in enumerate(reactions(model))
        ),
        concentrations = Dict(
            mid => exp(value(opt_model[:logcs][i])) for
            (i, mid) in enumerate(metabolites(model))
        ),
    )
end

"""
    max_min_driving_force_variability(
        model::MetabolicModel,
        reaction_standard_gibbs_free_energies::Dict{String,Float64},
        optimizer;
        workers =[myid()],
        optimal_objective_value = nothing,
        bounds = z -> (z, Inf),
        flux_solution::Dict{String,Float64} = Dict{String,Float64}(),
        proton_ids::Vector{String} = ["h_c", "h_e"],
        water_ids::Vector{String} = ["h2o_c", "h2o_e"],
        constant_concentrations::Dict{String,Float64} = Dict{String,Float64}(),
        concentration_ratios::Dict{Tuple{String,String},Float64} = Dict{
            Tuple{String,String},
            Float64,
        }(),
        concentration_lb = 1e-9,
        concentration_ub = 100e-3,
        T = _constants.T,
        R = _constants.R,
        small_flux_tol = 1e-6,
        modifications = [],
        ignore_reaction_ids = [],
    )

Perform a variant of flux variability analysis on a max min driving force type problem.
Arguments are forwarded to [`max_min_driving_force`](@ref). Calls [`screen`](@ref)
internally and possibly distributes computation across `workers`. If
`optimal_objective_value = nothing`, the function first performs regular max min driving
force analysis to find the max min driving force of the model and sets this to
`optimal_objective_value`. Then iteratively maximizes and minimizes the driving force across
each reaction, and then the concentrations while staying close to the original max min
driving force as specified in `bounds`.

The `bounds` is a user-supplied function that specifies the max min driving force bounds for
the variability optimizations, by default it restricts the flux objective value to the
precise optimum reached in the normal max min driving force analysis. It can return `-Inf`
and `Inf` in first and second pair to remove the limit. Use [`gamma_bounds`](@ref) and
[`objective_bounds`](@ref) for simple bounds.

Returns a matrix of solutions to [`max_min_driving_force`](@ref) additionally constrained as
described above, where the rows are in the order of the reactions and then the metabolites
of the `model`. For the reaction rows the first column is the maximum dG of that reaction,
and the second column is the minimum dG of that reaction subject to the above constraints.
For the metabolite rows, the first column is the maximum concentration, and the second column
is the minimum concentration subject to the constraints above.
"""
function max_min_driving_force_variability(
    model::MetabolicModel,
    reaction_standard_gibbs_free_energies::Dict{String,Float64},
    optimizer;
    workers = [myid()],
    optimal_objective_value = nothing,
    bounds = z -> (z, Inf),
    modifications = [],
    kwargs...,
)
    if isnothing(optimal_objective_value)
        initsol = max_min_driving_force(
            model,
            reaction_standard_gibbs_free_energies,
            optimizer;
            modifications,
            kwargs...,
        )
        mmdf = initsol.mmdf
    else
        mmdf = optimal_objective_value
    end

    lb, ub = bounds(mmdf)

    dgr_variants = [
        [[_mmdf_add_df_bound(lb, ub), _mmdf_dgr_objective(ridx, sense)]] for
        ridx = 1:n_reactions(model), sense in [MOI.MAX_SENSE, MOI.MIN_SENSE]
    ]
    concen_variants = [
        [[_mmdf_add_df_bound(lb, ub), _mmdf_concen_objective(midx, sense)]] for
        midx = 1:n_metabolites(model), sense in [MOI.MAX_SENSE, MOI.MIN_SENSE]
    ]

    return screen(
        model;
        args = [dgr_variants; concen_variants],
        analysis = (m, args) -> max_min_driving_force(
            m,
            reaction_standard_gibbs_free_energies,
            optimizer;
            modifications = [args; modifications],
            kwargs...,
        ),
        workers,
    )
end

"""
Helper function to change the objective to optimizing some dG.
"""
function _mmdf_dgr_objective(ridx, sense)
    (model, opt_model) -> begin
        @objective(opt_model, sense, opt_model[:dgrs][ridx])
    end
end

"""
Helper function to change the objective to optimizing some concentration.
"""
function _mmdf_concen_objective(midx, sense)
    (model, opt_model) -> begin
        @objective(opt_model, sense, opt_model[:logcs][midx])
    end
end

"""
Helper function to add a new constraint on the driving force.
"""
function _mmdf_add_df_bound(lb, ub)
    (model, opt_model) -> begin
        if lb == ub
            fix(opt_model[:mmdf], lb; force = true)
        else
            @constraint(opt_model, lb <= opt_model[:mmdf] <= ub)
        end
    end
end


text:
function setup_jld2dataset_test()
    datasets = [
        ("d1", rand(2, 10)),
        ("g1/d1", rand(10)),
        ("g1/d2", string.('a':'j')),
        ("g2/g1/d1", rand(Bool, 3, 3, 10))
    ]

    fid = jldopen("test.jld2", "w")
    for (path, data) in datasets
        fid[path] = data
    end
    close(fid)

    return first.(datasets), Tuple(last.(datasets))
end
cleanup_jld2dataset_test() = rm("test.jld2")

@testset "JLD2Dataset" begin
    paths, datas = setup_jld2dataset_test()
    @testset "Single path" begin
        dataset = JLD2Dataset("test.jld2", "d1")
        for i in 1:10
            @test @inferred(getobs(dataset, i)) == getobs(datas[1], i)
        end
        @test numobs(dataset) == 10
        close(dataset)
    end
    @testset "Multiple paths" begin
        dataset = JLD2Dataset("test.jld2", paths)
        for i in 1:10
            @test @inferred(getobs(dataset, i)) == getobs(datas, i)
        end
        @test numobs(dataset) == 10
        close(dataset)
    end
    cleanup_jld2dataset_test()
end

ext: jl


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Testing

@testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin

    # field problem
    block = Problem(Elasticity, "BLOCK", 2)
    block.properties.formulation = :plane_stress
    block.properties.finite_strain = true
    block.properties.geometric_stiffness = true

    nodes = Dict{Int, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])

    element = Element(Quad4, [1, 2, 3, 4])
    update!(element, "geometry", nodes)
    update!(element, "youngs modulus", 288.0)
    update!(element, "poissons ratio", 1/3)
    push!(block, element)

    # boundary conditions
    bc = Problem(Dirichlet, "bc", 2, "displacement")
    bel1 = Element(Seg2, [1, 2])
    bel2 = Element(Seg2, [3, 4])
    bel3 = Element(Seg2, [4, 1])
    update!([bel1, bel2, bel3], "geometry", nodes)
    update!(bel1, "displacement 2", 0.0)
    update!(bel2, "displacement 2", 0.5)
    update!(bel3, "displacement 1", 0.0)
    push!(bc, bel1, bel2, bel3)

    solver = NonlinearSolver("solve block problem")
    push!(solver, block, bc)
    solver()

    # from code aster
    eps_expected = [-2.08333312468287E-01, 6.25000000000000E-01, 0.0]
    sig_expected = [ 4.50685020821470E-06, 4.62857140373777E+02, 0.0]
    u3_expected =  [-2.36237356855269E-01, 5.00000000000000E-01]

    u3 = reshape(block.assembly.u, 2, 4)[:, 3]
    info("u3 = $u3")
    @test isapprox(u3, u3_expected, atol=1.0e-5)

#= TODO: to postprocess
    info("strain")
    for ip in get_integration_points(element)
        eps = ip("strain")
        @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] eps[1] eps[2] eps[3]
        @test isapprox(eps, eps_expected, atol=1.0e-5)
    end
    info("cauchy stress")
    for ip in get_integration_points(element)
        sig = ip("cauchy stress")
        @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3]
        @test isapprox(sig, sig_expected)
    end
    info("pk2 stress")
    for ip in get_integration_points(element)
        sig = ip("pk2 stress")
        @printf "%i | %8.3f %8.3f | %8.3f %8.3f %8.3f\n" ip.id ip.coords[1] ip.coords[2] sig[1] sig[2] sig[3]
        @test isapprox(sig, sig_expected)
    end
=#
end


text:
using Base.Sort
using LinearAlgebra

@doc raw"""
This is a traditional cubic spline. It has N+1 points on the abcissa, ``\tau``.
It has N cubic interpolants in `c`, which is size 4 x N.
At each interval, the polynomial is in Horner Form,

``P_i(x) = c_{1,i} + c_{2,i}(x-x_i) + c_{3,i}(x-x_i)^2 + c_{4,i}(x-x_i)^3``

This spline has one type for the abcissa and one for the ordinate coefficients.
It should be the case that one(T) * one(X) is of type T.

This can be a spline of another order. It's just a matter of the dimensions
of the coefficient array, `c`.
"""
struct CubicSpline{X,T}
    τ::AbstractVector{X}
    c::AbstractMatrix{T}
end

"""
((cs.c[4, i] * Δ + cs.c[3, i]) * Δ + cs.c[2, i]) * Δ + cs.c[1, i]
"""
function horner_in_interval(cs::CubicSpline{X,T}, i, x) where {X,T}
    Δ = x - cs.τ[i]
    m = size(cs.c, 1)
    if m > 0
        total = cs.c[m, i]
        for d = (m - 1):-1:1
            total *= Δ
            total += cs.c[d, i]
        end
        total
    else
        zero(T)
    end
end

function (cs::CubicSpline)(x::A) where {A <: Real}
    # Index of first greater-than-or-equal-to x.
    i = Sort.searchsortedfirst(cs.τ, x) - 1
    # For points off the sides, use the closest polynomial.
    i = max(min(length(cs.τ) - 1, i), 1)
    horner_in_interval(cs, i, x)
end

function evaluate!(cs::CubicSpline, x::AbstractVector, y::AbstractVector)
    # Index of first greater-than-or-equal-to x.
    k = 1
    i = Sort.searchsortedfirst(cs.τ, x[k]) - 1
    i = max(min(length(cs.τ) - 1, i), 1)
    while k <= length(x)
        y[k] = horner_in_interval(cs, i, x[k])
        k += 1
        while i + 1 < length(cs.τ) && cs.τ[i + 1] < x[k]
            i += 1
        end
    end
end

"""
Integrate a spline from a knot at `i` to the value `x`.
# (((A(1/4)*cs.c[4, i] * Δ + A(1/3)*cs.c[3, i]) * Δ + A(1/2)*cs.c[2, i]) * Δ + cs.c[1, i])*Δ
"""
function integral_in_interval(cs::CubicSpline{X,T}, i, x::A) where {A <: Real, X, T}
    Δ = x - cs.τ[i]
    m = size(cs.c, 1)
    if m > 0
        total = zero(T)
        for d = m:-1:1
            total += cs.c[d, i] / T(d)
            total *= Δ
        end
        total
    else
        zero(T)
    end
end

"""
Integrate a spline from `x1` to `x2`.
"""
function integrate(cs::CubicSpline, x1::A, x2::A) where {A <: Real}
    # Index of first greater-than-or-equal-to x.
    i = Sort.searchsortedfirst(cs.τ, x1) - 1
    # For points off the sides, use the closest polynomial.
    i = max(min(length(cs.τ) - 1, i), 1)
    j = Sort.searchsortedfirst(cs.τ, x2) - 1
    # For points off the sides, use the closest polynomial.
    j = max(min(length(cs.τ) - 1, j), 1)
    intervening_intervals = zero(A)
    for k in i:(j-1)
        intervening_intervals += integral_in_interval(cs, k, cs.τ[k + 1])
    end
    intervening_intervals + integral_in_interval(cs, j, x2) - integral_in_interval(cs, i, x1)
end

"""
Given a spline, create a new spline that is its derivative.
"""
function derivative(cs::CubicSpline{X,T}) where {X,T}
    m = size(cs.c, 1)
    c2 = similar(cs.c[2:m, :])
    for i = 1:size(cs.c, 2)
        for d = 1:(m - 1)
            c2[d, i] = T(d) * cs.c[d + 1, i]
        end
    end
    CubicSpline{X,T}(cs.τ, c2)
end

"""
Multiply two splines of order `m1` and `m2` to get a spline of order `m1+m2`.
This only works if the two splines have the same abcissa.
"""
function Base.:*(cs1::CubicSpline{X,T}, cs2::CubicSpline{X,T}) where {X,T}
    m1 = size(cs1.c, 1)
    m2 = size(cs2.c, 1)
    c3 = zeros(T, m1 + m2 - 1, size(cs1.c, 2))
    for i in 1:size(cs1.c, 2)
        for d1 = 1:m1
            for d2 = 1:m2
                c3[d1 + d2 - 1, i] += cs1.c[d1, i] * cs2.c[d2, i]
            end
        end
    end
    CubicSpline{X, T}(cs1.τ, c3)
end

@doc raw"""
τ are the abcissa from 1 to N+1. f are the known values of length N+1.
fp_bounds are derivatives at the first and last endpoint.
fp is the output, of length N + 1.
This is from Conte and deBoor Sec. 4.8, Eq. 4.60. The equation is this, with ``i``
running from 2 to N.

``(\Delta x_i) s_{i-1}+2(\Delta x_{i-1}+\Delta x_i)s_i +(\Delta x_{i-1}s_{i+1}) =
3(f[x_{i-1},x_i]\Delta x_i+f[x_i,x_{i+1}]\Delta x_{i-1})``

Here, ``\Delta x_i = x_{i+1} -x{i}`` and

``f[x_i,x_{i+1}] = \frac{f[x_{i+1}] - f[x_i]}{x_{i+1}-x_i}``

"""
function global_derivatives!(τ, f::AbstractVector{T}, fp) where {T <: Real}
    # This implementation builds the tridiagonal matrix and then solves it.
    # You could write a Gaussian elimination to generate terms on the fly, saving memory.
    # There are N-1 equations.
    N = length(τ) - 1
    # lower, mid, and upper of tridiagonal matrix.
    dl = zeros(T, N - 2)  # indices from 3 to N
    d = zeros(T, N - 1) # indices from 2 to N
    du = zeros(T, N - 2) # indices from 2 to N-1
    rhs = zeros(T, N - 1) # indices from 2 to N

    for il = 3:N
        dl[il - 2] = τ[il + 1] - τ[il]  # \Delta x_i
    end
    for im = 2:N
        d[im - 1] = T(2) * (τ[im + 1] - τ[im - 1])  # 2(\Delta x_{i-1} + \Delta x_i)
    end
    for il = 2:(N - 1)
        du[il - 1] = τ[il] - τ[il - 1]  # \Delta x_{i-1}
    end
    for il = 2:N
        fleft = (f[il] - f[il - 1]) / (τ[il] - τ[il - 1])  # f[x_{i-1},x_i]
        fright = (f[il + 1] - f[il]) / (τ[il + 1] - τ[il])  # f[x_i,x_{i+1}]
        # 3(f[x_{i-1},x_i]\Delta x_i+f[x_i,x_{i+1}]\Delta x_{i-1})
        rhs[il - 1] = T(3) * (fleft * (τ[il + 1] - τ[il]) + fright * (τ[il] - τ[il - 1]))
    end
    # The known derivatives at the boundaries get moved from the left-hand side to the RHS.
    rhs[1] -= fp[1] * (τ[3] - τ[2])  # i = 2, (\Delta x_i) s_{i-1}
    rhs[N - 1] -= fp[end] * (τ[N] - τ[N - 1])  # i = N, \Delta x_{i-1} s_{i+1}
    tmat = Tridiagonal(dl, d, du)
    fpans = tmat \ rhs
    fp[2:(end - 1)] .= fpans
    nothing
end

"""
Given an abcissa, τ, of length N+1 values f of length N+1, and derivatives s
of length N+1, this computes the 4xN matrix c of cubic coefficients.
Conte and deBoor Eq. 4.55
"""
function cubic_spline_coefficients!(τ, f::AbstractVector{T}, s, c) where {T <: Real}
    for i = 1:(length(τ) - 1)
        c[1, i] = f[i]
        c[2, i] = s[i]
        Δ = τ[i + 1] - τ[i]
        dd = (f[i + 1] - f[i]) / Δ
        c4p = (s[i + 1] + s[i] - T(2) * dd) / Δ
        c[3, i] = (dd - s[i]) / Δ - c4p
        c[4, i] = c4p / Δ
    end
    nothing
end

"""
"Accurate Monotonicity Preserving Cubic Interpolation" by James M Hyman. 1983.
"""
function deboor_swartz_criterion(s::T, sm1, sp1) where {T <: Real}
    smin = min(sm1, sp1)
    smax = max(sm1, sp1)
    # Equation 2.3
    if zero(T) < smin
        s = min(max(zero(T), s), T(3) * smin)
    elseif smax < zero(T)
        s = max(min(zero(T), s), T(3) * smax)
    elseif sm1 * sp1 <= zero(T)
        s = zero(T)
    else
        s  # Data are locally monotone and within bounds.
    end
end

"""
This comes from Hyman's paper but is modified by the R implementation of splinefun,
which Simon Wood wrote. I didn't understand from the paper the exact structure of
decisions to take in the function.
"""
function hyman_criterion(s::T, sm1, sp1) where {T <: Real}
    # Equation 2.6 "extends" equation 2.3.
    # The paper doesn't specify setting s=sp1 _before_ the if-then about the sign of s.
    σ = s
    if sm1 * sp1 > 0
        σ = sp1  # The paper seems to say this should be set to zero.
    end
    if σ >= 0
        s = min(max(0, s), T(3) * min(abs(sm1), abs(sp1)))
    else
        s = max(min(0, s), -T(3) * min(abs(sm1), abs(sp1)))
    end
    s
end

"""
Uses Hyman piecewise monotonicity.
"""
function project_to_monotonicity!(x, f, s)
    s[1] = hyman_criterion(s[1], (f[2] - f[1])/(x[2] - x[1]), (f[2] - f[1])/(x[2] - x[1]))
    for i in 2:(length(s)-1)
        s[i] = hyman_criterion(
            s[i],
            (f[i] - f[i-1])/(x[i] - x[i-1]),
            (f[i+1] - f[i])/(x[i+1] - x[i])
            )
    end
    s[end] = hyman_criterion(
        s[end],
        (f[end] - f[end-1])/(x[end]-x[end-1]),
        (f[end] - f[end-1])/(x[end]-x[end-1])
        )
    nothing
end

"""
This is how splinefun sets the coefficients. Seems to have the same
result as the code above.
"""
function hyman_coefficients!(τ, f::AbstractVector{T}, s, c) where {T <: Real}
    for i = 1:(length(τ) - 1)
        c[1, i] = f[i]
        c[2, i] = s[i]
        if i < length(τ)
            Δ = τ[i + 1] - τ[i]
            y = -(f[i + 1] - f[i])
            c[3, i] = -(3y + (2s[i] + s[i+1])Δ) / Δ^2
            c[4, i] = (2y/Δ + s[i] + s[i+1])/Δ^2
        else
            # That's correct. This doesn't get called, but R's splinefun sets these.
            # Leaving here until I figure out why we need a polynomial defined for the
            # region to the right of the last input point.
            Δ = τ[i] - τ[i-1]
            y = -(f[i] - f[i-1])
            c[3, i] = (3y + (s[i-1] + 2s[i])Δ) / Δ^2
            c[4, i] = c[4, i-1]
        end
    end
    nothing
end

# Options for endpoints
struct ZeroDerivativeEndpoints
end

struct FlatEndpoints
end

# Options for monotonicity.
struct FreeSlope
end

function project_slope!(::FreeSlope, τ, f, s)
end

struct Monotonic
end

function project_slope!(::Monotonic, τ, f, s)
    project_to_monotonicity!(τ, f, s)
end

# Options for determination of derivatives
struct GlobalDerivatives
end

function cubic_spline(
    τ::AbstractVector{X}, f::AbstractVector{T}, fp; slope=FreeSlope()
    ) where {X <: Real, T <: Real}

    cubic_spline(τ, f, convert(typeof(f), fp); slope = slope)
end

function cubic_spline(
    τ::AbstractVector{X}, f::AbstractVector{T}, fp::AbstractVector{T};
    slope=FreeSlope()
    ) where {X <: Real, T <: Real}

    N = length(τ) - 1
    s = zeros(T, N + 1)
    s[1] = fp[1]
    s[end] = fp[2]
    global_derivatives!(τ, f, s)
    project_slope!(slope, τ, f, s)
    c = zeros(T, 4, N)
    cubic_spline_coefficients!(τ, f, s, c)
    CubicSpline{X,T}(τ, c)
end


text:
const twopi_df64  = Double64(6.283185307179586, 2.4492935982947064e-16)
const onepi_df64  = Double64(3.141592653589793, 1.2246467991473532e-16)
const halfpi_df64 = Double64(1.5707963267948966, 6.123233995736766e-17)
const qrtrpi_df64 = Double64(0.7853981633974483, 3.061616997868383e-17)

const pi_1o1_t64_hi = pi_1o1_t64[1]
const pi_1o1_t64_md = pi_1o1_t64[2]
const pi_1o1_t64_lo = pi_1o1_t64[3]

const pi_1o4_t64_hi = pi_1o4_t64[1]
const pi_1o4_t64_md = pi_1o4_t64[2]
const pi_1o4_t64_lo = pi_1o4_t64[3]

function value_minus_pi(x::DoubleFloat{Float64})
    hi, lo = sub232(HI(x), LO(x), pi_1o1_t64_hi, pi_1o1_t64_md, pi_1o1_t64_lo)
    return DoubleFloat{Float64}(hi, lo)
end

function value_minus_qrtrpi(x::DoubleFloat{Float64})
    hi, lo = sub232(HI(x), LO(x), pi_1o4_t64_hi, pi_1o4_t64_md, pi_1o4_t64_lo)
    return DoubleFloat{Float64}(hi, lo)
end

function mod2pi(x::DoubleFloat{Float64})
    signbit(x) && return mod2pi_neg(x)
    x < twopi_df64 && return x
    w1 = mul323(inv_pi_2o1_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi_2o1_t64, w2)
    z = Double64(y)
    return z
end

function mod2pi_neg(x::DoubleFloat{Float64})
    m = mod2pi(-x)
    return Double64(sub322(pi_2o1_t64, HILO(m)))
end

mod2pi(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(mod2pi(DoubleFloat{Float64}(x)))
mod2pi(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(mod2pi(DoubleFloat{Float64}(x)))

function modpi(x::DoubleFloat{Float64})
    signbit(x) && return modpi_neg(x)
    x < onepi_df64 && return x
    w1 = mul323(inv_pi1o1_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi1o1_t64, w2)
    z = Double64(y)
    return z
end

function modpi_neg(x::DoubleFloat{Float64})
    m = modpi(-x)
    return Double64(sub322(pi1o1_t64, HILO(m)))
end
modpi(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modpi(DoubleFloat{Float64}(x)))
modpi(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modpi(DoubleFloat{Float64}(x)))

function modpipm(x::DoubleFloat{Float64})
    abs(x) < onepi_df64 && return x
    w1 = mul323(inv_pi1o1_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi1o1_t64, w2)
    z = Double64(y)
    return z
end

modpipm(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modpipm(DoubleFloat{Float64}(x)))
modpipm(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modpipm(DoubleFloat{Float64}(x)))

function modhalfpi(x::DoubleFloat{Float64})
    signbit(x) && return modhalfpi_neg(x)
    x < halfpi_df64 && return x
    w1 = mul323(inv_pi1o2_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi1o2_t64, w2)
    z = Double64(y)
    return z
end

function modhalfpi_neg(x::DoubleFloat{Float64})
    m = modhalfpi(-x)
    return Double64(sub322(pi1o2_t64, HILO(m)))
end

modhalfpi(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modhalfpi(DoubleFloat{Float64}(x)))
modhalfpi(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modhalfpi(DoubleFloat{Float64}(x)))

function modhalfpipm(x::DoubleFloat{Float64})
    abs(x) < halfpi_df64 && return x
    w1 = mul323(inv_pi1o2_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi1o2_t64, w2)
    z = Double64(y)
    return z
end

modhalfpipm(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modhalfpipm(DoubleFloat{Float64}(x)))
modhalfpipm(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modhalfpipm(DoubleFloat{Float64}(x)))

function modqrtrpi(x::DoubleFloat{Float64})
    signbit(x) && return modqrtrpi_neg(x)
    x < qrtrpi_df64 && return x
    w1 = mul323(inv_pi1o4_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi1o4_t64, w2)
    z = Double64(y)
    return z
end

function modqrtrpi_neg(x::DoubleFloat{Float64})
    m = modqrtrpi(-x)
    return Double64(sub322(pi1o4_t64, HILO(m)))
end

modqrtrpi(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modqrtrpi(DoubleFloat{Float64}(x)))
modqrtrpi(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modqrtrpi(DoubleFloat{Float64}(x)))

function modqrtrpipm(x::DoubleFloat{Float64})
    abs(x) < qrtrpi_df64 && return x
    w1 = mul323(inv_pi1o4_t64, HILO(x))
    w2 = add_2(w1[1] - trunc(Int,w1[1]), w1[2], w1[3])
    y = mul322(pi1o4_t64, w2)
    z = Double64(y)
    return z
end

modqrtrpipm(x::DoubleFloat{Float32}) = DoubleFloat{Float32}(modqrtrpipm(DoubleFloat{Float64}(x)))
modqrtrpipm(x::DoubleFloat{Float16}) = DoubleFloat{Float16}(modqrtrpipm(DoubleFloat{Float64}(x)))

#=
    rem2pi(x) =  x - 2π*round(x/(2π),r)
    rem2pi(x, RoundDown) == mod2pi(x)

    •    if r == RoundNearest, then the result is in the interval [-π, π]. This will
        generally be the most accurate result.

    •    if r == RoundToZero, then the result is in the interval [0, 2π] if x is positive,.
        or [-2π, 0] otherwise.

    •    if r == RoundDown, then the result is in the interval [0, 2π].

    •    if r == RoundUp, then the result is in the interval [-2π, 0].
=#

rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Down}) = mod2pi(x)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Up}) = -rem2pi(-x, RoundDown)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:Nearest}) = modpipm(x)
rem2pi(x::DoubleFloat{Float64}, rounding::RoundingMode{:ToZero}) =
    signbit(x) ? rem2pi(x, RoundUp) : rem2pi(x, RoundDown)

rem2pi(x::DoubleFloat{Float32}, rounding::RoundingMode) =
    DoubleFloat{Float32}(rem2pi(DoubleFloat{Float64}(x), rounding))
rem2pi(x::DoubleFloat{Float16}, rounding::RoundingMode) =
    DoubleFloat{Float16}(rem2pi(DoubleFloat{Float64}(x), rounding))


text:
"""
    Fixed{T} <: JointConstraint{T}

    fixed connection between two bodies
"""
Fixed(pbody::Node{T}, cbody::Node{T};
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T})) where T =
    Translational{T,3}(pbody, cbody;
        parent_vertex,
        child_vertex),
    Rotational{T,3}(pbody, cbody;
        axis_offset)

"""
    Prismatic{T} <: JointConstraint{T}

    one translational degree of freedom between two bodies
"""
Prismatic(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,1),
    tra_joint_limits=[szeros(T,0), szeros(T,0)]) where T =
    Translational{T,2}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,3}(pbody, cbody;
        axis_offset,
        spring,
        damper)

"""
    Planar{T} <: JointConstraint{T}

    two translational degree of freedom between two bodies
"""
Planar(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,2),
    tra_joint_limits=[szeros(T,0), szeros(T,0)]) where T =
    Translational{T,1}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,3}(pbody, cbody;
        axis_offset,
        spring,
        damper)

"""
    FixedOrientation{T} <: JointConstraint{T}

    three translational degree of freedom between two bodies
"""
FixedOrientation(pbody::Node{T}, cbody::Node{T};
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,3),
    tra_joint_limits=[szeros(T,0), szeros(T,0)]) where T =
    Translational{T,0}(pbody, cbody;
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,3}(pbody, cbody;
        axis_offset,
        spring,
        damper)

"""
    Revolute{T} <: JointConstraint{T}

    one rotational degree of freedom between two bodies
"""
Revolute(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    rot_spring_offset=szeros(T,1),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,3}(pbody, cbody;
        parent_vertex,
        child_vertex,
        spring,
        damper),
    Rotational{T,2}(pbody, cbody;
        axis,
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    Cylindrical{T} <: JointConstraint{T}

    one translational and one rotational degree of freedom between two bodies
"""
Cylindrical(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,1),
    rot_spring_offset=szeros(T,1),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,2}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,2}(pbody, cbody;
        axis,
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    PlanarAxis{T} <: JointConstraint{T}

    two translational and one rotational degree of freedom between two bodies
"""
PlanarAxis(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,2),
    rot_spring_offset=szeros(T,1),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,1}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,2}(pbody, cbody;
        axis,
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    FreeRevolute{T} <: JointConstraint{T}

    free translation with rotation along one axis
"""
FreeRevolute(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,3),
    rot_spring_offset=szeros(T,1),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,0}(pbody, cbody;
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,2}(pbody, cbody;
        axis,
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    Orbital{T} <: JointConstraint{T}

    two rotational degrees of freedom between two bodies
"""
Orbital(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    rot_spring_offset=szeros(T,2),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,3}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper),
    Rotational{T,1}(pbody, cbody;
        axis,
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    PrismaticOrbital{T} <: JointConstraint{T}

    one translational and two rotational degrees of freedom between two bodies
"""
PrismaticOrbital(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,1),
    rot_spring_offset=szeros(T,2),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,2}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,1}(pbody, cbody;
        axis,
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    PlanarOrbital{T} <: JointConstraint{T}

    two translational and two rotational degrees of freedom between two bodies
"""
PlanarOrbital(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,2),
    rot_spring_offset=szeros(T,2),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,1}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,1}(pbody, cbody;
        axis,
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    FreeOrbital{T} <: JointConstraint{T}

    three translational and two rotational degrees of freedom between two bodies
"""
FreeOrbital(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,3),
    rot_spring_offset=szeros(T,2),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,0}(pbody, cbody;
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,1}(pbody, cbody;
        axis,
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    Spherical{T} <: JointConstraint{T}

    three rotational degrees of freedom between two bodies
"""
Spherical(pbody::Node{T}, cbody::Node{T};
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    axis_offset=one(Quaternion{T}),
    spring=zero(T),
    damper=zero(T),
    rot_spring_offset=szeros(T,3),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,3}(pbody, cbody;
        parent_vertex,
        child_vertex,
        spring,
        damper),
    Rotational{T,0}(pbody, cbody;
        axis_offset,
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    CylindricalFree{T} <: JointConstraint{T}

    one translational and three rotational degrees of freedom between two bodies
"""
CylindricalFree(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    spring=zero(T), damper=zero(T),
    tra_spring_offset=szeros(T,1),
    rot_spring_offset=szeros(T,3),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,2}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,0}(pbody, cbody;
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    PlanarFree{T} <: JointConstraint{T}

    two translational and three rotational degrees of freedom between two bodies
"""
PlanarFree(pbody::Node{T}, cbody::Node{T}, axis;
    parent_vertex=szeros(T, 3),
    child_vertex=szeros(T, 3),
    spring=zero(T),
    damper= zero(T),
    tra_spring_offset=szeros(T,2),
    rot_spring_offset=szeros(T,3),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,1}(pbody, cbody;
        parent_vertex,
        child_vertex,
        axis,
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,0}(pbody, cbody;
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

"""
    Floating{T} <: JointConstraint{T}

    no restricted degrees of freedom between two bodies
"""
Floating(pbody::Node{T}, cbody::Node{T};
    spring=zero(T),
    damper=zero(T),
    tra_spring_offset=szeros(T,3),
    rot_spring_offset=szeros(T,3),
    rot_joint_limits=[szeros(T,0), szeros(T,0)],
    tra_joint_limits=[szeros(T,0), szeros(T,0)],
    spring_type=:linear) where T =
    Translational{T,0}(pbody, cbody;
        spring,
        damper,
        spring_offset=tra_spring_offset,
        joint_limits=tra_joint_limits),
    Rotational{T,0}(pbody, cbody;
        spring,
        damper,
        spring_offset=rot_spring_offset,
        joint_limits=rot_joint_limits,
        spring_type=spring_type)

function Prototype(joint_type::Symbol, pbody::Node{T}, cbody::Node{T}, axis;
        parent_vertex=szeros(T, 3),
        child_vertex=szeros(T, 3),
        axis_offset=one(Quaternion{T}),
        spring=zero(T),
        damper=zero(T),
        tra_spring_offset=nothing,
        rot_spring_offset=nothing,
        tra_joint_limits=[szeros(T,0), szeros(T,0)],
        rot_joint_limits=[szeros(T,0), szeros(T,0)],
        spring_type=:linear) where T

    N̄tra, N̄rot = nullspace_dimension(joint_type)
    (tra_spring_offset == nothing) && (tra_spring_offset = szeros(T,N̄tra))
    (rot_spring_offset == nothing) && (rot_spring_offset = szeros(T,N̄rot))
    (joint_type == :Fixed)            && (return            Fixed(pbody, cbody;       parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset))
    (joint_type == :Prismatic)        && (return        Prismatic(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset,                                      tra_joint_limits=tra_joint_limits))
    (joint_type == :Planar)           && (return           Planar(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset,                                      tra_joint_limits=tra_joint_limits))
    (joint_type == :FixedOrientation) && (return FixedOrientation(pbody, cbody;                     axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset,                                      tra_joint_limits=tra_joint_limits))
    (joint_type == :Revolute)         && (return         Revolute(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper,                                      rot_spring_offset=rot_spring_offset,                                    rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :Cylindrical)      && (return      Cylindrical(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :PlanarAxis)       && (return       PlanarAxis(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :FreeRevolute)     && (return     FreeRevolute(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :Orbital)          && (return          Orbital(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper,                                      rot_spring_offset=rot_spring_offset,                                    rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :PrismaticOrbital) && (return PrismaticOrbital(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :PlanarOrbital)    && (return    PlanarOrbital(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :FreeOrbital)      && (return      FreeOrbital(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :Spherical)        && (return        Spherical(pbody, cbody;       parent_vertex=parent_vertex, child_vertex=child_vertex, axis_offset=axis_offset, spring=spring, damper=damper,                                      rot_spring_offset=rot_spring_offset,                                    rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :CylindricalFree)  && (return  CylindricalFree(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex,                  spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :PlanarFree)       && (return       PlanarFree(pbody, cbody, axis; parent_vertex=parent_vertex, child_vertex=child_vertex,                  spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
    (joint_type == :Floating)         && (return         Floating(pbody, cbody;                                      spring=spring, damper=damper, tra_spring_offset=tra_spring_offset, rot_spring_offset=rot_spring_offset, tra_joint_limits=tra_joint_limits, rot_joint_limits=rot_joint_limits, spring_type=spring_type))
end

function nullspace_dimension(joint_type::Symbol)
    (joint_type == :Fixed)            && (return 0, 0)
    (joint_type == :Prismatic)        && (return 1, 0)
    (joint_type == :Planar)           && (return 2, 0)
    (joint_type == :FixedOrientation) && (return 3, 0)
    (joint_type == :Revolute)         && (return 0, 1)
    (joint_type == :Cylindrical)      && (return 1, 1)
    (joint_type == :PlanarAxis)       && (return 2, 1)
    (joint_type == :FreeRevolute)     && (return 3, 1)
    (joint_type == :Orbital)          && (return 0, 2)
    (joint_type == :PrismaticOrbital) && (return 1, 2)
    (joint_type == :PlanarOrbital)    && (return 2, 2)
    (joint_type == :FreeOrbital)      && (return 3, 2)
    (joint_type == :Spherical)        && (return 0, 3)
    (joint_type == :CylindricalFree)  && (return 1, 3)
    (joint_type == :PlanarFree)       && (return 2, 3)
    (joint_type == :Floating)         && (return 3, 3)
end


text:
function get_users()
    active_users = []
    user_paths = glob(joinpath(Pinguu.base_dir[2:end], "users/*"));
    for path in user_paths
        push!(active_users, split(path, "/")[end])
    end
    user_paths = glob(joinpath(Pinguu.base_dir[2:end], "invest/users/*"));
    for path in user_paths
        push!(active_users, split(path, "/")[end])
    end
    return unique(active_users)
end
return get_users


text:
const obsnames = ["Energy", "Energy^2", "Specific Heat",
                  "|Magnetization|", "Magnetization^2", "Magnetization^4",
                  "Susceptibility", "Connected Susceptibility", "Binder Ratio",
                 ]
const Ns = [2,8]
const Qs = [2,4]
const J = 1.0
const Ts = [0.3, 1.0, 3.0, 10.0]

"""
exact calculation for `Q` state Potts model on `N` site fully connected network
"""
function exact(Q, J, N, Ts)
    j = J/N
    Es = zeros(0)
    Ms = zeros(0)
    spins = zeros(Int, N)
    nst = Q^N
    for ist in 0:(nst-1)
        spins .= 1
        for i in 1:N
            spins[i] = mod(ist,Q)+1
            ist = div(ist,Q)
        end
        E = 0.0
        M = 0.0
        Qinv = 1.0/Q
        nq = count(s->s==1, spins)
        E -= 0.5j*(nq*(nq-1))
        M += nq*(1.0-Qinv)
        for q in 2:Q
            nq = count(s->s==q, spins)
            E -= 0.5j*(nq*(nq-1))
            M -= nq*Qinv
        end
        M /= N
        push!(Es, E)
        push!(Ms, M)
    end
    sorted = sortperm(Es,rev=true)
    E0 = Es[sorted[end]]

    mbetas = Ts .\ (-1.0)
    nT = length(Ts)
    res = Dict(n=>zeros(nT) for n in obsnames)
    res["Magnetization"] = zeros(nT)
    Zs = zeros(nT)
    for i in sorted
        E = Es[i] - E0
        zs = exp.(mbetas.*E)
        Zs .+= zs
        E += E0
        res["Energy"] .+= zs.*(E/N)
        res["Energy^2"] .+= zs.*(E/N)^2

        M = Ms[i]
        res["Magnetization"] .+= zs.*M
        res["|Magnetization|"] .+= zs.*abs(M)
        res["Magnetization^2"] .+= zs.*M^2
        res["Magnetization^4"] .+= zs.*M^4
    end
    res["Magnetization"] ./= Zs
    for n in obsnames
        res[n] ./= Zs
    end
    return res
end

betas = 1.0 ./ Ts

for N in Ns
    for Q in Qs
        res = exact(Q,J,N,Ts)
        res["Specific Heat"] = N.*(res["Energy^2"] .- res["Energy"].^2).*betas.^2
        res["Susceptibility"] = N.*(res["Magnetization^2"] .- res["Magnetization"].^2).*betas
        res["Connected Susceptibility"] = N.*(res["Magnetization^2"] .- res["|Magnetization|"].^2).*betas
        res["Binder Ratio"] = res["Magnetization^4"] ./ (res["Magnetization^2"].^2)
        open(@sprintf("Q_%d__J_%.1f__N_%d.dat", Q, J, N), "w") do io
            for i in 1:length(Ts)
                @printf(io, "%.15f", Ts[i])
                for n in obsnames
                    @printf(io, " %.15f", res[n][i])
                end
                println(io)
            end
        end
    end
end


text:
using DashBootstrapComponents, DashHtmlComponents

modal = html_div([
    dbc_button("Open", id="open-centered"),
    dbc_modal(
        [
            dbc_modalheader("Header"),
            dbc_modalbody("This modal is vertically centered"),
            dbc_modalfooter(
                dbc_button(
                    "Close",
                    id="close-centered",
                    className="ml-auto",
                    n_clicks=0,
                ),
            ),
        ],
        id="modal-centered",
        centered=true,
        is_open=false,
    ),
]);

callback!(
    app,
    Output("modal-centered", "is_open"),
    Input("open-centered", "n_clicks"),
    Input("close-centered", "n_clicks"),
    State("modal-centered", "is_open"),
) do n1, n2, is_open
    return n1 > 0 || n2 > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module TestArrays

using UnionArrays
using Test

function test()
    A = reshape(UnionVector(Any[1, 2.0]), 1, :)
    @test A isa UnionArray
    @test A[1] === 1
    @test A[2] === 2.0
    @test A[1, 1] === 1
    @test A[1, 2] === 2.0
    @test (A[1, 1] = 3.0) isa Any
    @test A[1, 1] === 3.0
    @test (A[1] = 3.0) isa Any
    @test A[1] === 3.0
end

end  # module


text:
```

ДАНО: Робот - в юго-западном углу поля, на котором расставлено некоторое количество маркеров

РЕЗУЛЬТАТ: Функция вернула значение средней температуры всех замаркированных клеток

```

using HorizonSideRobots

show!(originRobot)
include("../AbstractType.jl")
include("../functions.jl")

originRobot = Robot(animate = false, "10example.sit")
mod_CmR = MODcountmarkersRobot(originRobot,0,0)

function moveLikeSnake!(robot::MODcountmarkersRobot, moveSides::NTuple{2,HorizonSide}, borderSide::HorizonSide)::Nothing
    flag = true
    while flag
        for side in moveSides
            while !(isborder(robot.robot, side))
                move!(robot, side)
            end
            if !(isborder(robot.robot, borderSide))
                move!(robot, borderSide)
            else
                flag = false
            end
        end
    end
end

# master главная функция программа
function master!(robot::MODcountmarkersRobot)::Real
    moveLikeSnake!(robot, (Nord, Sud), Ost)
    return (robot.temperature / robot.markers)
end

master!(mod_CmR)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
# Analyse results, used to print posterior inference results

#using Plots
using PyPlot
using KernelDensity
using Distributions

using DataFrames

include(pwd()*"/utilities/posteriorinference.jl")

remove_missing_values(x) = reshape(collect(skipmissing(x)),7,:)

problem = "real data"
problem = "sim data small problem"

algorithm = "MCWM"
algorithm = "DA"
algorithm = "ADA"

if problem == "real data" && algorithm == "MCWM"
    jobname = "gp_training_7_par_lunarc_new_data_4_coresnew_data"
elseif problem == "real data" && algorithm == "DA"
    jobname = "_dagpest7new_datada_gp_mcmc"
elseif problem == "real data" && algorithm == "ADA"
    jobname = "_dagpest7new_dataada_gp_mcmc_dt"
elseif problem == "sim data small problem" && algorithm == "MCWM"
    jobname = "gp_training_7_par_lunarc_simdata_4_coressimdata"
elseif problem == "sim data small problem" && algorithm == "DA"
    jobname = "_dagpest7simdatada_gp_mcmc"
else #problem == "sim data scaled up problem" && algorithm = "ADA"
    jobname = "_dagpest7simdataada_gp_mcmc_dt"
end

if problem == "real data"
    plot_theta_true = false
else
    plot_theta_true = true
end

data_res = convert(Array,readtable("DWPSDE model/Results/"*problem*"/output_res"*jobname*".csv"))

M, N = size(data_res)

data_param = convert(Array,readtable("DWPSDE model/Results/"*problem*"/output_param"*jobname*".csv"))

theta_true = data_param[1:N-2]
burn_in = Int64(data_param[N-2+1])

data_prior_dist = convert(Array,readtable("DWPSDE model/Results/"*problem*"/output_prior_dist"*jobname*".csv"))

data_prior_dist_type = convert(Array,readtable("DWPSDE model/Results/"*problem*"/output_prior_dist_type"*jobname*".csv"))
data_prior_dist_type = data_prior_dist_type[2]

Z = convert(Array,readtable("DWPSDE model/Results/"*problem*"/data_used"*jobname*".csv"))
Z = Z[:,1]

if algorithm == "DA" || algorithm == "ADA"
  burn_in = 1
end

Theta = remove_missing_values(data_res[:,1:N-2]') # stor data in column-major order
loglik = data_res[:,N-1]
accept_vec = data_res[:,N]

accept_rate = sum(accept_vec)/M
nbr_acf_lags =  50

acf = ones(N-2,nbr_acf_lags+1)

lower_q_int_limit = 2.5
upper_q_int_limit = 97.5

# use L"$\log r$"

if N == 6
    title_vec_log = [ L"$\log \kappa$"; L"$\log \gamma$"; L"$\log c$"; L"$\log d$"]
    title_vec = [ L"$\kappa$"; L"Gamma"; "c"; "d"]
elseif N == 4
    title_vec_log = [ L"$\log c$"; L"$\log d$" ]
    title_vec = [ "c"; "d" ];
elseif N == 5
    title_vec_log = [ L"$\log A$";L"$\log c$"; L"$\log d$" ]
    title_vec = [ "A";"c"; "d" ]
elseif N == 8
    title_vec_log = [ L"$\log A$"; L"$\log c$"; L"$\log d$"; L"$\log p_1$"; L"$\log p_2$"; L"$\log sigma$"]
    title_vec = [  "A    "; "c    "; "d    "; "p_1  "; "p_1  "; "sigma"]
elseif N == 7
    title_vec_log = [ L"$\log \kappa$"; L"$\log Gamma$"; L"$\log c$"; L"$\log d$"; L"$\log sigma$"];
    title_vec = [ L"$\kappa$"; "Gamma"; "c"; "d"; "sigma"]
elseif N == 9
    title_vec_log = [ L"$\log \kappa$"; L"$\log \gamma$"; L"$\log c$"; L"$\log d$"; L"$\log p_1$"; L"$\log p_2$"; L"$\log \sigma$"]
    title_vec = [  L"$\kappa$"; L"$\gamma$"; L"$c$"; L"$d$"; L"$p_1$"; L"$p_1$"; L"$\sigma$"]
else
    title_vec_log = [ L"$\log \kappa$"; L"$\log Gamma$"; L"$\log A$"; L"$\log c$"; L"$\log d$"; L"$\log g$"; L"$\log p_1$"; L"$\log p_2$"; L"$\log sigma$"]
    title_vec = [ L"$\kappa$"; "Gamma"; "A"; "c"; "d"; "g"; "p_1"; "p_1"; "sigma"]
end

acceptance_rate = sum(accept_vec[burn_in:end])/length(accept_vec[burn_in:end])

# print info
println("---------------------------------------------------")
print("Problem:")
print(problem)
println("")
print("Algorithm:")
print(algorithm)
println("")
println("")

@printf "Accept rate: %.4f %% \n" round(acceptance_rate*100,2)

println("")

@printf "True parameter values:\n"
Base.showarray(STDOUT,round.(theta_true,2),false)
@printf "\n"

println("")

@printf "Posterior mean:\n"
Base.showarray(STDOUT,round.(mean(Theta[:,burn_in+1:end],2),2),false)
@printf "\n"

println("")

@printf "Posterior standard deviation:\n"
Base.showarray(STDOUT,round.(std(Theta[:,burn_in+1:end],2),2),false)
@printf "\n"

println("")

@printf "Posterior quantile intervals (2.5th and 97.5th quantiles as default):\n"
Base.showarray(STDOUT,round.(calcquantileint(Theta[:,burn_in+1:end],lower_q_int_limit,upper_q_int_limit),2),false)
@printf "\n"
println("---------------------------------------------------")

# plot trace plots

text_size = 15
label_size = 10

#on log-scale
PyPlot.figure(figsize=(10,20))
ax = axes()
for i = 1:N-2
    ax1 = PyPlot.subplot(N-2,1,i)
    PyPlot.plot(Theta[i,:])
    plot_theta_true == true ? PyPlot.plot(ones(size(Theta,2),1)*theta_true[i], "k") : 2
    PyPlot.ylabel(title_vec_log[i],fontsize=text_size)
    i < 7 ? ax1[:axes][:get_xaxis]()[:set_ticks]([]) : 2
end
PyPlot.xlabel("Iteration",fontsize=text_size)
ax[:tick_params]("both",labelsize = label_size)

# on non-log scale
PyPlot.figure()
for i = 1:N-2
    PyPlot.subplot(N-2,1,i)
    PyPlot.plot(exp.(Theta[i,:]))
    plot_theta_true == true ? PyPlot.plot(ones(size(Theta,2),1)*exp(theta_true[i]), "k") :
    PyPlot.ylabel(title_vec[i],fontsize=text_size)
end
PyPlot.xlabel("Iteration")

# plot trace plots after burn in

# on log-scale
PyPlot.figure(figsize=(10,20))
ax = axes()
x_axis = burn_in+1:size(Theta,2)
for i = 1:N-2
    ax1 = PyPlot.subplot(N-2,1,i)
    PyPlot.plot(x_axis, Theta[i,burn_in+1:end])
    plot_theta_true == true ? PyPlot.plot(x_axis, ones(length(x_axis),1)*theta_true[i], "k") : 2
    PyPlot.ylabel(title_vec_log[i],fontsize=text_size)
    i < 7 ? ax1[:axes][:get_xaxis]()[:set_ticks]([]) : 2
end
PyPlot.xlabel("Iteration",fontsize=text_size)
ax[:tick_params]("both",labelsize = label_size)

# on non-log scale
PyPlot.figure()
x_axis = burn_in+1:size(Theta,2)
for i = 1:N-2
    PyPlot.subplot(N-2,1,i)
    PyPlot.plot(x_axis, exp.(Theta[i,burn_in+1:end]))
    plot_theta_true == true ? PyPlot.plot(x_axis, ones(length(x_axis),1)*exp(theta_true[i]), "k") :
    PyPlot.ylabel(title_vec[i])
end
PyPlot.xlabel("Iteration")

# calc acf for each marginal chain
for i = 1:N-2
    acf[i,2:end] = autocor(Theta[i,burn_in:end],1:nbr_acf_lags)
end

# plot acf
PyPlot.figure()
for i = 1:N-2
    PyPlot.subplot(N-2,1,i)
    PyPlot.plot(0:nbr_acf_lags, acf[i,:])
    PyPlot.ylabel(title_vec_log[i])
end
PyPlot.xlabel("Lag")

# plot posterior
PyPlot.figure()
for i = 1:N-2
    PyPlot.subplot(N-2,1,i)
    h = kde(Theta[i,burn_in+1:end])
    PyPlot.plot(h.x,h.density, "b")

    if data_prior_dist_type == "Uniform"
        error("Uniform priors are not not implemented.")
    elseif data_prior_dist_type == "Normal"
        x_grid = (data_prior_dist[i,1]-4*data_prior_dist[i,2]):0.001:(data_prior_dist[i,1]+4*data_prior_dist[i,2])
        PyPlot.plot(x_grid, pdf(Normal(data_prior_dist[i,1],data_prior_dist[i,2]), x_grid) ,"g")
    end

    plot_theta_true == true ? PyPlot.plot((theta_true[i], theta_true[i]), (0, maximum(h.density)), "k") :

    PyPlot.ylabel(title_vec_log[i])
end

# plot posterior, non-log scale
PyPlot.figure()
for i = 1:N-2
    PyPlot.subplot(N-2,1,i)
    h = kde(exp.(Theta[i,burn_in+1:end]))
    PyPlot.plot(h.x,h.density, "b")
    plot_theta_true == true ? PyPlot.plot((exp(theta_true[i]), exp(theta_true[i])), (0, maximum(h.density)), "k") :
    PyPlot.ylabel(title_vec[i])
end

# plot log-lik
PyPlot.figure()
PyPlot.plot(loglik)
PyPlot.ylabel(L"$\log-likelhood")
PyPlot.xlabel("Iteration")

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
ans = 2

ans = 3+ans

println(ans)
println("hello world")


text:
using Test

m = m = test_page_model()
include("../src/components/TotalForcing.jl")

totalforcing = add_comp!(m, TotalForcing)

totalforcing[:f_CO2forcing] = readpagedata(m,"test/validationdata/f_co2forcing.csv")
totalforcing[:f_CH4forcing] = readpagedata(m,"test/validationdata/f_ch4forcing.csv")
totalforcing[:f_N2Oforcing] = readpagedata(m,"test/validationdata/f_n2oforcing.csv")
totalforcing[:f_lineargasforcing] = readpagedata(m,"test/validationdata/f_LGforcing.csv")
totalforcing[:exf_excessforcing] = readpagedata(m,"data/exf_excessforcing.csv")

run(m)

forcing=m[:TotalForcing, :ft_totalforcing]
forcing_compare=readpagedata(m,"test/validationdata/ft_totalforcing.csv")

@test forcing ≈ forcing_compare rtol=1e-3


text:
###Functions used for heat and moisture transfer calculation####

function mesh(E,dl)	#Asigning a material key to certain point
	id=1
	dim=Array{Int64}(undef,size(E)[1],2)
	for i=1:size(E)[1]
		xmin=minimum([minimum(E[i,2]),minimum(E[i,4])])
		xmax=maximum([maximum(E[i,2]),maximum(E[i,4])])
		ymin=minimum([minimum(E[i,3]),minimum(E[i,5])])
		ymax=maximum([maximum(E[i,3]),maximum(E[i,5])])
		x=xmin
		y=ymin
		eps=10e-7
		xd=(xmax-xmin)
		yd=(ymax-ymin)
		nx=xd/dl
		nx=round(nx)
		if abs(nx-round(nx))>eps
			println("ERROR! Incorrect subdivision!") #Subdivision is constant and divisible with element's xd and yd
			break
		end
		nx=floor(Int64,nx)
		ny=yd/dl
		ny=round(ny)
		if abs(ny-round(ny))>eps
			println("ERROR! Incorrect subdivision!")
			break
		end
		ny=floor(Int64,ny)
		dim[i,1]=nx
		dim[i,2]=ny
		y=ymin+dl/2
		Md=zeros(nx*ny,8)
		for j=1:ny
			x=xmin+dl/2
			for k=1:nx
				#Md=[Md;round(x,digits=3) round(y,digits=3) E[i,1]]
				Md[(j-1)*nx+k,1]=id
				Md[(j-1)*nx+k,2]=(j-1)*nx+k
				Md[(j-1)*nx+k,3]=x
				Md[(j-1)*nx+k,4]=y
				Md[(j-1)*nx+k,5]=E[i,1]
				Md[(j-1)*nx+k,6]=E[i,7]
				Md[(j-1)*nx+k,7]=E[i,8]
				Md[(j-1)*nx+k,8]=E[i,9]
				x=x+dl
				id=id+1
			end
			y=y+dl
		end
		writedlm("mesh/mesh$i.csv",  Md, ',')
	end
	M2=readdlm("mesh/mesh1.csv",',')
	for i=2:size(E)[1]
		m2=readdlm("mesh/mesh$i.csv",',')
		M2=[M2;m2]
	end
	writedlm("mesh/mesh.csv",  M2, ',')
	#for i=1:16
		#rm("mesh/mesh$i.csv", recursive=true)
	#end
	Md=copy(M2)
	Md=[["id_glob" "id_loc" "x_coord" "y_coord" "material" "lambda" "rho" "cp"];Md]
	Z0_vect=zeros((size(Md)[1]-1))
	Z0_vect=["z_coord";Z0_vect]
	writedlm("mesh/paraview_mesh.csv",  [Md Z0_vect], ',')
	writedlm("mesh/dim.csv",  dim, ',')
	M2=M2
end
####End of function mesh(E,dl)####
function nsort(E,M,dl)
	M_ele=Vector{Any}(undef,size(E)[1])
	for i=1:size(E)[1]
		M_ele[i,1]=readdlm("mesh/mesh$i.csv",',')
	end
	#for i=1:size(E)[1]
		#rm("mesh/mesh$i.csv", recursive=true)
	#end
	Nsort=M_ele
end
####End of function nsort(E,M,dl)####
function mappoint(E,Ns,x1,y1)	#Connection between nodes
	ind1=[1]
	for k=1:size(E)[1]
		if E[k,1]>99
			ind1[1]=k
			break
		end
	end
	ind1=ind1[1]#-1
	k=floor(Int64,ind1)
	el_id=[1]
	p=[x1 y1]
	for i=1:k
		for j=1:size(E)[1]
			if j!=i
				if (p[1]>E[j,2])&(p[1]<E[j,4])&(p[2]>E[j,3])&(p[2]<E[j,5])
					el_id[1]=j
				elseif (p[1]>E[j,2])&(p[1]<E[j,4])&(p[2]>E[j,3])&(p[2]<E[j,5])
					el_id[1]=j
				elseif (p[1]>E[j,2])&(p[1]<E[j,4])&(p[2]>E[j,3])&(p[2]<E[j,5])
					el_id[1]=j
				elseif (p[1]>E[j,2])&(p[1]<E[j,4])&(p[2]>E[j,3])&(p[2]<E[j,5])
					el_id[1]=j
				end
			end
		end
	end
	el_id=floor(Int64,el_id[1])
end
####End of function con2(E,M,dl,Ns)####
function elecon(xp,yp,E,Ns,dl,neig,mat)#(idp,E,Ns,dl,neig,mat)
	#xp=mat[idp,3]
	#yp=mat[idp,4]
	eps=10e-6
	glob_id=[1]
	dim=readdlm("mesh/dim.csv",',')
	#left neighbour
	if neig==2
		xl=xp-dl
		yl=yp
		el_id=mappoint(E,Ns,xl,yl)
		mat2=Ns[el_id]
		nx=floor(Int64,dim[el_id,1])
		ny=floor(Int64,dim[el_id,2])
		for zz=1:ny
			if (abs(xl-mat2[zz*nx,3])<eps)&(abs(yl-mat2[zz*nx,4])<eps)
				glob_id[1]=mat2[zz*nx,1]
			end
		end
	end
	#right
	if neig==3
		xr=xp+dl
		yr=yp
		el_id=mappoint(E,Ns,xr,yr)
		mat2=Ns[el_id]
		nx=floor(Int64,dim[el_id,1])
		ny=floor(Int64,dim[el_id,2])
		for zz=1:ny
			if (abs(xr-mat2[(zz-1)*nx+1,3])<eps)&(abs(yr-mat2[(zz-1)*nx+1,4])<eps)
				glob_id[1]=mat2[(zz-1)*nx+1,1]
			end
		end
	end
	#down
	if neig==4
		xd=xp
		yd=yp-dl
		el_id=mappoint(E,Ns,xd,yd)
		mat2=Ns[el_id]
		nx=floor(Int64,dim[el_id,1])
		ny=floor(Int64,dim[el_id,2])
		for zz=1:nx
			if (abs(xd-mat2[(ny-1)*nx+zz,3])<eps)&(abs(yd-mat2[(ny-1)*nx+zz,4])<eps)
				glob_id[1]=mat2[(ny-1)*nx+zz,1]
			end
		end
	end
	#up
	if neig==5
		xu=xp
		yu=yp+dl
		el_id=mappoint(E,Ns,xu,yu)
		mat2=Ns[el_id]
		nx=floor(Int64,dim[el_id,1])
		ny=floor(Int64,dim[el_id,2])
		for zz=1:nx
			if (abs(xu-mat2[zz,3])<eps)&(abs(yu-mat2[zz,4])<eps)
				glob_id[1]=mat2[zz,1]
			end
		end
	end
	glob_id=floor(Int64,glob_id[1])
end
####End of function elecon(x,y,M,mat)####
function con1(E,M,dl,Ns)	#Connection between nodes
	ind1=[1 1]
	for k=1:size(E)[1]
		if E[k,1]>99
			ind1[1]=Ns[k][1,1]-1
			ind1[2]=k
			break
		end
	end
	ind=floor(Int64,ind1[1])
	k=floor(Int64,ind1[2]-1)
	C=Array{Int64}(undef,ind,5)
	dim=readdlm("mesh/dim.csv",',')
	C[:,1]=M[1:ind,1]
	for i=1:k
		mat=Ns[i]
		nx=floor(Int64,dim[i,1])
		ny=floor(Int64,dim[i,2])
		for j=1:nx*ny
			if mod((j-1),nx)==0
				left=2
				C[floor(Int64,mat[j,1]),2]=elecon(mat[j,3],mat[j,4],E,Ns,dl,left,mat)
				C[floor(Int64,mat[j,1]),3]=M[floor(Int64,mat[j,1])+1,1]
				if j>1
					C[floor(Int64,mat[j,1]),4]=M[floor(Int64,mat[j,1])-nx,1]
				end
				if j<nx*ny#(nx*ny-1)
					C[floor(Int64,mat[j,1]),5]=M[floor(Int64,mat[j,1])+nx-1,1]+1
				end
			elseif mod(j,nx)==0
				right=3
				C[floor(Int64,mat[j,1]),3]=elecon(mat[j,3],mat[j,4],E,Ns,dl,right,mat)
				C[floor(Int64,mat[j,1]),2]=M[floor(Int64,mat[j,1])-1,1]
				if j<nx*ny
					C[floor(Int64,mat[j,1]),5]=M[floor(Int64,mat[j,1])+nx,1]
				end
				if j>nx
					C[floor(Int64,mat[j,1]),4]=M[floor(Int64,mat[j,1])-nx,1]
				end
			end
			if j<=nx
				down=4
				C[floor(Int64,mat[j,1]),4]=elecon(mat[j,3],mat[j,4],E,Ns,dl,down,mat)
			elseif j>=ny*nx-nx+1
				up=5
				C[floor(Int64,mat[j,1]),5]=elecon(mat[j,3],mat[j,4],E,Ns,dl,up,mat)
			end
			if (mod((j-1),nx)!=0)&(mod(j,nx)!=0)&(j>nx)&(j<ny*nx-nx)
				C[floor(Int64,mat[j,1]),2]=M[floor(Int64,mat[j,1])-1,1]
				C[floor(Int64,mat[j,1]),3]=M[floor(Int64,mat[j,1])+1,1]
				C[floor(Int64,mat[j,1]),4]=M[floor(Int64,mat[j,1])-nx,1]
				C[floor(Int64,mat[j,1]),5]=M[nx+floor(Int64,mat[j,1]),1]
			end
			if (j>1)&(j<nx)
				C[floor(Int64,mat[j,1]),2]=M[floor(Int64,mat[j,1])-1,1]
				C[floor(Int64,mat[j,1]),3]=M[floor(Int64,mat[j,1])+1,1]
				C[floor(Int64,mat[j,1]),5]=M[nx+floor(Int64,mat[j,1]),1]
			elseif (j>(ny*nx-nx+1))&(j<nx*ny)
				C[floor(Int64,mat[j,1]),2]=M[floor(Int64,mat[j,1])-1,1]
				C[floor(Int64,mat[j,1]),3]=M[floor(Int64,mat[j,1])+1,1]
				C[floor(Int64,mat[j,1]),4]=M[floor(Int64,mat[j,1])-nx,1]
			end
		end
	end
	#writedlm("mesh/Con.csv",  C, ',')
	C=C
end
####End of function con1(E,M,dl,Ns)####
function Kmat(M,C,dl,dt=0) 	#Conductivity matrix
	ind=Con[end,1]
	#K=Array{Any}(ind,5)
	K=Matrix{Any}(undef,ind,5)
	f=Matrix{Any}(undef,ind,2)
	if dt!=0
		for i=1:ind
			j=1
			if M[C[i,2],8]<0
				Kl=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kl=1/(1/(2*M[C[i,2],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,3],8]<0
				Kr=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kr=1/(1/(2*M[C[i,3],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,4],8]<0
				Kd=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kd=1/(1/(2*M[C[i,4],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,5],8]<0
				Ku=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Ku=1/(1/(2*M[C[i,5],6])+1/(2*M[C[i,1],6]))
			end
			a=M[C[i,1],7]*M[C[i,1],8]
			K[i,1]=[Con[i,1] a*dl^2/dt+Kl+Kr+Kd+Ku]
			f[i,1]=[1 0]
			f[i,2]=[1 0]
			if C[i,2]<=ind
				K[i,2]=[C[i,2] -Kl]
			else
				K[i,2]=[1 0]
				f[i,j]=[C[i,2] Kl]
				j=j+1
			end
			if C[i,3]<=ind
				K[i,3]=[C[i,3] -Kr]
			else
				K[i,3]=[1 0]
				f[i,j]=[C[i,3] Kr]
				j=j+1
			end
			if C[i,4]<=ind
				K[i,4]=[C[i,4] -Kd]
			else
				K[i,4]=[1 0]
				f[i,j]=[C[i,4] Kd]
				j=j+1
			end
			if C[i,5]<=ind
				K[i,5]=[C[i,5] -Ku]
			else
				K[i,5]=[1 0]
				f[i,j]=[C[i,5] Ku]
			end
		end
		K=[K f]
	else
		for i=1:ind
			j=1
			if M[C[i,2],8]<0
				Kl=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kl=1/(1/(2*M[C[i,2],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,3],8]<0
				Kr=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kr=1/(1/(2*M[C[i,3],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,4],8]<0
				Kd=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kd=1/(1/(2*M[C[i,4],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,5],8]<0
				Ku=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Ku=1/(1/(2*M[C[i,5],6])+1/(2*M[C[i,1],6]))
			end
			K[i,1]=[Con[i,1] (Kl+Kr+Kd+Ku)]
			f[i,1]=[1 0]
			f[i,2]=[1 0]
			if C[i,2]<=ind
				K[i,2]=[C[i,2] -Kl]
			else
				K[i,2]=[1 0]
				f[i,j]=[C[i,2] Kl]
				j=j+1
			end
			if C[i,3]<=ind
				K[i,3]=[C[i,3] -Kr]
			else
				K[i,3]=[1 0]
				f[i,j]=[C[i,3] Kr]
				j=j+1
			end
			if C[i,4]<=ind
				K[i,4]=[C[i,4] -Kd]
			else
				K[i,4]=[1 0]
				f[i,j]=[C[i,4] Kd]
				j=j+1
			end
			if C[i,5]<=ind
				K[i,5]=[C[i,5] -Ku]
			else
				K[i,5]=[1 0]
				f[i,j]=[C[i,5] Ku]
			end
		end
		K=[K f]
	end
	#writedlm("mesh/K.csv",  K, ',')
	#writedlm("mesh/f.csv",  f, ',')
	K=K
end
####End of function Kmat(M,con,dl)####
function modconjgrad(Km,M,dt=0,dl=0,x=0,in_con=0)
    if in_con==0
        b2=zeros(size(Km)[1],1)
	else
		b2=init(in_con,M,dl,dt)
    end
    if x==0
        x=ones(size(Km)[1],1)#rand(0:20,size(Km)[1],1)#ones(size(Km)[1],1)
    end
    b1=bcf(Km,M)
	b=b1+b2
    r=b-mult(Km,x)
    p=r
    rsold=(r'*r)[1]
	iter=1000#size(Km)[1] #iterations
	i=1
	j=1
    while i<iter#for i=1:10000#size(b)[1]
		@time begin
        Ap=mult(Km,p)
		#end
		#@time begin
        alpha=rsold/(p'*Ap)[1]
		#end
		#@time begin
        x=x+alpha*p
		#end
		#@time begin
		zz=mult(Km,p)
		#end
		#@time begin
		r=r-alpha*zz
		#end
		#@time begin
        rsnew=(r'*r)[1]
		println(sqrt(rsnew))
        if sqrt(rsnew)<1e-7
			println("Convergence 1e-7 achieved -- Number of iterations: $((j-1)*(iter-1)+i)")
			#M2=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]#M2=[["x_coord" "y_coord"];M[:,3:4]]#M=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]
			#Z0_vect=zeros(size(M2)[1]-1)#Z0_vect=zeros(size(M)[1])
			#Z0_vect=["z_coord";Z0_vect]#Z0_vect=["z_coord";Z0_vect]
			#x1=["temp";x]#x1=["temp";x;M[size(x)[1]+1:end,7]]
			#writedlm("results/paraview_results_Case2.csv",  [M2 Z0_vect x1], ',')#writedlm("results/paraview_results.csv",  [M2 Z0_vect x1], ',')
        	break
        end
        p=r+(rsnew/rsold)*p
        rsold=rsnew
		println(i)
		end
		x=x
		i=i+1
		if (i==(iter-1))&(sqrt(rsnew)<0.1)
			#x=x+rand(-5:5,size(Km)[1],1)
			i=1
			j=j+1
		end
		if i==iter#9999#size(b)[1]
			println("Convergence not achieved! Starting again!")
			i=1
			x=ones(size(Km)[1],1)
			b1=bcf(Km,M,const_bc,BC)
    		b=b1+b2
    		r=b-mult(Km,x)
    		p=r
    		rsold=(r'*r)[1]
			j=j+1
			#M=copy(M)
			#M=[["x_coord" "y_coord"];M[1:(size(x)[1]),3:4]]
			#Z0_vect=zeros((size(M)[1]-1))
			#Z0_vect=["z_coord";Z0_vect]
			#x1=["temp";x]
			#writedlm("results/paraview_results.csv",  [M Z0_vect x1], ',')
		end
		if j>20
			println("ERROR! Convergence was not achieved after $(j-1) iterations of $iter trials!")
			x=Array{Any}(undef,size(Km)[1],1)
			break
		end
    end
    x=x
	#writedlm("results/x.csv",  x, ',')
end
####End of function modconjgrad(Km,fm,x=0,in_con=0)####
function bcf(K,M)
	f=zeros(size(K)[1],1)
	pm=zeros(size(K)[1],5)

	for i=1:size(f)[1]
		#if floor(Int64,K[i,6][1])!=0
			pm[i,1]=i
			pm[i,2]=M[floor(Int64,K[i,6][1]),7]
			pm[i,3]=K[i,6][2]
		#end
		#if floor(Int64,K[i,7][1])!=0
			#pm[i,1]=i
			pm[i,4]=M[floor(Int64,K[i,7][1]),7]
			pm[i,5]=K[i,7][2]
		#end
	end
	for i in pm[:,1]
		#if i>0
			qq=floor(Int64,i)
			f[qq,1]=pm[qq,2]*pm[qq,3]+pm[qq,4]*pm[qq,5]
		#end
	end
	f=f
end
####End of function bcf(fm)####
function mult2(K,x)
	xres=zeros(size(K)[1],1)
	for i=1:size(xres)[1]
			if floor(Int64,K[i,1][1])!=0
				xres[i,1]=x[floor(Int64,K[i,1][1]),1]*K[i,1][2]
			end
			if floor(Int64,K[i,2][1])!=0
				xres[i,1]=xres[i,1]+x[floor(Int64,K[i,2][1]),1]*K[i,2][2]
			end
			if floor(Int64,K[i,3][1])!=0
				xres[i,1]=xres[i,1]+x[floor(Int64,K[i,3][1]),1]*K[i,3][2]
			end
			if floor(Int64,K[i,4][1])!=0
				xres[i,1]=xres[i,1]+x[floor(Int64,K[i,4][1]),1]*K[i,4][2]
			end
			if floor(Int64,K[i,5][1])!=0
				xres[i,1]=xres[i,1]+x[floor(Int64,K[i,5][1]),1]*K[i,5][2]
			end
	end
	xres=xres
end
####End of function mult(Km,p)####
###For defined vector b### !!!for max 5 interrelations!!! -->
function modconjgrad2(K,f,x=0,tol=1e-10)#M,x=0,in_con=0,const_bc=1,BC=0)
    if x==0
        x1=zeros(size(K)[1],1)
    end
    #b1=bcf(Km,M,const_bc,BC)
    #b=in_con+b1
    r=f-K*x1
	#println("r0=",r)
    p=copy(r)
    rsold=(r'*r)[1]
	@time begin
    for i=1:10000#size(b)[1]
		@time begin
		Ap=K*p
        alpha=rsold/(p'*Ap)[1]
		#println("Alpha=",alpha)
        x1=x1+alpha*p
		#println("x1=",x)
        r=r-alpha*(K*p)
		#println("r$i=",r)
        rsnew=(r'*r)[1]
		if sqrt(rsnew)<tol
			println("Number of iterations: $i")
            break
        end
        p=r+(rsnew/rsold)*p
		#println("p$i=",p)
        rsold=rsnew
		#println(i)
		if i==10000
			println("Convergence not achieved!")
		end
		end
	end
    end
    x1=x1
	#writedlm("results/x.csv",  x, ',')
end
####End of function modconjgrad(Km,fm,x=0,in_con=0)####
function mult(K,x)
	xres=zeros(size(K)[1],1)
	ind=size(xres)[1]
	KK=K[:,1:5]
	Threads.@threads for i=1:ind
		xres[i]=xres[i]+x[floor(Int64,K[i,1][1]),1]*K[i,1][2]
		xres[i]=xres[i]+x[floor(Int64,K[i,2][1]),1]*K[i,2][2]
		xres[i]=xres[i]+x[floor(Int64,K[i,3][1]),1]*K[i,3][2]
		xres[i]=xres[i]+x[floor(Int64,K[i,4][1]),1]*K[i,4][2]
		xres[i]=xres[i]+x[floor(Int64,K[i,5][1]),1]*K[i,5][2]
		# for j=1:5
		# 	xres[i]=xres[i]+x[floor(Int64,K[i,j][1]),1]*K[i,j][2]
		# end
	end
	xres=xres
end
####End of function mult(Km,p)####
function eqpla(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x=9999,y=9999)
	zz=[z1,z2,z3,z4]
	p=0
	if x4>998
		p1=Vector([x1,y1,z1])
		p2=Vector([x2,y2,z2])
		p3=Vector([x3,y3,z3])
		v1=p3-p1
		v2=p2-p1
		cp=cross(v1,v2)
		a=cp[1]
		b=cp[2]
		c=cp[3]
		d=dot(cp,p3)
		if (x!=9999)&(y!=9999)
			zx=-(-a*x-b*y+d)/c
			#println("zx = $zx")
		end
		zx=zx
	end
	if 11111 in zz
		k=0
		for i=1:4
			if (zz[i]!=11111)&(zz[i]!=999)
				#println(zz)
				p=p+zz[i]
				k=k+1
			end
		end
		if (x!=9999)&(y!=9999)
			zx=p/k
			#println("zx = $zx")
		end
		zx=zx
	else
		p11=Vector([x1,y1,z1])
		p21=Vector([x2,y2,z2])
		p31=Vector([x3,y3,z3])
		v11=p31-p11
		v21=p21-p11
		cp1=cross(v11,v21)
		a1=cp1[1]
		b1=cp1[2]
		c1=cp1[3]
		d1=dot(cp1,p31)
		p12=Vector([x2,y2,z2])
		p22=Vector([x3,y3,z3])
		p32=Vector([x4,y4,z4])
		v12=p32-p12
		v22=p22-p12
		cp2=cross(v12,v22)
		a2=cp2[1]
		b2=cp2[2]
		c2=cp2[3]
		d2=dot(cp2,p32)
		p13=Vector([x3,y3,z3])
		p23=Vector([x4,y4,z4])
		p33=Vector([x1,y1,z1])
		v13=p33-p13
		v23=p23-p13
		cp3=cross(v13,v23)
		a3=cp3[1]
		b3=cp3[2]
		c3=cp3[3]
		d3=dot(cp3,p33)
		p14=Vector([x4,y4,z4])
		p24=Vector([x1,y1,z1])
		p34=Vector([x2,y2,z2])
		v14=p34-p14
		v24=p24-p14
		cp4=cross(v14,v24)
		a4=cp4[1]
		b4=cp4[2]
		c4=cp4[3]
		d4=dot(cp4,p34)
		a=0.25*(abs(a1)+abs(a2)+abs(a3)+abs(a4))
		b=0.25*(abs(b1)+abs(b2)+abs(b3)+abs(b4))
		c=0.25*(abs(c1)+abs(c2)+abs(c3)+abs(c4))
		d=0.25*(abs(d1)+abs(d2)+abs(d3)+abs(d4))
		#println("equation of plane is: $a x + $b y + $c z  $d = 0 ")
		if (x!=9999)&(y!=9999)
			zx=-(-a*x-b*y+d)/c
			#println("zx = $zx")
		end
		zx=zx
	end
end
####End of function eqpla(x1,y1,z1,x2,y2,z2,x3,y3,z3)####
function res(X,M,Con,dl,xp,yp)
	eps=1e-5
	temp=[9999.]
	j=1
	x=[999. 999. 999. 999.]
	y=[999. 999. 999. 999.]
	z=[999. 999. 999. 999.]
	if (xp>maximum(M[:,3]))|(xp<minimum(M[:,3]))|(yp>maximum(M[:,4]))|(yp<minimum(M[:,4]))
		println("ERROR! Check if the coordinates ($xp, $yp) are inside the domain!")
		return
	end
	for i=1:size(X)[1]
		if (abs(xp-M[i,3])<eps)&(abs(yp-M[i,4])<eps)
			temp[1]=X[i,1]
			zx=temp[1]
			println("Temperature of coordinate ($xp, $yp) is $zx")
			break
		end
		if ((sqrt((xp-M[i,3])^2+(yp-M[i,4])^2))<(dl*sqrt(2)))&((abs(xp-M[i,3])<dl)&(abs(yp-M[i,4])<dl))
			if i<=size(Con)[1]
				if (M[Con[i,1],8]==-1)|(M[Con[i,2],8]==-1)|(M[Con[i,3],8]==-1)|(M[Con[i,4],8]==-1)|(M[Con[i,5],8]==-1)
					x[j]=M[i,3]
					y[j]=M[i,4]
					z[j]=11111. #Key for adiabatic BC
					#println(z) ##
					j=j+1
					# temp[1]=eqpla(x[1],y[1],z[1],x[2],y[2],z[2],x[3],y[3],z[3],x[4],y[4],z[4],xp,yp)
					# zx=temp[1]
					# println("Temperature of coordinate ($xp, $yp) is $zx")
				else #If coordinate is not at adiabatic BC
					#println(M[i,:])
					x[j]=M[i,3]
					y[j]=M[i,4]
					z[j]=X[i,1]
					j=j+1
				end
			else
				#println(M[i,:])
				if M[i,8]==-1
					bb=1
					while bb<=size(Con)[1]
						if i==Con[bb,2]
							z[j]=X[bb,1]
							j=j+1
							#println(z)
						end
						bb=bb+1
					end
					bb=1
					while bb<=size(Con)[1]
						if i==Con[bb,3]
							z[j]=X[bb,1]
							j=j+1
							#println(z)
						end
						bb=bb+1
					end
					bb=1
					while bb<=size(Con)[1]
						if i==Con[bb,4]
							z[j]=X[bb,1]
							j=j+1
							#println(z)
						end
						bb=bb+1
					end
					bb=1
					while bb<=size(Con)[1]
						if i==Con[bb,5]
							z[j]=X[bb,1]
							j=j+1
							#println(z)
						end
						bb=bb+1
					end
				else
				x[j]=M[i,3]
				y[j]=M[i,4]
				z[j]=X[i,1]
				#println(M[i,1])
				#println(z)
				j=j+1
				end
			end
		end
	end
	if temp[1]>9998
		#println(z)
		temp[1]=eqpla(x[1],y[1],z[1],x[2],y[2],z[2],x[3],y[3],z[3],x[4],y[4],z[4],xp,yp)
		zx=temp[1]
		#println("Temperature of coordinate ($xp, $yp) is $zx")
	end
	zx=zx
end
####End of function results(X,xp,yp)####
function init(in_con,M,dl,dt)
	incon=zeros(size(in_con)[1])
	for i=1:size(incon)[1]
		incon[i]=in_con[i]*dl^2*M[i,7]*M[i,8]/dt
	end
	incon=Vector(incon)
end
####End of function init(in_con,M,dl,dt)
function Kspar_ss(M,C,dl) 	#Conductivity matrix
	ind=Con[end,1]
	I=ones(Int64, 5*ind)
	J=ones(Int64, 5*ind)
	V=zeros(Float64, 5*ind)
	f=zeros(Float64, ind)
	j=0
		for i=1:ind
			if M[C[i,2],8]<0
				Kl=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kl=1/(1/(2*M[C[i,2],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,3],8]<0
				Kr=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kr=1/(1/(2*M[C[i,3],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,4],8]<0
				Kd=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kd=1/(1/(2*M[C[i,4],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,5],8]<0
				Ku=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Ku=1/(1/(2*M[C[i,5],6])+1/(2*M[C[i,1],6]))
			end
			I[i+j]=i
			J[i+j]=Con[i,1]
			V[i+j]=Kl+Kr+Kd+Ku
			if C[i,2]<=ind
				j=j+1
				I[i+j]=i
				J[i+j]=Con[i,2]
				V[i+j]=-Kl
			else
				f[i]=f[i]+Kl*M[C[i,2],7]
			end
			if C[i,3]<=ind
				j=j+1
				I[i+j]=i
				J[i+j]=Con[i,3]
				V[i+j]=-Kr
			else
				f[i]=f[i]+Kr*M[C[i,3],7]
			end
			if C[i,4]<=ind
				j=j+1
				I[i+j]=i
				J[i+j]=Con[i,4]
				V[i+j]=-Kd
			else
				f[i]=f[i]+Kd*M[C[i,4],7]
			end
			if C[i,5]<=ind
				j=j+1
				I[i+j]=i
				J[i+j]=Con[i,5]
				V[i+j]=-Ku
			else
				f[i]=f[i]+Ku*M[C[i,5],7]
			end
		end
		# println(I)
		# println(J)
		# println(V)
		# println(find)
		# println(fval)
		K=sparse(I,J,V)
		f=sparsevec(f)
		Sind=[K,f]

	#writedlm("mesh/K.csv",  K, ',')
	#writedlm("mesh/f.csv",  f, ',')
	Sind=Sind
end
####End of function Kspar_ss(M,C,dl)
function Kspar_t(M,C,dl,dt) 	#Conductivity matrix
	ind=Con[end,1]
	I=ones(Int64, 5*ind)
	J=ones(Int64, 5*ind)
	V=zeros(Float64, 5*ind)
	f=zeros(Float64, ind)
	j=0
		for i=1:ind
			if M[C[i,2],8]<0
				Kl=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kl=1/(1/(2*M[C[i,2],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,3],8]<0
				Kr=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kr=1/(1/(2*M[C[i,3],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,4],8]<0
				Kd=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Kd=1/(1/(2*M[C[i,4],6])+1/(2*M[C[i,1],6]))
			end
			if M[C[i,5],8]<0
				Ku=1e-10#dl/2*(1/(M[C[i,1],6]))
			else
				Ku=1/(1/(2*M[C[i,5],6])+1/(2*M[C[i,1],6]))
			end
			a=M[C[i,1],7]*M[C[i,1],8]
			I[i+j]=i
			J[i+j]=Con[i,1]
			V[i+j]=a*dl^2/dt+Kl+Kr+Kd+Ku
			if C[i,2]<=ind
				j=j+1
				I[i+j]=i
				J[i+j]=Con[i,2]
				V[i+j]=-Kl
			else
				f[i]=f[i]+Kl*M[C[i,2],7]
			end
			if C[i,3]<=ind
				j=j+1
				I[i+j]=i
				J[i+j]=Con[i,3]
				V[i+j]=-Kr
			else
				f[i]=f[i]+Kr*M[C[i,3],7]
			end
			if C[i,4]<=ind
				j=j+1
				I[i+j]=i
				J[i+j]=Con[i,4]
				V[i+j]=-Kd
			else
				f[i]=f[i]+Kd*M[C[i,4],7]
			end
			if C[i,5]<=ind
				j=j+1
				I[i+j]=i
				J[i+j]=Con[i,5]
				V[i+j]=-Ku
			else
				f[i]=f[i]+Ku*M[C[i,5],7]
			end
		end
		# println(I)
		# println(J)
		# println(V)
		# println(find)
		# println(fval)
		K=sparse(I,J,V)
		f=sparsevec(f)
		Sind=[K,f]

	#writedlm("mesh/K.csv",  K, ',')
	#writedlm("mesh/f.csv",  f, ',')
	Sind=Sind
end
####End of function Kspar_t(M,C,dl,dt)

function moist(M,x_t,x_m)
	δ=1.85e-10
	Mm=copy(M)
	for i=1:size(x_t)[1]
		if x_t[i]<0
			a=22.544
			θ=272.44
			psat=611*exp((a*x_t[i])/(θ-x_t[i]))
			μm=μ_m(x_m[i],M[i,5])
			Mm[i,6]=psat*δ/μm
			Mm[i,7]=dw_dϕ(x_m[i],M[i,5])
			Mm[i,8]=1
		else
			a=17.08
			θ=234.18
			psat=611*exp((a*x_t[i])/(θ-x_t[i]))
			μm=μ_m(x_m[i],M[i,5])
			Mm[i,6]=psat*δ/μm
			Mm[i,7]=dw_dϕ(x_m[i],M[i,5])
			Mm[i,8]=1
		end
	end
	for i=(size(x_t)[1]+1):floor(Int64,M[end,1])
		Mm[i,6]=1e-3
	end
	Mm=Mm
end
####End of function moist(M,x_t,x_m)


text:
numero = 8 * 9 / 2       #= Conta basica =#
outro_numero = (9 + 2) * 3

println(numero + outro_numero)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using IntervalSets

@testset "util" begin
    @testset "BAT.fromui" begin
        @test BAT._all_in_ui(ones(2, 2))
        @test !BAT._all_in_ui([0.4, 0.2, 1.1])

        @test BAT.fromui(0.5, -1.0, 1.0) ≈ 0.0
        @test BAT.fromui(0.5, -0.5, 5.0) ≈ 2.25
        @test BAT.fromui(0.1, -0.5, 1.0) ≈ -0.35
        @test_throws ArgumentError BAT.fromui(-0.1, -1.0, 1.0)
        @test_throws ArgumentError BAT.fromui(1.1, -0.5, 1.0)

        @test BAT.fromui(0.5, ClosedInterval(-0.5, 1.0)) ≈ 0.25

        inv_fromui = @inferred inv(BAT.fromui)
        @test inv_fromui  == BAT.inv_fromui
        @test inv(inv_fromui) == BAT.fromui

        @test inv_fromui(2.25, -0.5, 5.0) ≈ 0.5

        @test inv_fromui(-0.35, -0.5, 1.0) ≈ 0.1
        @test inv_fromui(5.0, -0.5, 5.0) ≈ 1.0
        @test inv_fromui(-0.5, -0.5, 5.0) ≈ 0.0
        @test_throws ArgumentError inv_fromui(-0.2, -0.1, 3)
        @test_throws ArgumentError inv_fromui(3.1, -0.1, 3)

        @test inv_fromui(-0.35, ClosedInterval(-0.5, 1.0)) ≈ 0.1
    end
end


text:
const FILE_PATH_PROPS = ["polygon_file", "source_file", "ground_file", "mask_file", "output_file", "habitat_file", "point_file", "reclass_file"]
const DEFAULT_CFG     = {
        "Version" => {
            "version" => "unknown"
        },
        "Connection scheme for raster habitat data" => {
            "connect_four_neighbors_only" => false,
            "connect_using_avg_resistances" => false,
        },
        "Short circuit regions (aka polygons)" => {
            "use_polygons" => false,
            "polygon_file" => "(Browse for a short-circuit region file)"
        },
        "Options for advanced mode" => {
            "source_file" => "(Browse for a current source file)",
            "ground_file" => "(Browse for a ground point file)",
            "ground_file_is_resistances" => true,
            "use_unit_currents" => false,
            "use_direct_grounds" => false,
            "remove_src_or_gnd" => "keepall"
        },
        "Mask file" => {
            "mask_file" => nothing,
            "use_mask" => false
        },
        "Calculation options" => {
            "preemptive_memory_release" => false,
            "low_memory_mode" => false,
            "parallelize" => false,          # can parallelize if true. It may be overridden internally to false based on os support and suitability of the problem.
            "max_parallel" => 0,            # passing 0 results in using all available cpus
            "print_timings" => false,
            "print_rusages" => false,
            "solver" => "cg+amg"
        },
        "Options for one-to-all and all-to-one modes" => {
            "use_variable_source_strengths" => false,
            "variable_source_file" => nothing
        },
        "Output options" => {
            "set_null_currents_to_nodata" => false,
            "output_file" => "(Choose a base name for output files)",
            "write_cum_cur_map_only" => false,
            "log_transform_maps" => false,
            "write_max_cur_maps" => false,
            "compress_grids" => false,
            "set_null_voltages_to_nodata" => false,
            "set_focal_node_currents_to_zero" => false,
            "write_volt_maps" => false,
            "write_cur_maps" => false
        },
        "Habitat raster or graph" => {
            "habitat_map_is_resistances" => true,
            "habitat_file" => "(Browse for a resistance file)"
        },
        "Circuitscape mode" => {
            "scenario" => "not entered",
            "data_type" => "raster"
        },
        "Options for pairwise and one-to-all and all-to-one modes" => {
            "use_included_pairs" => false,
            "included_pairs_file" => "(Browse for a file with pairs to include or exclude)",
            "point_file" => "(Browse for file with locations of focal points or regions)"
        },
        "Options for reclassification of habitat data" => {
            "use_reclass_table" => false,
            "reclass_file" => "(Browse for file with reclassification data)"
        },
        "Logging Options" => {
            "profiler_log_file" => nothing,      # file to log timing and rusage profiling results
            "log_file" => nothing,               # file to log regular log messages
            "log_level" => "INFO",           # one of FATAL, ERROR, WARN, INFO, DEBUG
            "screenprint_log" => false        # whether to print logs to console (stdout)
        }
    }

const CHECKS_AND_MESSAGES = {
        "scenario" => "Please choose a scenario",
        "habitat_file" => "Please choose a resistance file",
        "output_file" => "Please choose an output file name",
        "point_file" => "Please choose a focal node file",
        "source_file" => "Please enter a current source file",
        "ground_file" => "Ground point file does not exist!",
        "reclass_file" => "Please choose a file with reclassification data",
        "polygon_file" => "Please enter a short-circuit region file or uncheck this option in the Options menu"
    }

type CSCfg
    inifile::Inifile
    val_dict::IniFile.HTSS
    logger::Logger
    filename::String

    function CSCfg(cfgfile::String="", rel_to_abs::String="")
        ret = new()
        ret.logger = getlogger("circuitscape.cfg")
        debug(ret.logger, "created CSCfg with config file \"$cfgfile\"")
        defaults = IniFile.HTSS()
        for section_vals in values(DEFAULT_CFG)
            merge!(defaults, section_vals)
        end

        inifile = Inifile()
        !isempty(cfgfile) && read(inifile, cfgfile)
        inifile.defaults = defaults

        ret.inifile = inifile
        ret.val_dict = as_dict(ret, rel_to_abs)
        ret.filename = cfgfile
        ret
    end
end

function show(io::IO, cfg::CSCfg)
    println(io, "CSCfg with $(length(cfg.val_dict)) entries")
    if !isempty(cfg.filename)
        println(io, "file: $(cfg.filename)")
    end
end

approp(val) = val
function approp(val::String)
    ((val == "nothing") || (val == "None")) && (return nothing)
    try
        val = parsebool(val)
    catch
        try
            val = parseint(val)
        catch
            try
                val = parsefloat(val)
            catch
                # ignore
            end
        end
    end
    val
end

function parsebool(str::String)
    str = lowercase(str)
    (str == "true") && return true
    (str == "false") && return false
    error("invalid boolean representation")
end

function as_dict(cscfg::CSCfg, rel_to_abs::String)
    ret = IniFile.HTSS()
    for (section_name, htss) in DEFAULT_CFG
        for (key, val) in htss
            cfgval = get(cscfg.inifile, section_name, key, val)
            if key in FILE_PATH_PROPS
                if (cfgval == nothing) || (val == cfgval)
                    cfgval = ""
                elseif !isempty(rel_to_abs) && !isabspath(cfgval)
                    cfgval = joinpath(rel_to_abs, cfgval)
                end
            end
            ret[key] = approp(cfgval)
        end
    end
    ret
end

function write(io::IO, cscfg::CSCfg, pretty_print::Bool=true)
    write_space = false
    for (section_name, htss) in DEFAULT_CFG
        if pretty_print
            if write_space
                println(io, "")
            else
                write_space = true
            end
        end
        println(io, "[$section_name]")

        for (key, val) in htss
            cfgval = get(cscfg, key)
            println(io, "$key=$cfgval")
        end
    end
end

get(cscfg::CSCfg, key::String) = cscfg.val_dict[key]
set(cscfg::CSCfg, key::String, val::IniFile.INIVAL) = (cscfg.val_dict[key] = approp(val))
function get(cscfg::CSCfg, key::String, typ::Type)
    v = get(cscfg, key)
    vtyp = typeof(v)
    if isa(v, String) && !issubtype(vtyp, typ)
        lcv = lowercase(v)
        ((lcv == "none") || (lcv == "nothing")) && (return nothing)
        if (typ == Bool)
            (v == "false") && (return false)
            (v == "true") && (return true)
            error("Invalid boolean value $v")
        else
            return convert(typ, v)
        end
    end
    return v
end