text:
# Tests for validated_integ

using TaylorModels
# using LinearAlgebra: norm
using Test
# using Random

const _num_tests = 1_000

setformat(:full)

# NOTE: IntervalArithmetic v0.16.0 includes this function; but
# IntervalRootFinding is bounded to use v0.15.x
interval_rand(X::Interval{T}) where {T} = X.lo + rand(T) * (X.hi - X.lo)
interval_rand(X::IntervalBox) = interval_rand.(X)

function test_integ(fexact, t0, qTM, q0, δq0)
    normalized_box = symmetric_box(length(q0), Float64)
    # Time domain
    domt = domain(qTM[1])
    # Random time (within time domain) and random initial condition
    δt = rand(domt)
    δtI = (δt .. δt) ∩ domt
    q0ξ = interval_rand(δq0)
    q0ξB = IntervalBox([(q0ξ[i] .. q0ξ[i]) ∩ δq0[i] for i in eachindex(q0ξ)])
    # Box computed to overapproximate the solution at time δt
    q = evaluate.(evaluate.(qTM, δtI), (normalized_box,))
    # Box computed from the exact solution must be within q
    bb = all(fexact(t0+δtI, q0 .+ q0ξB) .⊆ q)
    # Display details if bb is false
    bb || @show(t0, domt, remainder.(qTM),
            δt, δtI, q0ξ, q0ξB, q,
            fexact(t0+δtI, q0 .+ q0ξB))
    return bb
end

@testset "Tests for `validated_integ`" begin
    @testset "falling_ball!" begin
        @taylorize function falling_ball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, t0, x0) = (x0[1] + x0[2]*(t-t0) - 0.5*(t-t0)^2, x0[2] - (t-t0))

        # Initial conditions
        tini, tend = 0.0, 10.0
        normalized_box = symmetric_box(2, Float64)
        q0 = [10.0, 0.0]
        δq0 = 0.25 * normalized_box
        X0 = IntervalBox(q0 .+ δq0)

        # Parameters
        abstol = 1e-20
        orderQ = 2
        orderT = 4
        ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))

        @testset "Forward integration 1" begin
            tTM, qv, qTM = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)

            @test length(qv) == length(qTM[1, :]) == length(tTM)

            end_idx = lastindex(tTM)
            # Random.seed!(1)
            for it = 1:_num_tests
                n = rand(2:end_idx)
                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], qTM[:,n], q0, δq0)
            end

            tTMf, qvf, qTMf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
                adaptive=false)
            @test length(qvf) == length(qv)
            @test qTM == qTMf

            # initializaton with a Taylor model
            X0tm = qTM[:, 1]
            tTM2, qv2, qTM2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
            @test qTM == qTM2
        end

        @testset "Forward integration 2" begin
            tTM, qv, qTM = validated_integ2(falling_ball!, X0,
                tini, tend, orderQ, orderT, abstol)

            @test length(qv) == length(qTM[1, :]) == length(tTM)

            # Random.seed!(1)
            end_idx = lastindex(tTM)
            for it = 1:_num_tests
                n = rand(2:end_idx)
                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], qTM[:,n], q0, δq0)
            end

            # initializaton with a Taylor model
            X0tm = qTM[:, 1]
            tTM2, qv2, qTM2 = validated_integ2(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
            @test qTM == qTM2
        end

        # Initial conditions
        tini, tend = 10.0, 0.0
        q0 = [10.0, 0.0]
        δq0 = IntervalBox(-0.25 .. 0.25, 2)
        X0 = IntervalBox(q0 .+ δq0)

        @testset "Backward integration 1" begin
            tTM, qv, qTM = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)

            @test length(qv) == length(qTM[1, :]) == length(tTM)

            # Random.seed!(1)
            end_idx = lastindex(tTM)
            for it = 1:_num_tests
                n = rand(2:end_idx)
                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], qTM[:,n], q0, δq0)
            end

            tTMf, qvf, qTMf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
                adaptive=false)
            @test length(qvf) == length(qv)
            @test all(qTM .== qTMf)

            # initializaton with a Taylor model
            X0tm = qTM[:, 1]
            tTM2, qv2, qTM2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
            @test qTM == qTM2

            tTM2f, qv2f, qTM2f = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol,
                adaptive=false)
            @test length(qv2f) == length(qv2)
            @test all(qTM .== qTM2f)
        end

        @testset "Backward integration 2" begin
            tTM, qv, qTM = validated_integ2(falling_ball!, X0,
            tini, tend, orderQ, orderT, abstol)

            @test length(qv) == length(qTM[1, :]) == length(tTM)

            # Random.seed!(1)
            end_idx = lastindex(tTM)
            for it = 1:_num_tests
                n = rand(2:end_idx)
                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], qTM[:,n], q0, δq0)
            end
        end
    end

    @testset "x_square!" begin
        @taylorize function x_square!(dx, x, p, t)
            dx[1] = x[1]^2
            nothing
        end

        exactsol(t, x0) = 1 / (1/x0[1] - t)

        tini, tend = 0., 0.45
        normalized_box = symmetric_box(1, Float64)
        abstol = 1e-15
        orderQ = 5
        orderT = 20
        q0 = [2.]
        δq0 = 0.0625 * normalized_box
        X0 = IntervalBox(q0 .+ δq0)
        ξ = set_variables("ξₓ", numvars=1, order=2*orderQ)

        @testset "Forward integration 1" begin
            tTM, qv, qTM = validated_integ(x_square!, X0, tini, tend, orderQ, orderT, abstol)

            @test length(qv) == length(qTM[1, :]) == length(tTM)

            # Random.seed!(1)
            end_idx = lastindex(tTM)
            for it = 1:_num_tests
                n = rand(1:end_idx)
                @test test_integ((t,x)->exactsol(t,x), tTM[n], qTM[:,n], q0, δq0)
            end

            tTMf, qvf, qTMf = validated_integ(x_square!, X0, tini, tend, orderQ, orderT, abstol,
                adaptive=false)
            @test length(qvf) == length(qv)
            @test all(qTMf .== qTM)

            # initializaton with a Taylor model
            X0tm = copy(qTM[:, 1])
            tTM2, qv2, qTM2 = validated_integ(x_square!, X0tm, tini, tend, orderQ, orderT, abstol)
            @test qTM == qTM2
        end

        @testset "Forward integration 2" begin
            tTM, qv, qTM = validated_integ2(x_square!, X0, tini, tend, orderQ, orderT, abstol)

            @test length(qv) == length(qTM[1, :]) == length(tTM)

            # Random.seed!(1)
            end_idx = lastindex(tTM)
            for it = 1:_num_tests
                n = rand(1:end_idx)
                @test test_integ((t,x)->exactsol(t,x), tTM[n], qTM[:,n], q0, δq0)
            end
        end
    end

    @testset "Pendulum with constant torque" begin
        @taylorize function pendulum!(dx, x, p, t)
            si = sin(x[1])
            aux = 2 *  si
            dx[1] = x[2]
            dx[2] = aux + 8*x[3]
            dx[3] = zero(x[1])
            nothing
        end
        # Conserved quantity
        ene_pendulum(x) = x[2]^2/2 + 2 * cos(x[1]) - 8 * x[3]

        # Initial conditions
        tini, tend = 0.0, 12.0
        q0 = [1.1, 0.1, 0.0]
        δq0 = IntervalBox(-0.1 .. 0.1, -0.1 .. 0.1, 0..0)
        X0 = IntervalBox(q0 .+ δq0)
        ene0 = ene_pendulum(X0)

        # Parameters
        abstol = 1e-10
        orderQ = 3
        orderT = 10
        ξ = set_variables("ξ", order=2*orderQ, numvars=length(q0))

        tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol);
        @test all(ene0 .⊆ ene_pendulum.(qv))

        tTM, qv, qTM = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
            validatesteps=32);
        @test all(ene0 .⊆ ene_pendulum.(qv))

        # Initial conditions 2
        q0 = [1.1, 0.1, 0.0]
        δq0 = IntervalBox(-0.1 .. 0.1, -0.1 .. 0.1, -0.01 .. 0.01)
        X0 = IntervalBox(q0 .+ δq0)
        ene0 = ene_pendulum(X0)

        tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol);
        @test all(ene0 .⊆ ene_pendulum.(qv))

        tTM, qv, qTM = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
            validatesteps=32);
        @test all(ene0 .⊆ ene_pendulum.(qv))
    end
end


text:
module CompactBases

import Base: eltype, axes, size, ==, getindex, checkbounds, copyto!, similar, show, step
import Base.Broadcast: materialize

using ContinuumArrays
import ContinuumArrays: Basis, ℵ₁, Derivative, Inclusion, @simplify, simplifiable, _simplify,
    AbstractQuasiArrayApplyStyle

using QuasiArrays
import QuasiArrays: AbstractQuasiMatrix, QuasiAdjoint, MulQuasiArray,
    PInvQuasiMatrix, InvQuasiMatrix, QuasiDiagonal,
    BroadcastQuasiArray, SubQuasiArray, LazyQuasiMatrix

using IntervalSets

using LazyArrays

using LinearAlgebra
import LinearAlgebra: Matrix, dot
using SparseArrays
using FillArrays

using OffsetArrays
using BandedMatrices
using BlockBandedMatrices

using FastGaussQuadrature

using Formatting
using RecipesBase

vandermonde(B::Basis) = B[locs(B),:]

include("distributions.jl")
include("restricted_bases.jl")
include("types.jl")
include("unit_vectors.jl")
include("intervals.jl")
include("knot_sets.jl")
include("quadrature.jl")
include("skewtridiag.jl")

include("materialize_dsl.jl")

include("coulomb_derivative.jl")

include("fornberg.jl")
include("finite_differences.jl")
include("fd_operators.jl")
include("fd_derivatives.jl")

include("fedvr.jl")
include("fedvr_operators.jl")
include("fedvr_derivatives.jl")

include("bsplines.jl")
include("bspline_derivatives.jl")

include("inner_products.jl")
include("densities.jl")
include("linear_operators.jl")

"""
    centers(B)

Return the locations of the mass centers of all basis functions of
`B`; for orthogonal bases such as finite-differences and FE-DVR, this
is simply [`locs`](@ref), i.e. the location of the quadrature nodes.
"""
centers(B::BasisOrRestricted) = locs(B)

"""
    metric(B)

Returns the metric or mass matrix of the basis `B`, equivalent to
`S=B'B`.
"""
metric(B::BasisOrRestricted) = B'B

"""
    metric_shape(B)

Returns the shape of the metric or mass matrix of the basis `B`,
i.e. `Diagonal`, `BandedMatrix`, etc.
"""
metric_shape(B::BasisOrRestricted) = typeof(metric(B))

function matrix_element_metric(B::BasisOrRestricted)
    # For non-orthogonal bases, we need to apply the metric inverse,
    # after computing the action of a matrix representation of a
    # linear operator on a ket. However, when computing the inner
    # product with the bra afterwards, we again use the metric. Here,
    # we short-circuit this by combining them into them identity
    # operator I, if the operator metric equals the metric. For
    # orthogonal bases, with the integration weights not built into the
    # coefficients, and an identity matrix operator metric, we /do/
    # need to use the metric for the inner product.
    S = metric(B)
    oS = operator_metric(B)
    oS == S ? I : (oS \ S)
end

export AbstractFiniteDifferences,
    FiniteDifferences, StaggeredFiniteDifferences, ImplicitFiniteDifferences,
    Derivative, dot, QuasiDiagonal, Inclusion, .., distribution,
    FEDVR, Derivative, @elem, dot,
    BSpline,
    ArbitraryKnotSet, LinearKnotSet, ExpKnotSet,
    order, numintervals, numfunctions, nonempty_intervals,
    centers, vandermonde,
    metric, metric_shape, matrix_element_metric,
    FunctionProduct, Density,
    LinearOperator, DiagonalOperator, ShiftAndInvert

end # module

ext: jl


text:
"""
Main module for `GigaSOM.jl` - Huge-scale, high-performance flow cytometry clustering

The documentation is here: http://LCSB-BioCore.github.io/GigaSOM.jl
"""

module GigaSOM

using CSV
using DataFrames
using Distances
using Distributed
using Distributions
using FCSFiles
using FileIO
using DistributedArrays
using NearestNeighbors
using Serialization
using StableRNGs

include("base/structs.jl")

include("base/dataops.jl")
include("base/distributed.jl")
include("base/trainutils.jl")

include("analysis/core.jl")
include("analysis/embedding.jl")

include("io/dio.jl")
include("io/input.jl")
include("io/process.jl")
include("io/splitting.jl")

# include visualization files
# include("visualization/plotting.jl")

#core
export initGigaSOM, trainGigaSOM, mapToGigaSOM

#trainutils
export linearRadius, expRadius, gaussianKernel, bubbleKernel, thresholdKernel, distMatrix

#embedding
export embedGigaSOM

# structs
export Som, LoadedDataInfo

#io/input
export readFlowset,
    readFlowFrame,
    loadFCS,
    loadFCSHeader,
    getFCSSize,
    loadFCSSizes,
    loadFCSSet,
    selectFCSColumns,
    distributeFCSFileVector,
    distributeFileVector,
    getCSVSize,
    loadCSV,
    loadCSVSizes,
    loadCSVSet

#io/splitting
export slicesof, vcollectSlice, collectSlice

#io/process
export cleanNames!, getMetaData, getMarkerNames

# plotting
export plotCounts, plotPCA

#dataops (higher-level operations on data)
export dcopy,
    dselect,
    dapply_cols,
    dapply_rows,
    dstat,
    dstat_buckets,
    dcount,
    dcount_buckets,
    dscale,
    dtransform_asinh,
    dmedian,
    dmedian_buckets,
    mapbuckets,
    catmapbuckets

#distributed data tools
export save_at,
    get_from,
    get_val_from,
    remove_from,
    distribute_array,
    distribute_darray,
    undistribute,
    distributed_transform,
    distributed_mapreduce,
    distributed_foreach,
    distributed_collect,
    distributed_export,
    distributed_import,
    distributed_unlink

end # module


text:
# Run some basic tests to check that all the functions are correct
using Base.Test, ProgressMeter
using AutomotiveDrivingModels, POMDPs, GridInterpolations, StatsBase, POMDPToolbox

# include all the relevant files in the right order
include("occluded_intersection_env.jl");
include("pomdp_types.jl")
include("spaces.jl")
include("transition.jl")
include("observation.jl")

#### HELPER FOR TESTING #####

# check that state indexing is consistent
function check_indexing(pomdp::SingleOIPOMDP, state_space::Vector{SingleOIState})
    for (i,s) in enumerate(state_space)
        @test stateindex(pomdp, s) == i
    end
end

# custom type for ADM
mutable struct PomdpAction
    a::SingleOIAction
    pomdp::SingleOIPOMDP
end
# custom propagate for testing
function AutomotiveDrivingModels.propagate(veh::Vehicle, a::PomdpAction, roadway::Roadway, ΔT::Float64)
    states, probs = ego_transition(pomdp, veh.state, a.a, ΔT)
    return sample(states, Weights(probs))
end
# custom model for testing
mutable struct EgoTestModel <: DriverModel{PomdpAction}
    a::PomdpAction
end
Base.rand(model::EgoTestModel) = model.a

# test 0 acceleration
function test_no_motion(pomdp::SingleOIPOMDP)
    ego = Vehicle(ego_state(pomdp, pomdp.s_start, 0.), VehicleDef(), 1)
    scene = Scene()
    push!(scene, ego)
    models = Dict{Int, DriverModel}()
    models[1] = EgoTestModel(PomdpAction(SingleOIAction(0.0), pomdp))
    nticks = 100
    rec = SceneRecord(nticks+1, 0.1)
    simulate!(rec, scene, pomdp.env.roadway, models, nticks)
    @test rec[0][1].state == ego.state
    @test rec[0][1].state.posF.s == pomdp.s_start
end

# test max acceleration
function test_reaching_goal(pomdp::SingleOIPOMDP)
    ego = Vehicle(ego_state(pomdp, pomdp.s_start, 0.), VehicleDef(), 1)
    scene = Scene()
    push!(scene, ego)
    models = Dict{Int, DriverModel}()
    models[1] = EgoTestModel(PomdpAction(SingleOIAction(pomdp.max_acc), pomdp))
    nticks = 100
    rec = SceneRecord(nticks+1, 0.1)
    simulate!(rec, scene, pomdp.env.roadway, models, nticks)
    @test rec[0][1].state.posF.s == pomdp.s_goal
end

# Test that transition sums up to 1
function test_transition(pomdp::SingleOIPOMDP)
    states_to_test = rand(ordered_states(pomdp), 100)
    @showprogress 0.5 "Checking transitions for 100 random states..." for (i, s) in enumerate(states_to_test)
        for (j, a) in enumerate(iterator(actions(pomdp)))
            prob = 0.
            d = transition(pomdp, s, a)
            for (k, sp) in enumerate(ordered_states(pomdp))
                prob += pdf(d, sp)
            end
            if !(prob ≈ 1.0)
                println(" Transition sums up to : ", prob)
                println("   For state action pair : ", s, a)
                break
            end
        end
    end
end

# Test that transition sums up to 1
function test_full_transition(pomdp::SingleOIPOMDP)
    @showprogress 0.5 "Checking transitions..." for (i, s) in enumerate(ordered_states(pomdp))
        for (j, a) in enumerate(iterator(actions(pomdp)))
            prob = 0.
            d = transition(pomdp, s, a)
            for (k, sp) in enumerate(ordered_states(pomdp))
                prob += pdf(d, sp)
            end
            if !(prob ≈ 1.0)
                println(" Transition sums up to : ", prob)
                println("   For state action pair : ", s, a)
                break
            end
        end
    end
end

#Test that observations sum up to 1
function test_observation(pomdp::SingleOIPOMDP)
    @showprogress 0.5 "Checking observations..." for (i, sp) in enumerate(ordered_states(pomdp))
        for (j, a) in enumerate(iterator(actions(pomdp)))
            prob = 0.
            d = observation(pomdp, a, sp)
            for (k, o) in enumerate(ordered_observations(pomdp))
                prob += pdf(d, o)
            end
            if !(prob ≈ 1.0)
                println( " Observation sums up to : ", prob)
                println("   For state action pair : ", sp, a)
                break
            end
        end
    end
end

# fix rng for the test
rng = MersenneTwister(3);

# initialize the environment and the pomdp
println("Testing Constructors")
params = EnvParams()
env = SimpleInterEnv(params);
pomdp = SingleOIPOMDP(env = env)

#### TEST STATE SPACES
println("Testing state space")
state_space = states(pomdp)
@test !isempty(state_space)
@test n_states(pomdp) == length(state_space)
check_indexing(pomdp, state_space)
println("SingleOIPOMDP with default parameters has ", n_states(pomdp), " states")

#test helpers
@test off_the_grid(pomdp, get_off_the_grid(pomdp))
s_grid = get_ego_s_grid(pomdp)
v_grid = get_v_grid(pomdp)
@test n_car_states(pomdp) == length(car_states(pomdp))

#### TEST ACTION SPACE
println("Testing action space")
action_space = actions(pomdp)
@test !isempty(action_space)
@test n_actions(pomdp) == length(action_space)

### TEST DISTRIBUTION

test_no_motion(pomdp)
test_reaching_goal(pomdp)

#TODO Test Car Transition

test_observation(pomdp)
test_transition(pomdp)
test_full_transition(pomdp)

# Check solver requirements

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using Robotlib, TotalLeastSquares, FillArrays
# using HomotopyContinuation, DynamicPolynomials
using Test, Random
RTR(R) = R'R
vecangle(g,gf) = acos(min(1,g⋅gf/norm(gf)/norm(g)))*180/pi

rcond(K) = /(svdvals(K)[[1,3]]...)
toR(r) = toOrthoNormal(reshape(real(r),3,3))
tor(r) = toR(r)[:]

function rqi(K,r,n=1)
    for i = 1:n
        r = (K-I)\r
        r = r/norm(r)
    end
    r
end

function eigenR(K)
    v = real(eigen(K).vectors[:,1])
    R = reshape(v,3,3)
    det(R) < 0 && (R .*= -1)
    toOrthoNormal(R)
end
# using SumOfSquares, PolyJuMP, Ipopt
# m = SOSModel(with_optimizer(Ipopt.Optimizer))

# N         = 1000
# Rf        = Robotlib.toOrthoNormal(randn(3,3))
# mf        = 1
# gf        = [0,0,1] + 0.3randn(3)
# POSES     = cat([Robotlib.toOrthoNormal(randn(3,3)) for i = 1:N]..., dims=3)
# forces    = hcat([Rf'*(POSES[1:3,1:3,i]'*gf) for i = 1:N]...)' |> copy
#
#
# function calibForceHomotopy(POSES,F,m0=1)
#     N = size(POSES,3)
#     @polyvar w[1:9]
#     @polyvar g[1:3]
#     @polyvar λ[1:6]
#     I = Robotlib.I3
#     A  = Array{eltype(F)}(undef, 3N, 9)
#     B  = Array{Any}(undef, 3N)
#     At = (F,i) -> [F[i,1]*I   F[i,2]*I    F[i,3]*I]
#     for i = 1:N
#         RA               = POSES[1:3,1:3,i]'
#         A[3(i-1)+1:3i,:] = At(F,i)
#         B[3(i-1)+1:3i]   = -RA*g
#     end
#     @info("Condition number of Gram matrix: ", cond(A'A))
#     R = reshape(w,3,3)
#     J = sum((A*w-B).^2) + sum(1:3) do i
#         λ[i]*(R[:,i]'R[:,i] - 1)
#     end + λ[4]*R[:,1]'R[:,2] +
#     λ[5]*R[:,1]'R[:,3] +
#     λ[6]*R[:,2]'R[:,3]
#
#     dJ     = DynamicPolynomials.differentiate(J, [w;g;λ])
#     result = solve(dJ)
#     reals  = realsolutions(result)
#     costfun = w-> begin
#     sum((A*vec(w[1:9])-[b(w[10:12]) for b in B]).^2)
# end
# # R = tomat(reals[minindex])
# result, costfun
# end
#
#
#
# result, J = calibForceHomotopy(POSES,forces)
# reals = realsolutions(result)
# tomat(x) = reshape(x[1:9],3,3) |> copy
#
# Rs = tomat.(reals)
# RTR(R) = R'R
# minval, minindex = findmin(J.(reals))
# R = Rs[minindex]

function calibForceIterative(POSES,F,g; trace=false)
    N = size(POSES,3)
    I = Robotlib.I3
    A  = Array{eltype(F)}(undef, 3N, 3)
    B = F'[:]
    local m
    trace && (Rg = [])
    for iter = 1:6
        Rf, m = Robotlib.Calibration.calibForce(POSES,F,g; offset=false, verbose=true)
        for i = 1:N
            A[3(i-1)+1:3i,:] = Rf'POSES[:,:,i]'
        end
        g = A\B
        trace && push!(Rg, (Rf, g))
    end
    trace && (return Rf,g,m,Rg)
    Rf,g,m
end

function calibForceIterative2(POSES,F,g; trace=false)
    N = size(POSES,3)
    I = Robotlib.I3
    A  = Array{eltype(F)}(undef, 3N, 3)
    for i = 1:N
        A[3(i-1)+1:3i,:] = POSES[:,:,i]'
    end
    A = factorize(A)
    local m
    trace && (Rg = [])
    for iter = 1:6
        Rf, m = Robotlib.Calibration.calibForce(POSES,F,g; offset=false, verbose=false)
        B = (Rf*F')[:]
        g = A\B
        trace && push!(Rg, (Rf, g))
    end
    trace && (return Rf,g,m,Rg)
    Rf,g,m
end

function calibForceIterative3(POSES,F; trace=false)
    N = size(POSES,3)
    I = Robotlib.I3
    Rf = I
    A  = Array{eltype(F)}(undef, 3N, 3)
    B  = Array{eltype(F)}(undef, 3N)
    local m, g
    trace && (Rg = [])
    for iter = 1:6
        for i = 1:N
            A[3(i-1)+1:3i,:] = Rf'POSES[:,:,i]'
            B[3(i-1)+1:3i]   = F[i,:]
        end
        g = A\B
        Rf, m = Robotlib.Calibration.calibForce(POSES,F,g; offset=false, verbose=false)
        trace && push!(Rg, (Rf, g))
    end
    trace && (return Rf,g,m,Rg)
    Rf,g,m
end
using FillArrays

function calibForceIterative4(POSES,forces; trace=false)
    N = size(POSES,3)
    I3 = Matrix{Float64}(I, 3, 3)
    r = vec(I3)
    # D = permutedims(POSES, (2,1,3))
    # D = reshape(D, 3N,3)
    D = reshape(POSES, 3,3N)'
    F = kron(forces, I3)
    FF = pinv(F)*D
    DD = pinv(D)*F
    K = FF*DD
    e = eigen(K)
    sort(e.values, by=abs,rev=true)
    trace && (Rg = [])
    for iter = 1:5
        r = K*r
        r = vec(toOrthoNormal(reshape(r,3,3)))
        # trace && push!(Rg, (Rf, g))
    end
    # trace && (return Rf,g,m,Rg)
    # Rf,g,m
    toOrthoNormal(reshape(r,3,3)), e
end

function calibForceEigen(POSES,forces)
    N = size(POSES,3)
    D = reshape(POSES, 3,3N)'
    F = kron(forces, Eye(3))
    K = F\D*(D\F)
    R̂ = eigenR(K)
    ĝ = D\F*vec(R̂)
    R̂,ĝ
end

##
σ = 1
traces = map(1:30) do mc
    N         = 100
    Rf        = Robotlib.toOrthoNormal(randn(3,3))
    gf        = 100randn(3)
    POSES     = cat([Robotlib.toOrthoNormal(randn(3,3)) for i = 1:N]..., dims=3)
    forces    = hcat([Rf'*(POSES[1:3,1:3,i]'*gf) for i = 1:N]...)' |> copy
    forces .+= σ*randn(size(forces))
    R,g,m, Rg = calibForceIterative3(POSES,forces, trace=true)
    R,g = calibForceEigen(POSES,forces)
    Rf,gf,Rg, R,g
end
##
errors = map(traces) do (Rf, gf, Rg)
    R,g = Rg[end]
    Rerr = Rangle(R,Rf, true)
    gerr = norm(g-gf) #< 1e-10 + √3*3σ/sqrt(N)
    Rerr,gerr
end

histogram([getindex.(errors, 1) getindex.(errors, 2)], layout=2)
##
errortraces = map(traces) do (Rf, gf, Rg)
    Rerr = [Rangle(R,Rf, true) for (R,g) in Rg]
    gerr = [norm(g-gf)/norm(gf) for (R,g) in Rg] #< 1e-10 + √3*3σ/sqrt(N)
    gaerr = [vecangle(g,gf) for (R,g) in Rg]
    Rerr,gerr,gaerr
end

default(size=(800,600), grid=true, linealpha=0.3, linecolor=:black)
scales = (yscale=:identity, xscale=:identity, legend=false)
ylims = (ylims=(0,30),)
plot(getindex.(errortraces, 1); layout=(1,3), subplot=1, scales..., title="\$R\$ angle [deg]")
hline!([180], subplot=1, l=(:black,:dash))
plot!(getindex.(errortraces, 2); subplot=2, scales..., title="\$g\$ relative error")
plot!(getindex.(errortraces, 3); subplot=3, scales..., title="\$g\$ angle [deg]")

##
opts = (fillalpha=0.3,)
histogram([e[1][end] for e in errortraces]; layout=(3,1), subplot=1, title="Rangle", opts...)
histogram!([Rangle(t[1],t[4])*180/pi for t in traces]; subplot=1, opts...)

histogram!([e[2][end] for e in errortraces]; subplot=2, title="g relative", opts...)
histogram!([norm(t[2]-t[5])/norm(t[2]) for t in traces]; subplot=2, opts...)

histogram!([e[3][end] for e in errortraces]; subplot=3, title="g angle", opts...)
histogram!([vecangle(t[2],t[5]) for t in traces]; subplot=3, opts...)

# hline!([√3*3σ/sqrt(N)], l=(:black,:dash, 3), sp=2)

##
using Test
A = 10randn(100,3)
B = 10randn(100,3)

@test 1/opnorm(A) < opnorm(pinv(A))
# @show sort(abs.(eigvals(A*pinv(A)*B*pinv(B))), rev=true)

##
D = [toOrthoNormal(randn(3,3)) for _ in 1:1000]
B = [D' for D in D]
D = reduce(vcat,D)
B = reduce(vcat,B)
opnorm(D)
opnorm(B*pinv(D))
sort(eigvals(B*pinv(D)), by=abs, rev=true)
svdvals(B*pinv(D))

##
σ = 1
N         = 1000
Rf        = Robotlib.toOrthoNormal(randn(3,3))
gf        = randn(3)
gf        = gf/norm(gf)
A         = [Robotlib.toOrthoNormal(randn(3,3)) for i = 1:N]
POSES     = reduce((x,y)->cat(x,y, dims=3), A)
D         = reshape(POSES, 3,3N)'
forces    = hcat([Rf'*(POSES[1:3,1:3,i]'*gf) for i = 1:N]...)' |> copy
forces  .+= σ*randn()
F = kron(forces, Eye(3))

# FF = F*pinv(F)
# DD = D*pinv(D)

K = F\D*(D\F)
# K = pinv(F)*D*(pinv(D)*F)
# eigvals(K)
# svdvals(K)
# opnorm(K)
# rcond(K)
Rh = eigenR(K)
norm((D\F*vec(Rh)) - gf)
# Rangle(Rh,Rf)*180/pi

##
mapslices(eigen(K).vectors, dims=1) do v
    vecangle(real.(v),vec(Rf))
end

r = vec(Eye(3)) |> Vector
for i = 1:5
    global r
    # r = tor(K*r)
    r .= K*r
    r ./= norm(r)
    @show norm(toR(r) - Rf)
end
println()

R = reshape(real(eigen(K).vectors[:,1]),3,3)
R'Rf


text:
#=
    Fuel cell modelling
 =#

mutable struct FuelCell <: AbstractConverter
     # Paramètres
     α_p::Float64
     η_H2_E::Float64
     η_H2_H::Float64
     lifetime::Int64
     nHoursMax::Float64
     bounds::NamedTuple{(:lb, :ub), Tuple{Float64, Float64}}
     # Initial conditions
     powerMax_ini::Float64
     soh_ini::Float64
     # Variables
     powerMax::AbstractArray{Float64,2}
     carrier::Vector{EnergyCarrier}
     soh::AbstractArray{Float64,3}
     # Eco
     cost::AbstractArray{Float64,2}
     # Inner constructor
     FuelCell(; α_p = 8/100,
             η_H2_E = 0.4,
             η_H2_H = 0.4,
             lifetime = 14,
             nHoursMax = 10000.,
             bounds = (lb = 0., ub = 50.),
             powerMax_ini = 0.,
             soh_ini = 1.) =
             new(α_p, η_H2_E, η_H2_H, lifetime, nHoursMax, bounds, powerMax_ini, soh_ini)
end

### Preallocation
function preallocate!(fc::FuelCell, nh::Int64, ny::Int64, ns::Int64)
     fc.powerMax = convert(SharedArray,zeros(ny+1, ns)) ; fc.powerMax[1,:] .= fc.powerMax_ini
     fc.carrier = [Electricity(), Heat(), Hydrogen()]
     fc.carrier[1].power = convert(SharedArray,zeros(nh, ny, ns))
     fc.carrier[2].power = convert(SharedArray,zeros(nh, ny, ns))
     fc.carrier[3].power = convert(SharedArray,zeros(nh, ny, ns))
     fc.soh = convert(SharedArray,zeros(nh+1, ny+1, ns)) ; fc.soh[1,1,:] .= fc.soh_ini
     fc.cost = convert(SharedArray,zeros(ny, ns))
     return fc
end

### Operation dynamic
function compute_operation_dynamics!(h::Int64, y::Int64, s::Int64, fc::FuelCell, decision::Float64, Δh::Int64)
    fc.soh[h+1,y,s], fc.carrier[1].power[h,y,s], fc.carrier[2].power[h,y,s], fc.carrier[3].power[h,y,s] =
    compute_operation_dynamics(fc, (powerMax = fc.powerMax[y,s], soh = fc.soh[h,y,s]), decision, Δh)
end

function compute_operation_dynamics(fc::FuelCell, state::NamedTuple{(:powerMax, :soh), Tuple{Float64, Float64}}, decision::Float64, Δh::Int64)
    # Power constraint and correction
    fc.α_p * state.powerMax <= decision && state.soh * fc.nHoursMax / Δh > 1. ? power_E = min(decision, state.powerMax) : power_E = 0.
    # Power conversion
    power_H2 = - power_E / fc.η_H2_E
    power_H = - power_H2 * fc.η_H2_H
    # SoH computation
    soh_next = state.soh - (power_E > 0.) * Δh / fc.nHoursMax
    return soh_next, power_E, power_H, power_H2
end

### Investment dynamic
function compute_investment_dynamics!(y::Int64, s::Int64, fc::FuelCell, decision::Union{Float64, Int64})
    fc.powerMax[y+1,s], fc.soh[1,y+1,s] = compute_investment_dynamics(fc, (powerMax = fc.powerMax[y,s], soh = fc.soh[end,y,s]), decision)
end

function compute_investment_dynamics(fc::FuelCell, state::NamedTuple{(:powerMax, :soh), Tuple{Float64, Float64}}, decision::Union{Float64, Int64})
    if decision > 1e-2
        powerMax_next = decision
        soh_next = 1.
    else
        powerMax_next = state.powerMax
        soh_next = state.soh
    end
    return powerMax_next, soh_next
end


text:
struct SAGA_prox_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg}
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                  # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{R}             # stepsize
    SAG::Bool               # to activate SAG version
end

mutable struct SAGA_prox_state{R<:Real,Tx}
    a::Array{Tx}            # table of alphas
    γ::R                    # stepsize
    av::Tx                  # the running average
    z::Tx                   # the latest iteration
    # some extra placeholders
    ind::Int                # running idx set
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx                # placeholder for gradients
    w::Tx                   # input of prox
    a_old::Tx        # placeholder for previous alpha
end

function SAGA_prox_state(a, γ::R, av::Tx, z::Tx) where {R,Tx}
    return SAGA_prox_state{R,Tx}(a, γ, av, z, 1, copy(av), copy(av), copy(av), a[1])
end

function Base.iterate(iter::SAGA_prox_iterable{R,C,Tx}) where {R,C,Tx}
    N = iter.N
    # updating the stepsize
    if iter.γ === nothing
        if iter.L === nothing
            @warn "smoothness parameter absent"
            return nothing
        else
            L_M = maximum(iter.L)
            γ = 1/(5*N*L_M)
        end
    else
        γ = iter.γ # provided γ
    end
    # computing the gradients and updating the table
    a = Vector{Tx}(undef, 0)
    av = zero(iter.x0)
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], iter.x0)
        av += ∇f/N
        push!(a, iter.x0) # table of x_i
    end
    #initializing the vectors
    z, ~ = prox(iter.g, (1 - γ) .* iter.x0, γ)
    state = SAGA_prox_state(a, γ, av, z)
    return state, state
end

function Base.iterate(iter::SAGA_prox_iterable{R}, state::SAGA_prox_state{R}) where {R}

    for i = 1:iter.N # for speed in implementation
        state.ind = rand(1:iter.N) # one random number (b=1)
        gradient!(state.∇f_temp, iter.F[state.ind], state.z)
        gradient!(state.temp, iter.F[state.ind], state.a[state.ind])

        @. state.w = state.z - state.γ * (state.∇f_temp - state.temp + state.av)

        state.ind = rand(1:iter.N) # one random number (b=1)
        state.a_old .= state.a[state.ind]
        state.a[state.ind] .= state.z

        gradient!(state.∇f_temp, iter.F[state.ind], state.a_old)
        gradient!(state.temp, iter.F[state.ind], state.a[state.ind])
        @. state.av -= (state.∇f_temp - state.temp) / iter.N

        prox!(state.z, iter.g, state.w, state.γ)
    end

    return state, state
end

solution(state::SAGA_prox_state) = state.z
#TODO: minibatch


text:
if false
using Test
using GraphPlot
using Neo4j
using DistributedFactorGraphs
using Pkg
using Dates
using UUIDs
using TimeZones
using SHA

include("testBlocks.jl")

testDFGAPI = CloudGraphsDFG
testDFGAPI = LightDFG
end
# Build a basic graph.

testDFGAPI = LightDFG

##==============================================================================
## DataEntry Blobs
##==============================================================================

dfg, verts, facs = connectivityTestGraph(testDFGAPI)

dataset1 = rand(UInt8, 1000)
dataset2 = rand(UInt8, 1000)

##==============================================================================
## InMemoryDataEntry
##==============================================================================
ade,adb = addData!(InMemoryDataEntry, dfg, :x1, :random, dataset1)
gde,gdb = getData(dfg, :x1, :random)
dde,ddb = deleteData!(dfg, :x1, :random)

@test ade == gde == dde
@test adb == gdb == ddb

# @test_throws ErrorException addData!(dfg, :x2, deepcopy(ade), dataset2)
ade2,adb2 = addData!(dfg, :x2, deepcopy(ade))

ade3,adb3 = updateData!(dfg, :x2, deepcopy(ade))

@test ade == ade2 == ade3
@test adb == adb2 == adb3

deleteData!(dfg, :x2, :random)

##==============================================================================
## FileDataEntry
##==============================================================================
ade,adb = addData!(FileDataEntry, dfg, :x1, :random, "/tmp/dfgFileEntryBlob", dataset1)
gde,gdb = getData(dfg, :x1, :random)
dde,ddb = deleteData!(dfg, :x1, :random)

@test ade == gde == dde
@test adb == gdb == ddb

@test_throws ErrorException addData!(dfg, :x2, deepcopy(ade), dataset2)

ade2,adb2 = addData!(dfg, :x2, deepcopy(ade), dataset1)
ade3,adb3  = updateData!(dfg, :x2, deepcopy(ade), dataset1)

@test ade == ade2 == ade3
@test adb == adb2 == adb3

deleteData!(dfg, :x2, :random)

##==============================================================================
## FolderStore
##==============================================================================

# Create a data store and add it to DFG
ds = FolderStore{Vector{UInt8}}(:filestore, "/tmp/dfgFolderStore")
addBlobStore!(dfg, ds)

ade,adb = addData!(dfg, :filestore, :x1, :random, dataset1)
gde,gdb = getData(dfg, :x1, :random)
dde,ddb = deleteData!(dfg, :x1, :random)

@test ade == gde == dde
@test adb == gdb == ddb

ade2,adb2 = addData!(dfg, :x2, deepcopy(ade), dataset1)
# ade3,adb3 = updateData!(dfg, :x2, deepcopy(ade), dataset1)

@test ade == ade2# == ade3
@test adb == adb2# == adb3

deleteData!(dfg, :x2, :random)

#test default folder store
dfs = FolderStore("/tmp/defaultfolderstore")
@test dfs.folder == "/tmp/defaultfolderstore"
@test dfs.key == :default_folder_store
@test dfs isa FolderStore{Vector{UInt8}}

##==============================================================================
## InMemoryBlobStore
##==============================================================================

# Create a data store and add it to DFG
ds = InMemoryBlobStore()
addBlobStore!(dfg, ds)

ade,adb = addData!(dfg, :default_inmemory_store, :x1, :random, dataset1)
gde,gdb = getData(dfg, :x1, :random)
dde,ddb = deleteData!(dfg, :x1, :random)

@test ade == gde == dde
@test adb == gdb == ddb

ade2,adb2 = addData!(dfg, :x2, deepcopy(ade), dataset1)
# ade3,adb3 = updateData!(dfg, :x2, deepcopy(ade), dataset1)

@test ade == ade2# == ade3
@test adb == adb2# == adb3

deleteData!(dfg, :x2, :random)

##==============================================================================
## Unimplemented store
##==============================================================================
struct TestStore{T} <: DFG.AbstractBlobStore{T} end

store = TestStore{Int}()

@test_throws ErrorException getDataBlob(store, ade)
@test_throws ErrorException addDataBlob!(store, ade, 1)
@test_throws ErrorException updateDataBlob!(store,  ade, 1)
@test_throws ErrorException deleteDataBlob!(store, ade)
@test_throws ErrorException listDataBlobs(store)

##==============================================================================
## Unimplemented Entry Blob Crud
##==============================================================================
struct NotImplementedDE <: AbstractDataEntry end

nde = NotImplementedDE()

@test_throws ErrorException getDataBlob(dfg, nde)
@test_throws ErrorException addDataBlob!(dfg, nde, 1)
@test_throws ErrorException updateDataBlob!(dfg,  nde, 1)
@test_throws ErrorException deleteDataBlob!(dfg, nde)
@test_throws ErrorException listDataBlobs(dfg)


text:
export R_coeff, Q_coeff, impedanceM1
include("constant.jl")
#using .constant

# plan-wave reflection coefficient
"""

    R_coeff(Zc)

Calculate plan-wave reflection coefficient

# Arguments
- `Zc`: complex number, ground characteristic impedance

# Output
- `𝛹`: deg, grazing angle
- `R_mag`: coefficient magnitude
- `R_phase`: deg, phase

# Example
```julia-repl
julia> Zc = 10 + 10im
julia> 𝛹, R_mag, R_phase = R_coeff(Zc)
```

"""
function R_coeff(Zc;len::Int64=1000)
    𝛹 = range(0.0,π/2, length = len)
    R_coeff = ( sin.(𝛹) .- 1/Zc ) ./ ( sin.(𝛹) .+ 1/Zc )
    THETA = Float32.(𝛹*180/π)
    RMAG = Float32.(abs.(R_coeff))
    RPHASE = Float32.(angle.(R_coeff)*180/π)
    return THETA, RMAG, RPHASE
end

#@btime R_coeff(Zc)
#𝛹, R_mag, R_phase =   R_coeff(Zc)
#plot(𝛹,R_phase)

# spherical-wave reflection coefficient
"""

    Q_coeff(Zc,R2,f)

Calculate spherical-wave reflection coefficient

# Arguments:

- `Zc`: complex number, characteristic impedance
- `R2`: m, length of reflection ray path
- `f`: Hz, frequency

# Outputs:
- `𝛹`: deg, grazing angle
- `Q_mag`: coefficient magnitude
- `Q_phase`: deg, phase

# Example:
```julia-repl
julia> Zc = 10 + 10im
julia> 𝛹, Q_mag, Q_phase = R_coeff(Zc)
```

"""
function Q_coeff(Zc,R2,f,len=100)
    λ = c0/f
    k = 2*π/λ
    𝛹 = range(0,π/2, length = len)

    Q_coeff = randn(ComplexF64, (length(𝛹), 1))
    for i=1:length(𝛹)
        d = sqrt(1im*k*R2/2) * (1/Zc + sin(𝛹[i]))
        Fd = 1 + 1im*sqrt(π)*d*erfcx(-1im*(d))
        Q_coeff[i] = (Zc*sin(𝛹[i]) - 1 + 2*Fd) / (Zc*sin(𝛹[i]) + 1)
    end

    return Float32.(𝛹*180/π), Float32.(abs.(Q_coeff[:,1])), Float32.(angle.(Q_coeff[:,1])*180/π)

end


text:
using ModelingToolkit, Test, SparseArrays
@variables a b c

# Auxiliary Functions and Constants
get_sparsity_pattern(h::Array{Expression}) = sparse(Int64.(map(~, h .=== ModelingToolkit.Constant(0))))
get_sparsity_pattern(h::SparseMatrixCSC{Expression,Int64}) = sparse(Int64.(map(~, h .=== ModelingToolkit.Constant(0))))
get_sparsity_pattern(h::SparseVector{Expression,Int64}) = sparse(Int64.(map(~, h .=== ModelingToolkit.Constant(0))))

input = [1, 2, 3]

# ===== Dense tests =====
# Arrays of Matrices
h_dense_arraymat = [[a 1; b 0], [0 0; 0 0], [a c; 1 0]] # empty array support required
function h_dense_arraymat_julia(x)
    a, b, c = x
    return [[a[1] 1; b[1] 0], [0 0; 0 0], [a[1] c[1]; 1 0]]
end
function h_dense_arraymat_julia!(out, x)
    a, b, c = x
    out[1] .= [a[1] 1; b[1] 0]
    out[2] .= [0 0; 0 0]
    out[3] .= [a[1] c[1]; 1 0]
end

h_dense_arraymat_str = ModelingToolkit.build_function(h_dense_arraymat, [a, b, c])
h_dense_arraymat_oop = eval(h_dense_arraymat_str[1])
h_dense_arraymat_ip! = eval(h_dense_arraymat_str[2])
out_1_arraymat = [Array{Int64}(undef, 2, 2) for i in 1:3]
out_2_arraymat = [similar(x) for x in out_1_arraymat]
julia_dense_arraymat = h_dense_arraymat_julia(input)
mtk_dense_arraymat = h_dense_arraymat_oop(input)
@test_broken julia_dense_arraymat == mtk_dense_arraymat
h_dense_arraymat_julia!(out_1_arraymat, input)
h_dense_arraymat_ip!(out_2_arraymat, input)
@test out_1_arraymat == out_2_arraymat

# Arrays of 1D Vectors
h_dense_arrayvec = [[a, 0, c], [0, 0, 0], [1, a, b]] # same for empty vectors, etc.
function h_dense_arrayvec_julia(x)
    a, b, c = x
    return [[a[1], 0, c[1]], [0, 0, 0], [1, a[1], b[1]]]
end
function h_dense_arrayvec_julia!(out, x)
    a, b, c = x
    out[1] .= [a[1], 0, c[1]]
    out[2] .= [0, 0, 0]
    out[3] .= [1, a[1], b[1]]
end

h_dense_arrayvec_str = ModelingToolkit.build_function(h_dense_arrayvec, [a, b, c])
h_dense_arrayvec_oop = eval(h_dense_arrayvec_str[1])
h_dense_arrayvec_ip! = eval(h_dense_arrayvec_str[2])
out_1_arrayvec = [Vector{Int64}(undef, 3) for i in 1:3]
out_2_arrayvec = [Vector{Int64}(undef, 3) for i in 1:3]
julia_dense_arrayvec = h_dense_arrayvec_julia(input)
mtk_dense_arrayvec = h_dense_arrayvec_oop(input)
@test_broken julia_dense_arrayvec == mtk_dense_arrayvec
mtk_dense_arrayvec = @test_broken h_dense_arrayvec_oop(input)
h_dense_arrayvec_julia!(out_1_arrayvec, input)
h_dense_arrayvec_ip!(out_2_arrayvec, input)
@test out_1_arrayvec == out_2_arrayvec

# Arrays of Arrays of Matrices
h_dense_arrayNestedMat = [[[a 1; b 0], [0 0; 0 0]], [[b 1; a 0], [b c; 0 1]]]
function h_dense_arrayNestedMat_julia(x)
    a, b, c = x
    return [[[a[1] 1; b[1] 0], [0 0; 0 0]], [[b[1] 1; a[1] 0], [b[1] c[1]; 0 1]]]
end
function h_dense_arrayNestedMat_julia!(out, x)
    a, b, c = x
    out[1][1] .= [a[1] 1; b[1] 0]
    out[1][2] .= [0 0; 0 0]
    out[2][1] .= [b[1] 1; a[1] 0]
    out[2][2] .= [b[1] c[1]; 0 1]
end

h_dense_arrayNestedMat_str = ModelingToolkit.build_function(h_dense_arrayNestedMat, [a, b, c])
h_dense_arrayNestedMat_oop = eval(h_dense_arrayNestedMat_str[1])
h_dense_arrayNestedMat_ip! = eval(h_dense_arrayNestedMat_str[2])
out_1_arrayNestedMat = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]] # avoid undef broadcasting issue
out_2_arrayNestedMat = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]]
julia_dense_arrayNestedMat = h_dense_arrayNestedMat_julia(input)
mtk_dense_arrayNestedMat = h_dense_arrayNestedMat_oop(input)
@test_broken julia_dense_arrayNestedMat == mtk_dense_arrayNestedMat
h_dense_arrayNestedMat_julia!(out_1_arrayNestedMat, input)
h_dense_arrayNestedMat_ip!(out_2_arrayNestedMat, input)
@test out_1_arrayNestedMat == out_2_arrayNestedMat

# ===== Sparse tests =====
# Array of Matrices
h_sparse_arraymat = sparse.([[a 1; b 0], [0 0; 0 0], [a c; 1 0]])
function h_sparse_arraymat_julia(x)
    a, b, c = x
    return [sparse([a[1] 1; b[1] 0]), sparse([0 0; 0 0]), sparse([a[1] c[1]; 1 0])] # necessary because sparse([]) is a SparseVector, not SparseMatrix
end
function h_sparse_arraymat_julia!(out, x)
    a, b, c = x
    out[1][1, 1] = a[1]
    out[1][1, 2] = 1
    out[1][2, 1] = b[1]
    out[2] = sparse([0 0; 0 0]) # no undef constructor for SparseMatrixCSC
    out[3][1, 1] = a[1]
    out[3][1, 2] = c[1]
    out[3][2, 1] = 1
end

h_sparse_arraymat_str = ModelingToolkit.build_function(h_sparse_arraymat, [a, b, c])
h_sparse_arraymat_oop = eval(h_sparse_arraymat_str[1])
h_sparse_arraymat_ip! = eval(h_sparse_arraymat_str[2])
h_sparse_arraymat_sparsity_patterns = map(get_sparsity_pattern, h_sparse_arraymat)
out_1_arraymat = [similar(h) for h in h_sparse_arraymat_sparsity_patterns]
out_2_arraymat = [similar(h) for h in h_sparse_arraymat_sparsity_patterns] # can't do similar() because it will just be #undef, with the wrong sparsity pattern
julia_sparse_arraymat = h_sparse_arraymat_julia(input)
mtk_sparse_arraymat = h_sparse_arraymat_oop(input)
@test_broken julia_sparse_arraymat == mtk_sparse_arraymat
h_sparse_arraymat_julia!(out_1_arraymat, input)
h_sparse_arraymat_ip!(out_2_arraymat, input)
@test out_1_arraymat == out_2_arraymat

# Array of 1D Vectors
h_sparse_arrayvec = sparse.([[a, 0, c], [0, 0, 0], [1, a, b]])
function h_sparse_arrayvec_julia(x)
    a, b, c = x
    return sparse.([[a[1], 0, c[1]], [0, 0, 0], [1, a[1], b[1]]])
end
function h_sparse_arrayvec_julia!(out, x)
    a, b, c = x
    out[1][1] = a[1]
    out[1][3] = c[1]
    out[2] = sparse([0, 0, 0]) # necessary because sparsity pattern is 3 elements with 0 stored, not 0 elements
    out[3][1] = 1
    out[3][2] = a[1]
    out[3][3] = b[1]
end

h_sparse_arrayvec_str = ModelingToolkit.build_function(h_sparse_arrayvec, [a, b, c])
h_sparse_arrayvec_oop = eval(h_sparse_arrayvec_str[1])
h_sparse_arrayvec_ip! = eval(h_sparse_arrayvec_str[2])
h_sparse_arrayvec_sparsity_patterns = map(get_sparsity_pattern, h_sparse_arrayvec)
out_1_arrayvec = [similar(h) for h in h_sparse_arrayvec_sparsity_patterns]
out_2_arrayvec = [similar(h) for h in h_sparse_arrayvec_sparsity_patterns]
julia_sparse_arrayvec = h_sparse_arrayvec_julia(input)
mtk_sparse_arrayvec = h_sparse_arrayvec_oop(input)
@test_broken julia_sparse_arrayvec == mtk_sparse_arrayvec
h_sparse_arrayvec_julia!(out_1_arrayvec, input)
h_sparse_arrayvec_ip!(out_2_arrayvec, input)
@test out_1_arrayvec == out_2_arrayvec

# Arrays of Arrays of Matrices
h_sparse_arrayNestedMat = [[[a 1; b 0], [0 0; 0 0]], [[b 1; a 0], [b c; 0 1]]]
function h_sparse_arrayNestedMat_julia(x)
    a, b, c = x
    return [sparse.([[a[1] 1; b[1] 0], [0 0; 0 0]]), sparse.([[b[1] 1; a[1] 0], [b[1] c[1]; 0 1]])]
end
function h_sparse_arrayNestedMat_julia!(out, x)
    a, b, c = x
    out[1][1][1, 1] = a[1]
    out[1][1][1, 2] = 1
    out[1][1][2, 1] = b[1]
    out[1][2] = sparse([0 0; 0 0])
    out[2][1][1, 1] = b[1]
    out[2][1][1, 2] = 1
    out[2][1][2, 1] = a[1]
    out[2][2][1, 1] = b[1]
    out[2][2][1, 2] = c[1]
    out[2][2][2, 2] = 1
end

h_sparse_arrayNestedMat_str = ModelingToolkit.build_function(h_sparse_arrayNestedMat, [a, b, c])
h_sparse_arrayNestedMat_oop = eval(h_sparse_arrayNestedMat_str[1])
h_sparse_arrayNestedMat_ip! = eval(h_sparse_arrayNestedMat_str[2])
h_sparse_arrayNestedMat_sparsity_patterns = [map(get_sparsity_pattern, h) for h in h_sparse_arrayNestedMat]
out_1_arrayNestedMat = [[similar(h_sub) for h_sub in h] for h in h_sparse_arrayNestedMat_sparsity_patterns]
out_2_arrayNestedMat = [[similar(h_sub) for h_sub in h] for h in h_sparse_arrayNestedMat_sparsity_patterns]
julia_sparse_arrayNestedMat = h_sparse_arrayNestedMat_julia(input)
mtk_sparse_arrayNestedMat = h_sparse_arrayNestedMat_oop(input)
@test_broken julia_sparse_arrayNestedMat == mtk_sparse_arrayNestedMat
h_sparse_arrayNestedMat_julia!(out_1_arrayNestedMat, input)
h_sparse_arrayNestedMat_ip!(out_2_arrayNestedMat, input)
@test out_1_arrayNestedMat == out_2_arrayNestedMat

# Additional Tests
# Returning 0-element structures (corresponding to empty Jacobians)
# Arrays of Matrices
h_empty = [[a b; c 0], Array{Expression,2}(undef, 0,0)]
h_empty_str = ModelingToolkit.build_function(h_empty, [a, b, c])
h_empty_oop = eval(h_empty_str[1])
h_empty_ip! = eval(h_empty_str[2])
@test_broken h_empty_oop(input) == [[1 2; 3 0], Array{Int64,2}(undef,0,0)]
out = [Matrix{Int64}(undef, 2, 2), Matrix{Int64}(undef, 0, 0)]
h_empty_ip!(out, input) # should just not fail

# Array of Vectors
h_empty_vec = [[a, b, c, 0], Vector{Expression}(undef,0)]
h_empty_vec_str = ModelingToolkit.build_function(h_empty_vec, [a, b, c])
h_empty_vec_oop = eval(h_empty_vec_str[1])
h_empty_vec_ip! = eval(h_empty_vec_str[2])
@test_broken h_empty_vec_oop(input) == [[1, 2, 3, 0], Vector{Int64}(undef,0)]
out = [Vector{Int64}(undef, 4), Vector{Int64}(undef, 0)]
h_empty_vec_ip!(out, input) # should just not fail

# Arrays of Arrays of Matrices
h_emptyNested = [[[a b; c 0]], Array{Array{Expression, 2}}(undef, 0)] # emptyNested array of arrays
h_emptyNested_str = ModelingToolkit.build_function(h_emptyNested, [a, b, c])
h_emptyNested_oop = eval(h_emptyNested_str[1])
h_emptyNested_ip! = eval(h_emptyNested_str[2])
@test_broken h_emptyNested_oop(input) == [[[1 2; 3 0]], Array{Array{Int64, 2}}(undef, 0)]
out = [[[1 2;3 4]], Array{Array{Int64,2},1}(undef, 0)]
h_emptyNested_ip!(out, input) # should just not fail


text:
# # Blended Conditional Gradients

# The FW and AFW algorithms, and their lazy variants share one feature:
# they attempt to make primal progress over a reduced set of vertices. The AFW algorithm does this through
# away steps (which do not increase the cardinality of the active set), and the lazy variants do this through the
# use of previously exploited vertices. A third strategy that one can follow is to explicitly _blend_ Frank-Wolfe
# steps with gradient descent steps over the convex hull of the active set (note that this can be done without
# requiring a projection oracle over ``C``, thus making the algorithm projection-free). This results in the _Blended Conditional Gradient_
# (BCG) algorithm, which attempts to make as much progress as
# possible through the convex hull of the current active set ``S_t`` until it automatically detects that in order to
# make further progress it requires additional calls to the LMO.

# See also Blended Conditional Gradients: the unconditioning of conditional gradients, Braun et al, 2019, https://arxiv.org/abs/1805.07311

using FrankWolfe
using LinearAlgebra
using Random
using SparseArrays

n = 1000
k = 10000

Random.seed!(41)

matrix = rand(n, n)
hessian = transpose(matrix) * matrix
linear = rand(n)
f(x) = dot(linear, x) + 0.5 * transpose(x) * hessian * x
function grad!(storage, x)
    return storage .= linear + hessian * x
end
L = eigmax(hessian)

# We run over the probability simplex and call the LMO to get an initial feasible point:

lmo = FrankWolfe.ProbabilitySimplexOracle(1.0);
x00 = FrankWolfe.compute_extreme_point(lmo, zeros(n))

target_tolerance = 1e-5

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_accel_simplex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    hessian=hessian,
    emphasis=FrankWolfe.memory,
    L=L,
    accelerated=true,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_simplex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    hessian=hessian,
    emphasis=FrankWolfe.memory,
    L=L,
    accelerated=false,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_convex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    L=L,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

data = [trajectoryBCG_accel_simplex, trajectoryBCG_simplex, trajectoryBCG_convex]
label = ["BCG (accel simplex)", "BCG (simplex)", "BCG (convex)"]
FrankWolfe.plot_trajectories(data, label, xscalelog=true)

matrix = rand(n, n)
hessian = transpose(matrix) * matrix
linear = rand(n)
f(x) = dot(linear, x) + 0.5 * transpose(x) * hessian * x + 10
function grad!(storage, x)
    return storage .= linear + hessian * x
end
L = eigmax(hessian)

lmo = FrankWolfe.KSparseLMO(100, 100.0)
x00 = FrankWolfe.compute_extreme_point(lmo, zeros(n))

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_accel_simplex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    hessian=hessian,
    emphasis=FrankWolfe.memory,
    L=L,
    accelerated=true,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_simplex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    hessian=hessian,
    emphasis=FrankWolfe.memory,
    L=L,
    accelerated=false,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_convex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    L=L,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

data = [trajectoryBCG_accel_simplex, trajectoryBCG_simplex, trajectoryBCG_convex]
label = ["BCG (accel simplex)", "BCG (simplex)", "BCG (convex)"]
FrankWolfe.plot_trajectories(data, label, xscalelog=true)


text:
#Written By Johnathan Bizzano

include("../NeuralNetTest.jl")

const sample_size = 100
const iterations = 100
const scalefactor = .8
const crossoverrate = .6
const population_size = 30

const IsVerbose = true
const InputType = Float64
const StorageType = Float64

f(x) = MathConstants.e ^ sin(x)

function collect_data()
    inputs = zeros(Double, sample_size)
    outputs = zeros(Double, sample_size)
    for i in 1:sample_size
        inputs[i] = i
        outputs[i] = f(i)
    end
    TrainAndTestSplit(Data.ReadArrayToDataSet(inputs, outputs); trainFrac=.8, sort_randomly=true)
end

function test()
    designer = NeuralNet.ModelDesigner(input_shape=(1), error=Error.Meanabserr())
    rng_init = Initializer.RNG(-1.0:.01:1.0)
    push!(designer, Layers.Dense(100, rng_init))
    push!(designer, Layers.Dense(100, rng_init))
    push!(designer, Layers.Dense(1, rng_init))

    neuralnet = designer()
    nettrainer = trainer(neuralnet, optimizer=Optimizer.de_rand_1_bin(; population_size = population_size,
                     iterations = iterations, scalefactor = scalefactor, crossoverrate = crossoverrate, IsVerbose = IsVerbose))

    data = collect_data()
    println("Pre Training Test: f(3.0) = $(f(3.0)). Neural Net:$(neuralnet([3.0])[1])")
    train!(nettrainer, data[1]; epochs=1, testset=data[2], IsVerbose = IsVerbose)
    println("After Training Test: f(3.0) = $(f(3.0)). Neural Net:$(neuralnet([3.0])[1])")
end

test()


text:
using Base.Test
using AmplNLReader

function exercise_ampl_model(nlp :: AmplModel)
  show(STDOUT, nlp)
  print(STDOUT, nlp)

  # Perform dummy scaling.
  varscale(nlp, ones(nlp.meta.nvar))
  conscale(nlp, ones(nlp.meta.ncon))

  f = obj( nlp, nlp.meta.x0)
  g = grad(nlp, nlp.meta.x0)
  c = cons(nlp, nlp.meta.x0)
  J = jac( nlp, nlp.meta.x0)
  H = hess(nlp, nlp.meta.x0, y=ones(nlp.meta.ncon,))

  @printf("f(x0) = %f\n", f)
  @printf("∇f(x0) = "); display(g'); @printf("\n")
  @printf("c(x0) = ");  display(c'); @printf("\n")
  for j = 1 : nlp.meta.ncon
    @printf("c_%d(x0) = %8.1e\n", j, jth_con(nlp, nlp.meta.x0, j));
    @printf("∇c_%d(x0) =", j)
    display(jth_congrad(nlp, nlp.meta.x0, j)); @printf("\n")
    @printf("sparse ∇c_%d(x0) =", j)
    display(jth_sparse_congrad(nlp, nlp.meta.x0, j)); @printf("\n")
  end
  @printf("J(x0) = \n");      display(J); @printf("\n")
  @printf("∇²L(x0,y0) = \n"); display(H); @printf("\n")

  e = ones(nlp.meta.nvar)
  je = jprod(nlp, nlp.meta.x0, e)
  @printf("J(x0) * e = \n"); display(je); @printf("\n")
  jte = jtprod(nlp, nlp.meta.x0, ones(nlp.meta.ncon))
  @printf("J(x0)' * e = \n"); display(jte); @printf("\n")
  je = rand(nlp.meta.ncon)
  jprod!(nlp, nlp.meta.x0, e, je)
  @printf("J(x0) * e = \n"); display(je); @printf("\n")
  jte = rand(nlp.meta.nvar)
  jtprod!(nlp, nlp.meta.x0, ones(nlp.meta.ncon), jte)
  @printf("J(x0)' * e = \n"); display(jte); @printf("\n")
  he = hprod(nlp, nlp.meta.x0, e)
  @printf("∇²L(x0,y0) * e = \n"); display(he); @printf("\n")
  hprod!(nlp, nlp.meta.x0, e, he)
  @printf("∇²L(x0,y0) * e = \n"); display(he); @printf("\n")
  for j = 1 : nlp.meta.ncon
    Hje = jth_hprod(nlp, nlp.meta.x0, e, j)
    @printf("∇²c_%d(x0) * e = ", j); display(Hje'); @printf("\n")
  end

  ghje = ghjvprod(nlp, nlp.meta.x0, g, e)
  @printf "(∇f(x0), ∇²c_j(x0) * e) = "; display(ghje'); @printf("\n")

  write_sol(nlp, "And the winner is...", rand(nlp.meta.nvar), rand(nlp.meta.ncon))
  reset!(nlp)
end

path = dirname(@__FILE__)
#Tangi:
@printf("Tangi's test")
hs10 = AmplModel(joinpath(path, "hs10.nl"))
exercise_ampl_model(hs10)
amplmodel_finalize(hs10)

bard1 = AmplModel(joinpath(path, "flp2.nl"))
@show obj(bard1,bard1.meta.x0) cons(bard1,bard1.meta.x0) bard1.meta.x0
test=ones(4)
@show obj(bard1,test) cons(bard1,test) test
@show bard1.meta.lcon bard1.meta.ucon
@show jac(bard1,bard1.meta.x0)
amplmodel_finalize(bard1)

bard1 = AmplModel(joinpath(path, "flp2r.nl"))
@show obj(bard1,bard1.meta.x0) cons(bard1,bard1.meta.x0) bard1.meta.x0
test=ones(4)
@show obj(bard1,test) cons(bard1,test) test
@show bard1.meta.lcon bard1.meta.ucon
@show jac(bard1,bard1.meta.x0)
amplmodel_finalize(bard1)

bard1 = AmplModel(joinpath(path, "bard1_fail.nl"))
@show obj(bard1,bard1.meta.x0) cons(bard1,bard1.meta.x0) bard1.meta.x0
test=ones(bard1.meta.nvar)
@show obj(bard1,test) cons(bard1,test) test
@show bard1.meta.lcon bard1.meta.ucon
@show jac(bard1,bard1.meta.x0)
amplmodel_finalize(bard1)

@test_throws AmplException obj(bard1, bard1.meta.x0)


text:
function applylungics!(system::CVSystem,old=Dict("a"=>0),restart="no")
    if restart == "no"
        system.lungs.Pp[1] = 8*mmHgToPa;
        system.lungs.Pa[1,1] = 7.5*mmHgToPa;
        system.lungs.Va[1,1] = (system.lungs.Pa[1,1]*system.lungs.Ca[1] +
            system.lungs.V0[1]);
        system.lungs.Pa[1,2] = 7*mmHgToPa;
        system.lungs.Va[1,2] = (system.lungs.Pa[1,2]*system.lungs.Ca[2] +
            system.lungs.V0[2]);
        system.lungs.Pa[1,3] = 6.5*mmHgToPa;
        system.lungs.Va[1,3] = (system.lungs.Pa[1,3]*system.lungs.Ca[3] +
            system.lungs.V0[3]);
        system.lungs.Pv[1,1] = 4.5*mmHgToPa;
        system.lungs.Vv[1,1] = (system.lungs.Pv[1,1]*system.lungs.Cv[1] +
            system.lungs.V0[4]);
        system.lungs.Pv[1,2] = 4.5*mmHgToPa;
        system.lungs.Vv[1,2] = (system.lungs.Pv[1,2]*system.lungs.Cv[2] +
            system.lungs.V0[5]);
    elseif restart == "yes"
        system.lungs.Pp[1] = old["Pp"][end];
        system.lungs.Pa[1,1] = old["Pa"][end,1];
        system.lungs.Va[1,1] = old["Va"][end,1];
        system.lungs.Qa[1,1] = old["Qa"][end,1];
        system.lungs.Pa[1,2] = old["Pa"][end,2];
        system.lungs.Va[1,2] = old["Va"][end,2];
        system.lungs.Qa[1,2] = old["Qa"][end,2];
        system.lungs.Pa[1,3] = old["Pa"][end,3];
        system.lungs.Va[1,3] = old["Va"][end,3];
        system.lungs.Qa[1,3] = old["Qa"][end,3];
        system.lungs.Pv[1,1] = old["Pv"][end,1];
        system.lungs.Vv[1,1] = old["Vv"][end,1];
        system.lungs.Qv[1,1] = old["Qv"][end,1];
        system.lungs.Pv[1,2] = old["Pv"][end,2];
        system.lungs.Vv[1,2] = old["Vv"][end,2];
        system.lungs.Qv[1,2] = old["Qv"][end,2];
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
import TSne, Plots

function tsne(X)
    X = convert(Array{Float64,2}, X);
    Y = TSne.tsne(X, 3);

    return Y
end


text:
#module Markdown2HTML
# Markdown to HTML writer, Modified from Julia Base.Markdown html writer
using Markdown: MD, Header, Code, Paragraph, BlockQuote, Footnote,
      Admonition, List, HorizontalRule, Bold, Italic, Image, Link, LineBreak,
      LaTeX, isordered

function tohtml(io::IO, m::MIME"text/html", x)
    show(io, m, x)
end

function tohtml(io::IO, m::MIME"text/plain", x)
    htmlesc(io, sprint(show, m, x))
end

function tohtml(io::IO, m::MIME"image/png", img)
    print(io, """<img src="data:image/png;base64,""")
    print(io, stringmime(m, img))
    print(io, "\" />")
end

function tohtml(m::MIME"image/svg+xml", img)
    show(io, m, img)
end

# AbstractDisplay infrastructure

function bestmime(val)
    for mime in ("text/html", "image/svg+xml", "image/png", "text/plain")
        showable(mime, val) && return MIME(Symbol(mime))
    end
    error("Cannot render $val to Markdown.")
end

tohtml(io::IO, x) = tohtml(io, bestmime(x), x)

# Utils

function withtag(f, io::IO, tag, attrs...)
    print(io, "<$tag")
    for (attr, value) in attrs
        print(io, " ")
        htmlesc(io, attr)
        print(io, "=\"")
        htmlesc(io, value)
        print(io, "\"")
    end
    f === nothing && return print(io, " />")

    print(io, ">")
    f()
    print(io, "</$tag>")
end

tag(io::IO, tag, attrs...) = withtag(nothing, io, tag, attrs...)

const _htmlescape_chars = Dict('<'=>"<",   '>'=>">",
                               '"'=>""", '&'=>"&",
                               # ' '=>" ",
                               )
for ch in "'`!\$%()=+{}[]"
    _htmlescape_chars[ch] = "&#$(Int(ch));"
end

function htmlesc(io::IO, s::AbstractString)
    # s1 = replace(s, r"&(?!(\w+|\#\d+);)" => "&")
    for ch in s
        print(io, get(_htmlescape_chars, ch, ch))
    end
end
function htmlesc(io::IO, s::Symbol)
    htmlesc(io, string(s))
end
function htmlesc(io::IO, xs::Union{AbstractString,Symbol}...)
    for s in xs
        htmlesc(io, s)
    end
end
function htmlesc(s::Union{AbstractString,Symbol})
    sprint(htmlesc, s)
end

# Block elements

function html(io::IO, content::Vector)
    for md in content
        html(io, md)
        println(io)
    end
end

html(io::IO, md::MD) = html(io, md.content)

function html(io::IO, header::Header{l}) where l
    withtag(io, "h$l") do
        htmlinline(io, header.text)
    end
end

function html(io::IO, code::Code)
    withtag(io, :pre) do
        maybe_lang = !isempty(code.language) ? Any[:class=>"language-$(code.language)"] : []
        withtag(io, :code, maybe_lang...) do
            htmlesc(io, code.code)
            # TODO should print newline if this is longer than one line ?
        end
    end
end

function html(io::IO, md::Paragraph)
    withtag(io, :p) do
        htmlinline(io, md.content)
    end
end

function html(io::IO, md::BlockQuote)
    withtag(io, :blockquote) do
        println(io)
        html(io, md.content)
    end
end

function html(io::IO, f::Footnote)
    withtag(io, :div, :class => "footnote", :id => "footnote-$(f.id)") do
        withtag(io, :p, :class => "footnote-title") do
            print(io, f.id)
        end
        html(io, f.text)
    end
end

function html(io::IO, md::Admonition)
    withtag(io, :div, :class => "admonition $(md.category)") do
        withtag(io, :p, :class => "admonition-title") do
            print(io, md.title)
        end
        html(io, md.content)
    end
end

function html(io::IO, md::List)
    maybe_attr = md.ordered > 1 ? Any[:start => string(md.ordered)] : []
    withtag(io, isordered(md) ? :ol : :ul, maybe_attr...) do
        for item in md.items
            println(io)
            withtag(io, :li) do
                html(io, item)
            end
        end
        println(io)
    end
end

function html(io::IO, md::HorizontalRule)
    tag(io, :hr)
end

function html(io::IO, tex::LaTeX)
    withtag(io, :p, :class => "math") do
        write(io, string("\\[\n", tex.formula, "\n\\]"))
    end
end

function html(io::IO, comment::Comment)
    write(io, "\n<!-- $(comment.text) -->\n")
end

html(io::IO, x) = tohtml(io, x)

# Inline elements

function htmlinline(io::IO, content::Vector)
    for x in content
        htmlinline(io, x)
    end
end

function htmlinline(io::IO, code::Code)
    withtag(io, :code) do
        htmlesc(io, code.code)
    end
end

function htmlinline(io::IO, tex::LaTeX)
    withtag(io, :span, :class => "math") do
        write(io, string("\$", tex.formula, "\$"))
    end
end

function htmlinline(io::IO, md::Union{Symbol,AbstractString})
    htmlesc(io, md)
end

function htmlinline(io::IO, md::Bold)
    withtag(io, :strong) do
        htmlinline(io, md.text)
    end
end

function htmlinline(io::IO, md::Italic)
    withtag(io, :em) do
        htmlinline(io, md.text)
    end
end

function htmlinline(io::IO, md::Image)
    tag(io, :img, :src=>md.url, :alt=>md.alt)
end

function htmlinline(io::IO, f::Footnote)
    withtag(io, :a, :href => "#footnote-$(f.id)", :class => "footnote") do
        print(io, "[", f.id, "]")
    end
end

function htmlinline(io::IO, link::Link)
    withtag(io, :a, :href=>link.url) do
        htmlinline(io, link.text)
    end
end

function htmlinline(io::IO, br::LineBreak)
    tag(io, :br)
end

function htmlinline(io::IO, comment::Comment)
    write(io, "<!-- $(comment.text) -->")
end

htmlinline(io::IO, x) = tohtml(io, x)

# API

html(md) = sprint(html, md)

#end


text:
using Test, WeakRefStrings, Parsers
import Parsers: SENTINEL, OK, EOF, OVERFLOW, QUOTED, DELIMITED, INVALID_QUOTED_FIELD, ESCAPED_STRING, NEWLINE, SUCCESS, peekbyte, incr!, checksentinel, checkdelim, checkcmtemptylines

@testset "InlineString basics" begin

y = "abcdef"
x = InlineString(y)
x, overflow = WeakRefStrings.addcodeunit(x, UInt8('g'))
@test !overflow
@test x == "abcdefg"
x, overflow = WeakRefStrings.addcodeunit(x, UInt8('g'))
@test overflow

x = InlineString("abc")
@test x == InlineString7(x) == InlineString15(x) == InlineString31(x) == InlineString63(x)
@test x == InlineString127(x) == InlineString255(x)
y = InlineString7(x)
@test InlineString3(y) == x
@test_throws ArgumentError InlineString3(InlineString("abcd"))
@test_throws ArgumentError InlineString1(InlineString("ab"))
x = InlineString("a")
y = InlineString7(x)
@test x == y
@test InlineString1(y) == x
@test InlineString1(x) == x

@test promote_type(InlineString1, InlineString3) === InlineString3
@test promote_type(InlineString1, InlineString255) === InlineString255
@test promote_type(InlineString31, InlineString127) === InlineString127
@test promote_type(InlineString255, InlineString7) === InlineString255
@test promote_type(InlineString63, InlineString15) === InlineString63

# Ensure we haven't caused ambiguity with Base.
# https://discourse.julialang.org/t/having-trouble-implementating-a-tables-jl-row-table-when-using-badukgoweiqitools-dataframe-tbl-no-longer-works/63622/1
@test promote_type(Union{}, String) == String

end # @testset

@testset "InlineString operations" begin
    for y in ("",  "🍕", "a", "a"^3, "a"^7, "a"^15, "a"^31, "a"^63, "a"^127, "a"^255)
        x = InlineString(y)
        @show typeof(x)
        @test codeunits(x) == codeunits(y)
        @test sizeof(x) == sizeof(y)
        @test ncodeunits(x) == ncodeunits(y)
        @test length(x) == length(y)
        @test codeunit(x) == UInt8
        @test lastindex(x) == lastindex(y)
        @test isempty(x) == isempty(y)
        @test String(x) === y
        @test Symbol(x) == Symbol(y)
        @test Vector{UInt8}(x) == Vector{UInt8}(y)
        @test Array{UInt8}(x) == Array{UInt8}(y)
        @test isascii(x) == isascii(y)
        @test x * x == y * y
        @test x^5 == y^5
        @test string(x) == string(y)
        @test join([x, x]) == join([y, y])
        @test reverse(x) == reverse(y)
        y != "" && @test startswith(x, "a") == startswith(y, "a")
        y != "" && @test endswith(x, "a") == endswith(y, "a")
        y != "" && @test findfirst(==(x[1]), x) === findfirst(==(x[1]), y)
        y != "" && @test findlast(==(x[1]), x) === findlast(==(x[1]), y)
        y != "" && @test findnext(==(x[1]), x, 1) === findnext(==(x[1]), y, 1)
        y != "" && @test findprev(==(x[1]), x, length(x)) === findprev(==(x[1]), y, length(x))
        @test lpad(x, 12) == lpad(y, 12)
        @test rpad(x, 12) == rpad(y, 12)
        y != "" && @test replace(x, x[1] => 'a') == replace(y, x[1] => 'a')
        r1 = match(Regex(x), x)
        r2 = match(Regex(y), y)
        @test r1 === r2 || r1.match == r2.match
        for i = 1:ncodeunits(x)
            @test codeunit(x, i) == codeunit(y, i)
            @test isvalid(x, i) == isvalid(y, i)
            @test thisind(x, i) == thisind(y, i)
            @test nextind(x, i) == nextind(y, i)
            @test prevind(x, i) == prevind(y, i)
            @test iterate(x, i) == iterate(y, i)
        end
        for i = 1:length(x)
            if isvalid(x, i)
                @test x[i] == y[i]
            end
        end
        @test x == x
        @test x == y
        @test y == x
        @test hash(x) == hash(y)
        @test cmp(x, x) == 0
        @test cmp(x, y) == 0
        @test cmp(y, x) == 0
    end
end

@testset "InlineString parsing" begin
testcases = [
    ("", InlineString7(""), NamedTuple(), OK | EOF),
    (" ", InlineString7(" "), NamedTuple(), OK | EOF),
    (" \"", InlineString7(), NamedTuple(), OK | QUOTED | EOF | INVALID_QUOTED_FIELD), # invalid quoted
    (" \"\" ", InlineString7(), NamedTuple(), OK | QUOTED | EOF), # quoted
    (" \" ", InlineString7(), NamedTuple(), OK | QUOTED | INVALID_QUOTED_FIELD | EOF), # invalid quoted
    (" \" \" ", InlineString7(" "), NamedTuple(), OK | QUOTED | EOF), # quoted
    ("NA", InlineString7(), (; sentinel=["NA"]), EOF | SENTINEL), # sentinel
    ("\"\"", InlineString7(), NamedTuple(), OK | QUOTED | EOF), # same e & cq
    ("\"\",", InlineString7(), NamedTuple(), OK | QUOTED | EOF | DELIMITED), # same e & cq
    ("\"\"\"\"", InlineString7("\""), NamedTuple(), OK | QUOTED | ESCAPED_STRING | EOF), # same e & cq
    ("\"\\", InlineString7(), (; escapechar=UInt8('\\')), OK | QUOTED | INVALID_QUOTED_FIELD | EOF), # \\ e, invalid quoted
    ("\"\\\"\"", InlineString7("\""), (; escapechar=UInt8('\\')), OK | QUOTED | ESCAPED_STRING | EOF), # \\ e, valid
    ("\"\"", InlineString7(), (; escapechar=UInt8('\\')), OK | QUOTED | EOF), # diff e & cq
    ("\"a", InlineString7(), NamedTuple(), OK | QUOTED | INVALID_QUOTED_FIELD | EOF), # invalid quoted
    ("\"a\"", InlineString7("a"), NamedTuple(), OK | QUOTED | EOF), # quoted
    ("\"a\" ", InlineString7("a"), NamedTuple(), OK | QUOTED | EOF), # quoted
    ("\"a\",", InlineString7("a"), NamedTuple(), OK | QUOTED | EOF | DELIMITED), # quoted
    ("a,", InlineString7("a"), NamedTuple(), OK | EOF | DELIMITED),
    ("a__", InlineString7("a"), (; delim="__"), OK | EOF | DELIMITED),
    ("a,", InlineString7("a"), (; ignorerepeated=true), OK | EOF | DELIMITED),
    ("a__", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | EOF | DELIMITED),
    ("a\n", InlineString7("a"), (; ignorerepeated=true), OK | NEWLINE | EOF),
    ("a\r", InlineString7("a"), (; ignorerepeated=true), OK | NEWLINE | EOF),
    ("a\r\n", InlineString7("a"), (; ignorerepeated=true), OK | NEWLINE | EOF),
    ("a", InlineString7("a"), (; ignorerepeated=true), OK | EOF),
    ("a,,\n", InlineString7("a"), (; ignorerepeated=true), OK | DELIMITED | NEWLINE | EOF),
    ("a\n", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | NEWLINE | EOF),
    ("a\r", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | NEWLINE | EOF),
    ("a\r\n", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | NEWLINE | EOF),
    ("a", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | EOF),
    ("a____\n", InlineString7("a"), (; delim="__", ignorerepeated=true), OK | DELIMITED | NEWLINE | EOF),
    ("a\n", InlineString7("a"), NamedTuple(), OK | NEWLINE | EOF),
    ("a\r", InlineString7("a"), NamedTuple(), OK | NEWLINE | EOF),
    ("a\r\n", InlineString7("a"), NamedTuple(), OK | NEWLINE | EOF),
    ("abcdefg", InlineString7("abcdefg"), (; delim=nothing), OK | EOF),
    ("", InlineString7(), (; sentinel=missing), SENTINEL | EOF),
]

for (i, case) in enumerate(testcases)
    println("testing case = $i")
    buf, check, opts, checkcode = case
    res = Parsers.xparse(InlineString7, buf; opts...)
    @test check === res.val
    @test checkcode == res.code
end

res = Parsers.xparse(InlineString1, "")
@test Parsers.overflow(res.code)
res = Parsers.xparse(InlineString1, "ab")
@test Parsers.overflow(res.code)
res = Parsers.xparse(InlineString1, "b")
@test res.val === InlineString("b")

end # @testset


text:
###################################################
# generalized Bloch Hamiltonians that can take any
# Green's function as an argument.
###################################################
"""
    apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p, t)

Apply the generalized Bloch Hamiltonian to `m` and write the resulting derivative wrt. time into `∂m∂t`.

# Arguments
- `∂m∂t::Vector{<:Number}`: Vector describing to derivative of `m` wrt. time; this vector has to be of the same size as `m`, but can contain any value, which is replaced by `H * m`
- `m::Vector{<:Number}`: Vector the spin ensemble state of the form `[xf, yf, zf, zs, 1]` if now gradient is calculated or of the form `[xf, yf, zf, zs, 1, ∂xf/∂θ1, ∂yf/∂θ1, ∂zf/∂θ1, ∂zs/∂θ1, 0, ..., ∂xf/∂θn, ∂yf/∂θn, ∂zf/∂θn, ∂zs/∂θn, 0]` if n derivatives wrt. `θn` are calculated
- `mfun`: History fuction; can be initialized with `mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5n + 5)` for n gradients, and is then updated by the delay differential equation solvers
- `p::NTuple{9,10, or 11, Any}`: `(ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g)`, with
    -`ω1::Number`: Rabi frequency in rad/s (rotation about the y-axis)
    -`B1::Number`: B1 scaling normalized so that `B1=1` corresponds to a perfectly calibrated RF field
    -`ω0::Number`: Larmor or off-resonance frequency in rad/s
    -`m0s::Number`: Fractional semi-solid spin pool size in the range of 0 to 1
    -`R1f::Number`: Longitudinal spin relaxation rate of the free pool in 1/seconds
    -`R2f::Number`: Trasversal spin relaxation rate of the free pool in 1/seconds
    -`Rx::Number`: Exchange rate between the two pools in 1/seconds
    -`R1s::Number`: Longitudinal spin relaxation rate of the semi-solid pool in 1/seconds
    -`T2s::Number`: Trasversal spin relaxation time of the semi-solid pool in seconds
    -`g::Function`: Green's function of the form `G(κ) = G((t-τ)/T2s)`
    or `(ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, zs_idx, g)` with
    - `zs_idx::Integer`: Index to be used history function to be used in the Green's function; Default is 4 (zs), and for derivatives 9, 14, ... are used
    or `(ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s, grad_list)` with
    - `dG_o_dT2s_x_T2s::Function`: Derivative of the Green's function wrt. T2s, multiplied by T2s; of the form `dG_o_dT2s_x_T2s(κ) = dG_o_dT2s_x_T2s((t-τ)/T2s)`
    - `grad_list::Vector{<:grad_param}`: List of gradients to be calucualted; any subset of `[grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]`; length of the vector must be n (cf. arguments `m` and `∂m∂t`); ; the derivative wrt. to apparent `R1a = R1f = R1s` can be calculated with `grad_R1a()`
- `t::Number`: Time in seconds

Optional:
- `pulsetype=:normal`: Use default for a regular RF-pulse; the option `pulsetype=:inversion` should be handled with care as it is only inteded to calculate the saturation of the semi-solid pool and its derivative.

# Examples
```jldoctest
julia> using DifferentialEquations

julia> α = π/2;

julia> TRF = 100e-6;

julia> ω1 = α/TRF;

julia> B1 = 1;

julia> ω0 = 0;

julia> m0s = 0.2;

julia> R1f = 1/3;

julia> R2f = 15;

julia> R1s = 2;

julia> T2s = 10e-6;

julia> Rx = 30;

julia> G = interpolate_greens_function(greens_superlorentzian, 0, TRF / T2s);

julia> m0 = [0; 0; 1-m0s; m0s; 1];

julia> mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5);

julia> sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G)))
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 9-element Vector{Float64}:
 0.0
 1.375006182301112e-7
 1.512506800531223e-6
 8.04255964308422e-6
 2.1078484518172075e-5
 3.9114139737297395e-5
 6.268790292353144e-5
 9.147705095836351e-5
 0.0001
u: 9-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.8, 0.2, 1.0]
 [0.001727880603076341, 0.0, 0.7999981340131751, 0.19999953350448, 1.0]
 [0.019004717382235078, 0.0, 0.7997742277135814, 0.19994357804868362, 1.0]
 [0.10079108788864381, 0.0, 0.7936248155456308, 0.19842287319977142, 1.0]
 [0.26002573409211854, 0.0, 0.7565529847116781, 0.18981913439205117, 1.0]
 [0.4610419429673234, 0.0, 0.6537242534410707, 0.16937688960872904, 1.0]
 [0.666173809234659, 0.0, 0.44261243675880635, 0.13589316093594178, 1.0]
 [0.7923116716967052, 0.0, 0.10713152453821989, 0.09390269485075242, 1.0]
 [0.7994211188440913, 0.0, 0.0004403374356386769, 0.0821480968384837, 1.0]

julia> using Plots

julia> plot(sol, labels=["xf" "yf" "zf" "zs" "1"], xlabel="t (s)", ylabel="m(t)");

julia> dG_o_dT2s_x_T2s = interpolate_greens_function(dG_o_dT2s_x_T2s_superlorentzian, 0, TRF / T2s);

julia> grad_list = [grad_R2f(), grad_m0s()];

julia> m0 = [0; 0; 1-m0s; m0s; 1; zeros(5*length(grad_list))];

julia> mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? 0.0 : zeros(5 + 5*length(grad_list));

julia> sol = solve(DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G, dG_o_dT2s_x_T2s, grad_list)));

julia> plot(sol);

```
"""
function apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p::NTuple{11,Any}, t)
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, zs_idx, g = p

    ∂m∂t[1] = - R2f * m[1] - ω0  * m[2] + B1 * ω1 * m[3]
    ∂m∂t[2] =   ω0  * m[1] - R2f * m[2]
    ∂m∂t[3] = - B1 * ω1  * m[1] - (R1f + Rx * m0s) * m[3] + Rx * (1 - m0s) * m[4] + (1 - m0s) * R1f * m[5]

    if ω0 == 0
        xs = 0
        ys = quadgk(x -> g((t - x) / T2s) * mfun(p, x; idxs=zs_idx), eps(), t, order=100)[1]
    else
        xs = sin(ω0 * t) * quadgk(x -> g((t - x) / T2s) * mfun(p, x; idxs=zs_idx) * sin(ω0 * x), eps(), t, order=100)[1]
        ys = cos(ω0 * t) * quadgk(x -> g((t - x) / T2s) * mfun(p, x; idxs=zs_idx) * cos(ω0 * x), eps(), t, order=100)[1]
    end

    ∂m∂t[4] = -B1^2 * ω1^2 * (xs + ys) + Rx * m0s  * m[3] - (R1s + Rx * (1 - m0s)) * m[4] + m0s * R1s * m[5]
    return ∂m∂t
end

function apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p::NTuple{10,Any}, t)
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g = p
    return apply_hamiltonian_gbloch!(∂m∂t, m, mfun, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, 4, g), t)
end

# Version for an isolated semi-solid pool
function apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p::NTuple{6,Any}, t)
    ω1, B1, ω0, R1s, T2s, g = p

    yt = cos(ω0 * t) * quadgk(x -> g((t - x) / T2s) * mfun(p, x)[1] * cos(ω0 * x), 0, t, order=100)[1]
    xt = sin(ω0 * t) * quadgk(x -> g((t - x) / T2s) * mfun(p, x)[1] * sin(ω0 * x), 0, t, order=100)[1]

    ∂m∂t[1] = -B1^2 * ω1^2 * (xt + yt) + R1s * (1 - m[1])
end

function apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p::NTuple{12,Any}, t; pulsetype=:normal)
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s, grad_list = p

    # Apply Hamiltonian to M
    u_v1 = @view m[1:5]
    du_v1 = @view ∂m∂t[1:5]
    apply_hamiltonian_gbloch!(du_v1, u_v1, mfun, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, 4, g), t)

    # Apply Hamiltonian to all derivatives and add partial derivatives
    for i = 1:length(grad_list)
        du_v = @view ∂m∂t[5 * i + 1:5 * (i + 1)]
        u_v  = @view m[5 * i + 1:5 * (i + 1)]
        apply_hamiltonian_gbloch!(du_v, u_v, mfun, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, (5i + 4), g), t)

        if pulsetype==:normal || isa(grad_list[i], grad_T2s) || isa(grad_list[i], grad_B1)
            add_partial_derivative!(du_v, u_v1, x -> mfun(p, x; idxs=4), (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s), t, grad_list[i])
        end
    end
    return ∂m∂t
end

function apply_hamiltonian_gbloch_inversion!(∂m∂t, m, mfun, p, t)
    apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p, t; pulsetype=:inversion)
end

###################################################
# Bloch-McConnel model to simulate free precession
###################################################
function apply_hamiltonian_freeprecession!(∂m∂t, m, p::NTuple{6,Any}, t)
    ω0, m0s, R1f, R2f, Rx, R1s = p

    ∂m∂t[1] = - R2f * m[1] - ω0  * m[2]
    ∂m∂t[2] =   ω0  * m[1] - R2f * m[2]
    ∂m∂t[3] = - (R1f + Rx * m0s) * m[3] + Rx * (1 - m0s)  * m[4] + (1 - m0s) * R1f * m[5]
    ∂m∂t[4] =   Rx * m0s  * m[3] - (R1s + Rx * (1 - m0s)) * m[4] + m0s  * R1s * m[5]
    return ∂m∂t
end

function apply_hamiltonian_freeprecession!(∂m∂t, m, p::NTuple{7,Any}, t)
    ω0, m0s, R1f, R2f, Rx, R1s, grad_list = p

    # Apply Hamiltonian to M
    u_v1 = @view m[1:5]
    du_v1 = @view ∂m∂t[1:5]
    apply_hamiltonian_freeprecession!(du_v1, u_v1, (ω0, m0s, R1f, R2f, Rx, R1s), t)

    # Apply Hamiltonian to M and all its derivatives
    for i = 1:length(grad_list)
        du_v = @view ∂m∂t[5 * i + 1:5 * (i + 1)]
        u_v  = @view m[5 * i + 1:5 * (i + 1)]
        apply_hamiltonian_freeprecession!(du_v, u_v, (ω0, m0s, R1f, R2f, Rx, R1s), t)

        add_partial_derivative!(du_v, u_v1, undef, (0.0, 1.0, ω0, m0s, R1f, R2f, Rx, R1s, undef, undef, undef), t, grad_list[i])
    end
    return ∂m∂t
end

###################################################
# implementatoin of the partial derivates for
# calculationg th gradient
###################################################
function add_partial_derivative!(∂m∂t, m, mfun, p, t, grad_type::grad_m0s)
    #  _,  _,  _,   _, R1f,   _, Rx, R1s,   _, _,  _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, _, dG_o_dT2s_x_T2s = p

    ∂m∂t[3] -= Rx * m[3] + Rx * m[4] + R1f
    ∂m∂t[4] += Rx * m[3] + Rx * m[4] + R1s
    return ∂m∂t
end

function add_partial_derivative!(∂m∂t, m, mfun, p, t, grad_type::grad_R1a)
    # _, _, _, m0s, _, _, _, _, _, _, _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, _, dG_o_dT2s_x_T2s = p

    ∂m∂t[3] += - m[3] + (1 - m0s)
    ∂m∂t[4] += - m[4] + m0s
    return ∂m∂t
end

function add_partial_derivative!(∂m∂t, m, mfun, p, t, grad_type::grad_R1f)
    # _, _, _, m0s, _, _, _, _, _, _, _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, _, dG_o_dT2s_x_T2s = p

    ∂m∂t[3] += - m[3] + (1 - m0s)
    return ∂m∂t
end

function add_partial_derivative!(∂m∂t, m, mfun, p, t, grad_type::grad_R1s)
    # _, _, _, m0s, _, _, _, _, _, _, _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, _, dG_o_dT2s_x_T2s = p

    ∂m∂t[4] += - m[4] + m0s
    return ∂m∂t
end

function add_partial_derivative!(∂m∂t, m, mfun, p, t, grad_type::grad_R2f)
    ∂m∂t[1] -= m[1]
    ∂m∂t[2] -= m[2]
    return ∂m∂t
end

function add_partial_derivative!(∂m∂t, m, mfun, p, t, grad_type::grad_Rx)
    # _, _, _, m0s, _, _, _, _, _, _, _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, TRF, dG_o_dT2s_x_T2s = p

    ∂m∂t[3] += - m0s * m[3] + (1 - m0s) * m[4]
    ∂m∂t[4] +=   m0s * m[3] - (1 - m0s) * m[4]
    return ∂m∂t
end

# version for gBloch
function add_partial_derivative!(∂m∂t, m, mfun, p::Tuple{Any,Any,Any,Any,Any,Any,Any,Any,Any,Fun,Fun}, t, grad_type::grad_T2s)
    # ω1, B1, ω0,   _,   _,   _,  _,   _, T2s, _, dG_o_dT2s_x_T2s = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s = p

    yt = cos(ω0 * t) * quadgk(x -> dG_o_dT2s_x_T2s((t - x) / T2s) * mfun(x) * cos(ω0 * x), 0, t, order=100)[1]
    xt = sin(ω0 * t) * quadgk(x -> dG_o_dT2s_x_T2s((t - x) / T2s) * mfun(x) * sin(ω0 * x), 0, t, order=100)[1]

    ∂m∂t[4] -= B1^2 * ω1^2 * (xt + yt)/T2s
    return ∂m∂t
end

# version for free precession (does nothing)
function add_partial_derivative!(∂m∂t, m, mfun, p::Tuple{Any,Any,Any,Any,Any,Any,Any,Any,Any,UndefInitializer,UndefInitializer}, t, grad_type::grad_T2s)
    return ∂m∂t
end

# version for Graham's model
function add_partial_derivative!(∂m∂t, m, mfun, p::Tuple{Any,Any,Any,Any,Any,Any,Any,Any,Any,Number,Any}, t, grad_type::grad_T2s)
    # ω1, B1,  _,   _,   _,   _,  _,   _, T2s, TRF, _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, TRF, dG_o_dT2s_x_T2s = p

    df_PSD = (τ) -> quadgk(ct -> 8 / τ * (exp(-τ^2 / 8 * (3 * ct^2 - 1)^2) - 1) / (3 * ct^2 - 1)^2 + sqrt(2π) * erf(τ / sqrt(8) * abs(3 * ct^2 - 1)) / abs(3 * ct^2 - 1), 0.0, 1.0, order=100)[1]

    ∂m∂t[4] -= df_PSD(TRF / T2s) * B1^2 * ω1^2 * m[4]
    return ∂m∂t
end

# version for gBloch model
function add_partial_derivative!(∂m∂t, m, mfun, p::Tuple{Any,Any,Any,Any,Any,Any,Any,Any,Any,Fun,Fun}, t, grad_type::grad_ω0)
    # ω1, B1, ω0, _, _, _, T2s, _, g, _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s = p

    ∂m∂t[1] -= m[2]
    ∂m∂t[2] += m[1]

    yt = -sin(ω0 * t) * t * quadgk(x -> g((t - x) / T2s) * mfun(x) * cos(ω0 * x), 0, t, order=100)[1]
    yt += cos(ω0 * t) * quadgk(x -> g((t - x) / T2s) * mfun(x) * (-x) * sin(ω0 * x), 0, t, order=100)[1]

    xt = cos(ω0 * t) * t * quadgk(x -> g((t - x) / T2s) * mfun(x) * sin(ω0 * x), 0, t, order=100)[1]
    xt = sin(ω0 * t) * quadgk(x -> g((t - x) / T2s) * mfun(x) * x * cos(ω0 * x), 0, t, order=100)[1]

    ∂m∂t[4] -= B1^2 * ω1^2 * (xt + yt)
    return ∂m∂t
end

# version for free precession & Graham's model
function add_partial_derivative!(∂m∂t, m, mfun, p, t, grad_type::grad_ω0)
    ∂m∂t[1] -= m[2]
    ∂m∂t[2] += m[1]
    return ∂m∂t
end

# version for gBloch (using ApproxFun)
function add_partial_derivative!(∂m∂t, m, mfun, p::Tuple{Any,Any,Any,Any,Any,Any,Any,Any,Any,Fun,Any}, t, grad_type::grad_B1)
    # ω1, B1, ω0, _, _, _, T2s, _, g, _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g, dG_o_dT2s_x_T2s = p

    ∂m∂t[1] += ω1 * m[3]
    ∂m∂t[3] -= ω1 * m[1]

    yt = cos(ω0 * t) * quadgk(x -> g((t - x) / T2s) * mfun(x) * cos(ω0 * x), 0, t, order=100)[1]
    xt = sin(ω0 * t) * quadgk(x -> g((t - x) / T2s) * mfun(x) * sin(ω0 * x), 0, t, order=100)[1]

    ∂m∂t[4] -= 2 * B1 * ω1^2 * (xt + yt)
    return ∂m∂t
end

# version for free precession (does nothing)
function add_partial_derivative!(∂m∂t, m, mfun, p::Tuple{Any,Any,Any,Any,Any,Any,Any,Any,Any,UndefInitializer,UndefInitializer}, t, grad_type::grad_B1)
    return ∂m∂t
end

# version for Graham
function add_partial_derivative!(∂m∂t, m, mfun, p::Tuple{Any,Any,Any,Any,Any,Any,Any,Any,Any,Number,Any}, t, grad_type::grad_B1)
    # ω1, B1, _, _, _, _, T2s, _, TRF, _ = p
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, TRF, dG_o_dT2s_x_T2s = p

    f_PSD = (τ) -> quadgk(ct -> 1.0 / abs(1 - 3 * ct^2) * (4 / τ / abs(1 - 3 * ct^2) * (exp(- τ^2 / 8 * (1 - 3 * ct^2)^2) - 1) + sqrt(2π) * erf(τ / 2 / sqrt(2) * abs(1 - 3 * ct^2))), 0.0, 1.0, order=100)[1]

    ∂m∂t[1] += ω1 * m[3]
    ∂m∂t[3] -= ω1 * m[1]
    ∂m∂t[4] -= f_PSD(TRF / T2s) * 2 * B1 * ω1^2 * T2s * m[4]
    return ∂m∂t
end

###################################################
# Implementation for comparison: the super-Lorentzian
# Green's function is hard coded, which allows to
# use special solvers for the double integral
###################################################
function apply_hamiltonian_gbloch_superlorentzian!(∂m∂t, m, mfun, p::NTuple{11,Any}, t)
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, zs_idx, N = p

    gt = (t, T2s, ct) -> exp(- (t / T2s)^2 * (3 * ct^2 - 1)^2 / 8)

    function fy!(x, y, gt, mfun, p, T2s, zs_idx, t)
        for i = 1:size(x, 2)
            y[i] = gt(t - x[2,i], T2s, x[1,i]) * mfun(p, x[2,i]; idxs=zs_idx)
        end
    end

    dy1 = Cubature.pcubature_v((x, y) -> fy!(x, y, gt, mfun, p, T2s, zs_idx, t), [0.0, max(0.0, t - N * T2s)], [1.0, t])[1]

    if t > (N * T2s)
        dy2 = T2s * sqrt(2π / 3) * Cubature.pcubature(x -> mfun(p, x[1]; idxs=zs_idx) / (t - x[1]), [0.0], [t - N * T2s])[1]

        ∂m∂t[4] = -B1^2 * ω1^2 * ((dy1) + (dy2))
    else
        ∂m∂t[4] = -B1^2 * ω1^2 * (dy1)
    end

    ∂m∂t[1] = - R2f * m[1] - ω0  * m[2] + B1 * ω1 * m[3]
    ∂m∂t[2] =   ω0  * m[1] - R2f * m[2]
    ∂m∂t[3] = - B1 * ω1  * m[1] - (R1f + Rx * m0s) * m[3] +       Rx * (1 - m0s)  * m[4] + (1 - m0s) * R1f * m[5]
    ∂m∂t[4] +=             +       Rx * m0s  * m[3] - (R1s + Rx * (1 - m0s)) * m[4] +      m0s  * R1s * m[5]
    ∂m∂t[5] = 0.0
    return ∂m∂t
end

function apply_hamiltonian_gbloch_superlorentzian!(∂m∂t, m, mfun, p::NTuple{10,Any}, t)
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, N = p
    return apply_hamiltonian_gbloch_superlorentzian!(∂m∂t, m, mfun, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, 4, N), t)
end

###################################################
# Graham's spectral model
###################################################
function apply_hamiltonian_graham_superlorentzian!(∂m∂t, m, p::NTuple{10,Any}, t)
    ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s = p

    f_PSD(τ) = quadgk(ct -> 1 / abs(1 - 3 * ct^2) * (4 / τ / abs(1 - 3 * ct^2) * (exp(- τ^2 / 8 * (1 - 3 * ct^2)^2) - 1) + sqrt(2π) * erf(τ / 2 / sqrt(2) * abs(1 - 3 * ct^2))), 0.0, 1.0, order=100)[1]
    Rrf = f_PSD(TRF / T2s) * B1^2 * ω1^2 * T2s

    return apply_hamiltonian_linear!(∂m∂t, m, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, Rrf), t)
end

function apply_hamiltonian_graham_superlorentzian!(∂m∂t, m, p::NTuple{11,Any}, t)
    ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s, grad_list = p

    # Apply Hamiltonian to M
    u_v1 = @view m[1:5]
    du_v1 = @view ∂m∂t[1:5]
    apply_hamiltonian_graham_superlorentzian!(du_v1, u_v1, (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s), t)

    # Apply Hamiltonian to M and all its derivatives
    for i = 1:length(grad_list)
        du_v = @view ∂m∂t[5 * i + 1:5 * (i + 1)]
        u_v  = @view m[5 * i + 1:5 * (i + 1)]
        apply_hamiltonian_graham_superlorentzian!(du_v, u_v, (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s), t)

        add_partial_derivative!(du_v, u_v1, undef, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, TRF, undef), t, grad_list[i])
    end
    return ∂m∂t
end

function apply_hamiltonian_graham_superlorentzian_inversionpulse!(∂m∂t, m, p::NTuple{11,Any}, t)
    ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s, grad_list = p

    # Apply Hamiltonian to M
    u_v1 = @view m[1:5]
    du_v1 = @view ∂m∂t[1:5]
    apply_hamiltonian_graham_superlorentzian!(du_v1, u_v1, (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s), t)

    # Apply Hamiltonian to M and all its derivatives
    for i = 1:length(grad_list)
        du_v = @view ∂m∂t[5 * i + 1:5 * (i + 1)]
        u_v  = @view m[5 * i + 1:5 * (i + 1)]
        apply_hamiltonian_graham_superlorentzian!(du_v, u_v, (ω1, B1, ω0, TRF, m0s, R1f, R2f, Rx, R1s, T2s), t)

        if isa(grad_list[i], grad_B1) || isa(grad_list[i], grad_T2s)
            add_partial_derivative!(du_v, u_v1, undef, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, TRF, undef), t, grad_list[i])
        end
    end
    return ∂m∂t
end

function apply_hamiltonian_linear!(∂m∂t, m, p::NTuple{9,Any}, t)
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, Rrf = p

    apply_hamiltonian_freeprecession!(∂m∂t, m, (ω0, m0s, R1f, R2f, Rx, R1s), t)

    ∂m∂t[1] += B1 * ω1 * m[3]
    ∂m∂t[3] -= B1 * ω1 * m[1]
    ∂m∂t[4] -= Rrf * m[4]
    return ∂m∂t
end

function apply_hamiltonian_linear!(∂m∂t, m, p::NTuple{10,Any}, t)
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, Rrf_d, grad_list = p
    Rrf = Rrf_d[1]

     # Apply Hamiltonian to M
    u_v1 = @view m[1:5]
    du_v1 = @view ∂m∂t[1:5]
    apply_hamiltonian_linear!(du_v1, u_v1, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, Rrf), t)

     # Apply Hamiltonian to M and all its derivatives
    for i = 1:length(grad_list)
        du_v = @view ∂m∂t[5 * i + 1:5 * (i + 1)]
        u_v  = @view m[5 * i + 1:5 * (i + 1)]
        apply_hamiltonian_linear!(du_v, u_v, (ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, Rrf), t)

        add_partial_derivative!(du_v, u_v1, undef, (ω1, B1, ω0, m0s, R1f, R2f, 0, Rx, R1s, Rrf_d, undef), t, grad_list[i])
    end
    return ∂m∂t
end

##################################################################
# Sled's model
##################################################################
"""
    apply_hamiltonian_sled!(∂m∂t, m, p, t)

Apply Sled's Hamiltonian to `m` and write the resulting derivative wrt. time into `∂m∂t`.

# Arguments
- `∂m∂t::Vector{<:Number}`: Vector of length 1 describing to derivative of `m` wrt. time; this vector can contain any value, which is replaced by `H * m`
- `m::Vector{<:Number}`: Vector of length 1 describing the `zs` magnetization
- `p::NTuple{6 or 10, Any}`: `(ω1, B1, ω0, R1s, T2s, g)` for a simulating an isolated semi-solid pool or `(ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g)` for simulating a coupled spin system; with
    -`ω1::Number`: Rabi frequency in rad/s (rotation about the y-axis)
    -`B1::Number`: B1 scaling normalized so that `B1=1` corresponds to a perfectly calibrated RF field
    -`ω0::Number`: Larmor or off-resonance frequency in rad/s
    -`R1f::Number`: Longitudinal spin relaxation rate of the free pool in 1/seconds
    -`R2f::Number`: Trasversal spin relaxation rate of the free pool in 1/seconds
    -`R1s::Number`: Longitudinal spin relaxation rate of the semi-solid in 1/seconds
    -`Rx::Number`: Exchange rate between the two pools in 1/seconds
    -`T2s::Number`: Trasversal spin relaxation time in seconds
    -`g::Function`: Green's function of the form `G(κ) = G((t-τ)/T2s)`
- `t::Number`: Time in seconds

# Examples
```jldoctest
julia> using DifferentialEquations

julia> α = π/2;

julia> TRF = 100e-6;

julia> ω1 = α/TRF;

julia> B1 = 1;

julia> ω0 = 0;

julia> R1s = 2;

julia> T2s = 10e-6;

julia> G = interpolate_greens_function(greens_superlorentzian, 0, TRF / T2s);

julia> m0 = [1];

julia> sol = solve(ODEProblem(apply_hamiltonian_sled!, m0, (0, TRF), (ω1, 1, ω0, R1s, T2s, G)), Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 3-element Vector{Float64}:
 0.0
 7.475414666720001e-5
 0.0001
u: 3-element Vector{Vector{Float64}}:
 [1.0]
 [0.6313928231811974]
 [0.48953654496619214]

julia> using Plots

julia> plot(sol, labels=["zs"], xlabel="t (s)", ylabel="m(t)");

```
"""
function apply_hamiltonian_sled!(d∂m∂t, m, p::NTuple{6,Any}, t)
    ω1, B1, ω0, R1s, T2s, g = p

    xs = sin(ω0 * t) * quadgk(x -> g((t - x) / T2s) * sin(ω0 * x), 0, t, order=100)[1]
    ys = cos(ω0 * t) * quadgk(x -> g((t - x) / T2s) * cos(ω0 * x), 0, t, order=100)[1]

    d∂m∂t[1] = -B1^2 * ω1^2 * (xs + ys) * m[1] + R1s * (1 - m[1])
end

function apply_hamiltonian_sled!(∂m∂t, m, p::NTuple{10,Any}, t)
    ω1, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, g = p

    ∂m∂t[1] = - R2f * m[1] - ω0  * m[2] + B1 * ω1 * m[3]
    ∂m∂t[2] =   ω0  * m[1] - R2f * m[2]
    ∂m∂t[3] = - B1 * ω1  * m[1] - (R1f + Rx * m0s) * m[3] + Rx * (1 - m0s) * m[4] + (1 - m0s) * R1f * m[5]

    if ω0 == 0
        xs = 0
        ys = quadgk(x -> g((t - x) / T2s), 0, t, order=100)[1]
    else
        xs = sin(ω0 * t) * quadgk(x -> g((t - x) / T2s) * sin(ω0 * x), 0, t, order=100)[1]
        ys = cos(ω0 * t) * quadgk(x -> g((t - x) / T2s) * cos(ω0 * x), 0, t, order=100)[1]
    end

    ∂m∂t[4] = -B1^2 * ω1^2 * (xs + ys) * m[4] + Rx * m0s  * m[3] - (R1s + Rx * (1 - m0s)) * m[4] + m0s * R1s * m[5]
    return ∂m∂t
end


text:
module Strings

export escapehtml, tocameldash, iso8859_1_to_utf8

using ..IOExtras
import ..replace

"""
escapeHTML(i::String)

Returns a string with special HTML characters escaped: &, <, >, ", '
"""
function escapehtml(i::AbstractString)
    # Refer to http://stackoverflow.com/a/7382028/3822752 for spec. links
    o = replace(i, "&" =>"&")
    o = replace(o, "\""=>""")
    o = replace(o, "'" =>"'")
    o = replace(o, "<" =>"<")
    o = replace(o, ">" =>">")
    return o
end

"""
    tocameldash(s::String)

Ensure the first character and characters that follow a '-' are uppercase.
"""
function tocameldash(s::String)
    toUpper = UInt8('A') - UInt8('a')
    v = Vector{UInt8}(bytes(s))
    upper = true
    for i = 1:length(v)
        @inbounds b = v[i]
        if upper
            islower(b) && (v[i] = b + toUpper)
        else
            isupper(b) && (v[i] = lower(b))
        end
        upper = b == UInt8('-')
    end
    return String(v)
end

tocameldash(s::AbstractString) = tocameldash(String(s))

@inline islower(b::UInt8) = UInt8('a') <= b <= UInt8('z')
@inline isupper(b::UInt8) = UInt8('A') <= b <= UInt8('Z')
@inline lower(c::UInt8) = c | 0x20

"""
    iso8859_1_to_utf8(bytes)

Convert from ISO8859_1 to UTF8.
"""
function iso8859_1_to_utf8(bytes::Vector{UInt8})
    io = IOBuffer()
    for b in bytes
        if b < 0x80
            write(io, b)
        else
            write(io, 0xc0 | b >> 6)
            write(io, 0x80 | b & 0x3f)
        end
    end
    return String(take!(io))
end

end # module Strings