text:
#=
Feature extraction methods. Each method needs
1) Type with fields necessary for algorithm
2) function "feature" to operate on sort with type field defined above
3) any other necessary functions for extraction algorithm

=#

export FeatureTime

featureprepare(f::Feature,sort::Sorting)=nothing

#=
Temporal Waveform
=#

mutable struct FeatureTime <: Feature
end

FeatureTime(M::Int64,N::Int64)=FeatureTime()

function feature(f::FeatureTime,sort::Sorting)
    #sort.features=sub(sort.waveform,sort.dims)
    count=1
    #@inbounds for i in sort.dims
        #sort.features[count]=sort.waveforms[i]
        #count+=1
    #end
    @inbounds for i in sort.dims
        sort.features[count]=sort.p_temp[sort.cent-sort.s.win2+i-1]
        count+=1
    end
    nothing
end

mysize(feature::FeatureTime,wavelength::Int64)=wavelength

feature_type(feature::FeatureTime,in_type::DataType)=in_type

#=
online PCA

type FeaturePCA <: Feature
    oPCA::OnlineStats.OnlinePCA
end

function FeaturePCA()
    FeaturePCA(OnlineStats.OnlinePCA(window,4))
end

function FeaturePCA(win::Int64,dims::Int64)
    FeaturePCA(OnlineStats.OnlinePCA(win,dims))
end

function feature{D<:Detect,C<:Cluster,A<:Align,F<:FeaturePCA,R<:Reduction}(sort::Sorting{D,C,A,F,R})
    OnlineStats.update!(sort.f.oPCA,sort.waveforms[:,sort.numSpikes])
    sort.features[:]=sort.f.oPCA.V*sort.waveforms[:,sort.numSpikes]
    nothing
end

function mysize(feature::FeaturePCA,wavelength::Int64)
    feature.oPCA.k
end

=#

#=
Wavelet
=#
#=
type FeatureWPD <: Feature
end

function FeatureWPD(M::Int64,N::Int64)
    FeatureWPD()
end

function feature(f::FeatureWPD,sort::Sorting)
    #a=2^i where i = 1:L and L=log2(N) where N is signal length

end

function mysize(feature::FeatureWPD,wavelength::Int64)
end
=#

#=
Integral Transform

Zviagintsev et al 2006
=#
#=
type FeatureIT <: Feature
    N1::Int64
    N2::Int64
    wavemean::Float64
    itr::Int64
    tbeta::Int64
    talpha::Int64
end

FeatureIT()=FeatureIT(0,0,0.0,0,0,0)

FeatureIT(M::Int64,N::Int64)=FeatureIT()

function feature(f::FeatureIT,sort::Sorting)

    sort.features[:]=zeros(Float64,length(sort.features))
    for i=sort.f.talpha:(sort.f.talpha+sort.f.N1)
        sort.features[1]+=sort.waveforms[i,sort.numSpikes]
    end
    sort.features[1]=sort.features[1]/sort.f.N1

    for i=sort.f.tbeta:(sort.f.tbeta+sort.f.N2)
        sort.features[2]+=sort.waveforms[i,sort.numSpikes]
    end
    sort.features[2]=sort.features[2]/sort.f.N2

    nothing
end

mysize(feature::FeatureIT,wavelength::Int64)=2

function featureprepare(f::FeatureIT,sort::Sorting)

    sort.f.wavemean=mean(sort.waveforms[:,sort.numSpikes])
    tempN1=0
    tempN2=0
    thisval=0
    lastval=0
    N1=0
    N2=0
    talpha=0
    talphat=0
    tbeta=0
    tbetat=0

    for i=1:size(sort.waveforms,1)

        thisval=sign(sort.waveforms[i,sort.numSpikes]-sort.f.wavemean)
        if thisval==lastval
            if thisval==1
                tempN2+=1
                if tempN2>N2
                    N2=tempN2
                    tbeta=tbetat
                end
            else
                tempN1+=1
                if tempN1>N1
                    N1=tempN1
                    talpha=talphat
                end
            end
        else
            if thisval==1
                tempN2=1
                tbetat=i
            else
                tempN1=1
                talphat=i
            end
        end
        lastval=thisval
    end

    if sort.f.itr==0
        sort.f.N1=N1
        sort.f.N2=N2
        sort.f.talpha=talpha
        sort.f.tbeta=tbeta
    else
        sort.f.N1=round((sort.f.itr/(sort.f.itr+1))*sort.f.N1+(1/(sort.f.itr+1))*N1)
        sort.f.N2=round((sort.f.itr/(sort.f.itr+1))*sort.f.N2+(1/(sort.f.itr+1))*N2)
        sort.f.talpha=round((sort.f.itr/(sort.f.itr+1))*sort.f.talpha+(1/(sort.f.itr+1))*talpha)
        sort.f.tbeta=round((sort.f.itr/(sort.f.itr+1))*sort.f.tbeta+(1/(sort.f.itr+1))*tbeta)
    end

    if sort.f.N1+sort.f.talpha>size(sort.waveforms,1)
        sort.f.N1=size(sort.waveforms,1)-sort.f.talpha-1
    end
    if sort.f.N2+sort.f.tbeta>size(sort.waveforms,1)
        sort.f.N2=size(sort.waveforms,1)-sort.f.tbeta-1
    end

    sort.f.itr+=1

    nothing

end
=#
#=
Discrete Derivatives
=#
#=
type FeatureDD <: Feature
    inds::Array{Int64,2}
end

FeatureDD()=FeatureDD(ones(Int64,10,2))

FeatureDD(M::Int64,N::Int64)=FeatureDD(ones(Int64,N,2))

function feature(f::FeatureDD,sort::Sorting)
    for i=1:length(sort.dims)
        sort.features[i]=sort.waveform[sort.f.inds[i,2]]-sort.waveform[sort.f.inds[i,2]-sort.f.inds[i,1]]
    end
    nothing
end

function mysize(feature::FeatureDD,wavelength::Int64)
    sizeN=0
    for i=1:length(DD_inds)
        sizeN+=wavelength-DD_inds[i]
    end
    sizeN
end

function featureprepare(f::FeatureDD,sort::Sorting)
    counter=1
    counterdim=1
    for i in DD_inds
        for j=(i+1):size(sort.waveform,1)
            sort.fullfeature[counter]=sort.waveform[j]-sort.waveform[j-i]
            if counter==sort.dims[counterdim]
                sort.f.inds[counterdim,1]=i
                sort.f.inds[counterdim,2]=j
                counterdim+=1
            end
            counter+=1
        end
    end

    sort.features[:]=sort.fullfeature[sort.dims]
    nothing
end
=#
#=
#=
Curvature
=#
#=
type FeatureCurv <: Feature
end

FeatureCurv(M::Int64,N::Int64)=FeatureCurv()

function feature(f::FeatureCurv,sort::Sorting)
    V1=0.0
    V2=0.0

    for i in sort.dims
        V1=sort.waveform[i]-sort.waveform[i-1]
        V2=sort.waveform[i+1]-2*sort.waveform[i]+sort.waveform[i-1]
        sort.feature[i-1]=V2/(1+V1^2)^1.5
    end

    nothing
end

mysize(feature::FeatureCurv,wavelength::Int64)=wavelength-2

function featureprepare(f::FeatureCurv,sort::Sorting)
    V1=0.0
    V2=0.0

    for i=2:size(sort.waveform,1)-1
        V1=sort.waveform[i]-sort.waveform[i-1]
        V2=sort.waveform[i+1]-2*sort.waveform[i]+sort.waveform[i-1]
        sort.fullfeature[i-1]=V2/(1+V1^2)^1.5
    end

    nothing
end
=#
=#


text:
function Base.:*{N,ES}(lhs::Valid{N,ES}, rhs::Valid{N,ES})
    (isempty(lhs)    || isempty(rhs))    && (return Valid{N,ES}(∅))
    (isallreals(lhs) || isallreals(rhs)) && (return Valid{N,ES}(ℝp))

    #println("multiplying $lhs with $rhs")

    if roundsinf(lhs)
        infmul(lhs, rhs)
    elseif roundsinf(rhs)
        infmul(rhs, lhs)
    elseif containszero(lhs)
        zeromul(lhs, rhs)
    elseif containszero(rhs)
        zeromul(rhs, lhs)
    else
        stdmul(rhs, lhs)
    end
end

const __LHS_POS_RHS_POS = 0
const __LHS_NEG_RHS_POS = 1
const __LHS_POS_RHS_NEG = 2
const __LHS_NEG_RHS_NEG = 3

function min_not_inf{N,ES,mode1,mode2}(x::Sigmoid{N,ES,mode1}, y::Sigmoid{N,ES,mode2})
    isfinite(x) || return (mode1 == :inward_ulp) ? x : y
    isfinite(y) || return (mode2 == :inward_ulp) ? y : x
    min(x, y)
end

function max_not_inf{N,ES,mode1,mode2}(x::Sigmoid{N,ES,mode1}, y::Sigmoid{N,ES,mode2})
    isfinite(x) || return (mode1 == :inward_ulp) ? x : y
    isfinite(y) || return (mode2 == :inward_ulp) ? y : x
    max(x, y)
end

function Base.min{N,ES}(x::Sigmoid{N,ES,:inward_exact}, y::Sigmoid{N,ES,:inward_ulp})
    rx = reinterpret(Vnum{N,ES}, x)
    ry = reinterpret(Vnum{N,ES}, y)
    (rx < zero(Vnum{N,ES})) && (rx == ry) && return x
    (rx < ry) ? x : y
end
function Base.max{N,ES}(x::Sigmoid{N,ES,:inward_exact}, y::Sigmoid{N,ES,:inward_ulp})
    rx = reinterpret(Vnum{N,ES}, x)
    ry = reinterpret(Vnum{N,ES}, y)
    (rx > zero(Vnum{N,ES})) && (rx == ry) && return x
    (rx > ry) ? x : y
end

function Base.min{N,ES}(x::Sigmoid{N,ES,:outward_exact}, y::Sigmoid{N,ES,:outward_ulp})
    rx = reinterpret(Vnum{N,ES}, x)
    ry = reinterpret(Vnum{N,ES}, y)
    (rx > zero(Vnum{N,ES})) && (rx == ry) && return x
    (rx < ry) ? x : y
end
function Base.max{N,ES}(x::Sigmoid{N,ES,:outward_exact}, y::Sigmoid{N,ES,:outward_ulp})
    rx = reinterpret(Vnum{N,ES}, x)
    ry = reinterpret(Vnum{N,ES}, y)
    (rx < zero(Vnum{N,ES})) && (rx == ry) && return x
    (rx > ry) ? x : y
end

Base.min{N,ES}(x::Sigmoid{N,ES,:inward_ulp},y::Sigmoid{N,ES,:inward_exact})   = min(y, x)
Base.max{N,ES}(x::Sigmoid{N,ES,:inward_ulp},y::Sigmoid{N,ES,:inward_exact})   = max(y, x)
Base.min{N,ES}(x::Sigmoid{N,ES,:outward_ulp},y::Sigmoid{N,ES,:outward_exact}) = min(y, x)
Base.max{N,ES}(x::Sigmoid{N,ES,:outward_ulp},y::Sigmoid{N,ES,:outward_exact}) = max(y, x)

function infmul{N,ES}(lhs::Valid{N,ES}, rhs::Valid{N,ES})
    if containszero(rhs)
        return Valid{N,ES}(ℝp)
    elseif containszero(lhs)  #lhs contains zero AND infinity
        roundsinf(rhs) && return Valid{N,ES}(ℝp)

        #at this juncture, the value lhs must round both zero and infinity, and
        #the value rhs must be a standard, nonflipped double interval that is only on
        #one side of zero.

        # (100, 1) * (3, 4)     -> (300, 4)    (l * l, u * u)
        # (-1, -100) * (3, 4)   -> (-4, -300)  (l * u, u * l)
        # (100, 1) * (-4, -3)   -> (-4, -300)  (u * l, l * u)
        # (-1, -100) * (-4, -3) -> (300, 4)    (u * u, l * l)

        _state = rounds_positive(lhs) * 1 + nonpositive(rhs) * 2

        if _state == __LHS_POS_RHS_POS
            res = ((@d_lower lhs) * (@d_lower rhs)) → ((@d_upper lhs) * (@d_upper rhs))
        elseif (_state == __LHS_NEG_RHS_POS)
            res = ((@d_lower lhs) * (@d_upper rhs)) → ((@d_upper lhs) * (@d_lower rhs))
        elseif (_state == __LHS_POS_RHS_NEG)
            res = ((@d_upper lhs) * (@d_lower rhs)) → ((@d_lower lhs) * (@d_upper rhs))
        else   #state == 3
            res = ((@d_upper lhs) * (@d_upper rhs)) → ((@d_lower lhs) * (@d_lower rhs))
        end

        (@s prev(res.lower)) <= (@s res.upper) && (return Valid{N,ES}(ℝp))
        return res

    elseif roundsinf(rhs)  #now we must check if rhs rounds infinity.

        lower1 = (@d_lower lhs) * (@d_lower rhs)
        lower2 = (@d_upper lhs) * (@d_upper rhs)

        upper1 = (@d_lower lhs) * (@d_upper rhs)
        upper2 = (@d_upper lhs) * (@d_lower rhs)

        min_not_inf(lower1, lower2) → max_not_inf(upper1, upper2)
    else
        #the last case is if lhs rounds infinity but rhs is a "well-behaved" value.
        #canonical example:
        # (2, -3) * (5, 7) -> (10, -15)
        # (2, -3) * (-7, -5) -> (15, -10)

        if (rhs.lower > zero(Vnum{N,ES}))
            ((@d_lower lhs) * (@d_lower rhs)) → ((@d_upper lhs) * (@d_lower rhs))
        else
            ((@d_upper lhs) * (@d_upper rhs)) → ((@d_lower lhs) * (@d_upper rhs))
        end
    end
end

__simple_roundszero{T <: Valid}(v::T) = ((@s v.lower) < 0) & ((@s v.upper) > 0)

function zeromul{N,ES}(lhs::Valid{N,ES}, rhs::Valid{N,ES})
    #lhs and rhs guaranteed to not cross infinity.  lhs guaranteed to contain zero.
    if __simple_roundszero(rhs)
        # when rhs spans zero, we have to check four possible endpoints.
        lower1 = (@d_lower lhs) * (@d_upper rhs)
        lower2 = (@d_upper lhs) * (@d_lower rhs)
        upper1 = (@d_lower lhs) * (@d_lower rhs)
        upper2 = (@d_upper lhs) * (@d_upper rhs)

        return min_not_inf(lower1, lower2) → max_not_inf(upper1, upper2)
        # in the case where the rhs doesn't span zero, we must only multiply by the
        # extremum.
    else
        if nonpositive(rhs)
            return ((@d_lower rhs) * (@d_upper lhs)) → ((@d_lower rhs) * (@d_lower lhs))
        else
            return ((@d_lower lhs) * (@d_upper rhs)) → ((@d_upper lhs) * (@d_upper rhs))
        end
    end
end

function stdmul{N,ES}(lhs::Valid{N,ES}, rhs::Valid{N,ES})
  #both values are "reasonable."
  _state = nonpositive(lhs) * 1 + nonpositive(rhs) * 2

  if _state == __LHS_POS_RHS_POS
    ((@d_lower lhs) * (@d_lower rhs)) → ((@d_upper lhs) * (@d_upper rhs))
  elseif _state == __LHS_NEG_RHS_POS
    ((@d_lower lhs) * (@d_upper rhs)) → ((@d_upper lhs) * (@d_lower rhs))
  elseif _state == __LHS_POS_RHS_NEG
    ((@d_upper lhs) * (@d_lower rhs)) → ((@d_lower lhs) * (@d_upper rhs))
  else #__LHS_NEG_RHS_NEG
    ((@d_upper lhs) * (@d_upper rhs)) → ((@d_lower lhs) * (@d_lower rhs))
  end
end

ext: jl

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
module EchoStateNetworks

using Random
using LinearAlgebra

abstract type ESNArch end
struct NoTeacherForcing <: ESNArch end
struct TeacherForcing <: ESNArch end
export EchoStateNetwork
export train!, predict!

mutable struct EchoStateNetwork{T<:AbstractFloat}
    Nr::Integer     # the number of reservoir neurons
    Ni::Integer     # the number of input dimension
    No::Integer     # the number of output dimension
    Nl::Integer     # the number of reservoir layers
    sparsity::T
    spectral_radius::Array{T}
    noise_level::T
    leaking_rate::Array{T}
    teacher_forcing::Bool
    input_scaling::T
    inner_scaling::T
    Wr # reservoir weights
    Wi  # input weights
    Wo::Matrix{T}   # readout weights
    Wf::Matrix{T}   # feedback weights
    state
    input
    output
    prediction
    bias_scaling::T
    rng::AbstractRNG
    function EchoStateNetwork{T}(;Ni::Integer=1,
        No::Integer=1,
        Nr::Integer=100,
        Nl::Integer=1,
        sparsity::AbstractFloat=0.95,
        spectral_radius::Array{T}=repeat(0.95,Nl),
        noise_level::AbstractFloat=0.001,
        leaking_rate::Array{T}=repeat(0.95,Nl),
        teacher_forcing::Bool=true,
        input_scaling::T = 1,
        inner_scaling::T = 1,
        bias_scaling::T =0.01,
        rng::AbstractRNG=MersenneTwister(rand(UInt32)),
       ) where T<:AbstractFloat
        esn = new()
        esn.Ni = Ni
        esn.No = No
        esn.Nr = Nr
        esn.Nl = Nl
        esn.sparsity = sparsity
        esn.spectral_radius = spectral_radius
        esn.teacher_forcing = teacher_forcing
        esn.noise_level = noise_level
        esn.leaking_rate = leaking_rate
        esn.input_scaling = input_scaling
        esn.inner_scaling = inner_scaling
        esn.rng = rng

            init_weights!(esn)

        return esn
    end
end

function init_weights!(esn::EchoStateNetwork{T}) where T<:AbstractFloat
    # init reservoir weight matrix
    esn.Wr = []

    for num_layer in 1:esn.Nl
        push!(esn.Wr, rand(esn.rng, T, esn.Nr, esn.Nr) .- T(0.5))
        esn.Wr[num_layer] = (1 - esn.leaking_rate[num_layer])*Matrix{Float64}(I, esn.Nr, esn.Nr)+ esn.leaking_rate[num_layer]*esn.Wr[num_layer]
        # reduce connections based on `sparsity`
        Ns = round(Int,length(esn.Wr[num_layer])*esn.sparsity)# 向上去数， 9.5取10的这种意思

        for i in [rand(esn.rng, 1:length(esn.Wr[num_layer])) for i in 1:Ns]# 中括号里面的意思是生成Wr个0-1的随机数，然后取前Ns个
            Temp_Wr = esn.Wr[num_layer]
            Temp_Wr[i] = zero(T)
            esn.Wr[num_layer] = Temp_Wr #这个意思是将其中sparsity里面的的数赋值为0，进行矩阵稀疏化
        end # end for i
        ## rescale the matrix to fit the `spectral radius`
        esn.Wr[num_layer] *= esn.spectral_radius[num_layer]/maximum(abs.(eigvals(esn.Wr[num_layer])))#所有元素自乘后面这一大串
        esn.Wr[num_layer] *= esn.inner_scaling
    end # end for num_layer

    # init input weight matrix
    esn.Wi =[]
    for num_layer in 1:esn.Nl
        if num_layer==1
            push!(esn.Wi,rand(esn.rng, T, esn.Nr, esn.Ni+1) .- T(0.5))
            esn.Wi[num_layer] *= esn.input_scaling
        else
            push!(esn.Wi,rand(esn.rng, T, esn.Nr, esn.Nr+1) .- T(0.5))
            esn.Wi[num_layer] *=esn.input_scaling
        end
    end
    # init feedback weight matrix 这里没有进行更改，因为用不到feedback
    esn.Wf = rand(esn.rng, T, esn.Nr, esn.No) .- T(0.5)

    return esn
end # end for function init_weights!

function activate(f::Function, x::Array)
    return map(f, x)
end

function update(esn::EchoStateNetwork{T}, state::Vector{T},
                input::Vector{T}, output::Vector{T}, arch::NoTeacherForcing, i::Integer) where T<:AbstractFloat
    return activate(tanh, esn.Wi[i]*[one(T)*esn.bias_scaling;input] + esn.Wr[i]*state)
    + esn.noise_level*(rand(esn.rng, T, esn.Nr)-T(0.5))
end

function update(esn::EchoStateNetwork{T}, state::Vector{T},
                input::Vector{T}, output::Vector{T}, arch::TeacherForcing, i::Integer) where T<:AbstractFloat
    return activate(tanh, esn.Wi[i]*[one(T)*esn.bias_scaling;input] + esn.Wr[i]*state + esn.Wf*output)
    + esn.noise_level*(rand(esn.rng, T, esn.Nr)-T(0.5))
end

function reservoir_states(esn::EchoStateNetwork{T}, inputs, outputs::Matrix{T}) where T<:AbstractFloat
    u = []
    Nd = size(inputs, 2)
    states = zeros(T, (esn.Nr, Nd+1,esn.Nl))
    for t in 2:(Nd+1)
        arch = esn.teacher_forcing ? TeacherForcing() : NoTeacherForcing()
        for i in 1:esn.Nl
            if i==1
                u = inputs[:,t-1,i]
                #从第一个输入开始向网络输入数据
            else
                u = states[:,t,i-1]
                #从第二个i-1的第t时刻数据作为>2层数的输入数据
            end
            states[:,t,i] = (one(T)-esn.leaking_rate[i])*states[:,t-1,i] + esn.leaking_rate[i]*update(esn, states[:,t-1,i], u, outputs[:,t-1], arch,i)
        end
    end
    return states
end

function train!(esn::EchoStateNetwork{T}, inputs::Matrix{T}, outputs::Matrix{T};
                                  discard::Integer=min(div(size(inputs,2),10), 100), reg::AbstractFloat=1e-8) where T<:AbstractFloat
    #计算每一层的state
    @assert(size(inputs, 1) == esn.Ni)
    @assert(size(outputs, 1) == esn.No)
    @assert(size(inputs, 2) == size(outputs, 2))
    Nd = size(inputs,2)
    inputs_for_reservior = []
    #第一维是维度数ni，第二维是时刻数t，第三维是层数nl
    push!(inputs_for_reservior,inputs)
    esn.state = reservoir_states(esn, inputs_for_reservior[1], outputs)
    for i in 1:(esn.Nl-1)
        push!(inputs_for_reservior,esn.state[:,:,i])
    end

    # extended system states
    X = zeros(esn.Nl*esn.Nr,Nd)
    for t in 1:size(inputs,2)
        for num_layer in 1:esn.Nl
            X[1+(num_layer-1)*esn.Nr:esn.Nr+(num_layer-1)*esn.Nr,t] = esn.state[:,t+1,num_layer]
        end
    end
    # discard initial transient
    Xe = vcat(X[:,discard+1:end],ones(1,Nd-discard))
    tXe = Xe'
    O = outputs[:,discard+1:end]

    # calc output weight matrix
    esn.Wo = O*tXe*pinv(Xe*tXe + reg*Matrix{T}(I, size(Xe,1), size(Xe,1)))

    # store last states
    return esn
end

function predict!(esn::EchoStateNetwork{T}, inputs::Matrix{T}, outputs::Matrix{T} ; cont::Bool=true) where T<:AbstractFloat
    Nd = size(inputs, 2)
    states = zeros(T, (esn.Nr, Nd+1,esn.Nl))
    if cont
        for i in esn.Nl
            states[:,1,i] = esn.state[:,end,i]
        end
    else
        inputs = hcat(zeros(T,esn.Ni), inputs)
        outputs = hcat(zeros(T,esn.No), outputs)
    end

    arch = esn.teacher_forcing ? TeacherForcing() : NoTeacherForcing()
    for t = 2:Nd+1
        for i in esn.Nl
            if i ==1
                u = inputs[:,t-1]
            else
                u = states[:,t,i-1]
            end
            states[:,t,i] = (one(T)-esn.leaking_rate[i])*states[:,t-1,i] + esn.leaking_rate[i]*update(esn, states[:,t-1,i], u, outputs[:,t-1], arch,i)
        end
    end

    X = zeros(esn.Nl*esn.Nr,Nd)
    for t in 1:size(inputs,2)
        for num_layer in 1:esn.Nl
            X[1+(num_layer-1)*esn.Nr:esn.Nr+(num_layer-1)*esn.Nr,t] = states[:,t+1,num_layer]
        end
    end
    Xi = vcat(X,ones(1,Nd).*esn.bias_scaling)
    esn.prediction = esn.Wo*Xi
    return esn.prediction
end

end # module

max_stars_repo_name: Ziqiang-IRCN/Echo-State-Network

max_issues_repo_name: Ziqiang-IRCN/Echo-State-Network

max_forks_repo_name: Ziqiang-IRCN/Echo-State-Network


text:
module Similitude

#   This file is part of Similitude.jl
#   It is licensed under the MIT license
#   Similitude Copyright (C) 2020 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com
#     _     _         _  _  _                   _
#    | |   (_)       (_)| |(_) _               | |
#     \ \   _  ____   _ | | _ | |_   _   _   _ | |  ____
#      \ \ | ||    \ | || || ||  _) | | | | / || | / _  )
#  _____) )| || | | || || || || |__ | |_| |( (_| |( (/ /
# (______/ |_||_|_|_||_||_||_| \___) \____| \____| \____)

import Base: @pure
import UnitSystems
import UnitSystems: UnitSystem, universe, Coupling, logdb, expdb, dB
import UnitSystems: Systems, Dimensionless, Constants, Physics, Convert, Derived
const dir = dirname(pathof(UnitSystems))

using LinearAlgebra
import AbstractTensors: TupleVector, Values, value

include("group.jl")
include("dimension.jl")
include("constant.jl")

dimtext(u) = isq
printdims(io::IO,x::Group{T,dims},u::UnitSystem) where T = printdims(io,x,dimtext(normal(u)))
printdims(io::IO,x::Group{T,vals},u::UnitSystem) where T = printdims(io,x,basis)
printdims(io::IO,x::Group{T,dims}) where T = printdims(io,x,isq)
printdims(io::IO,x::Group{T,vals}) where T = printdims(io,x,basis)

logdb(x::Dimension{D}) where D = Dimension{logdb(D)}()
logdb(x::Quantity{D,U}) where {D,U} = Quantity{logdb(D),U}(logdb(x.v))

@pure UnitSystems.unit(x::Dimension,y=1) = x
@pure UnitSystems.unit(x::AbelianGroup,y=1) = x
@pure UnitSystems.unit(x::ConvertUnit,y=1) = x

UnitSystems.unit(x::Group,y=1) = x
UnitSystems.unit(x::Quantity,y=1) = x

for unit ∈ (Systems...,Dimensionless...,Constants...,Physics...,Derived...,Convert...)
    unit ∉ (:length,:time,:angle) && @eval export $unit
end

@pure ratio_calc(D::ExpGroup{B},U,S) where B = B^ratio_calc(value(D),U,S)
@pure ratio_calc(D::LogGroup{B},U,S) where B = log(B,ratio_calc(value(D),U,S))
@pure ratio_calc(D::LogGroup{ℯ},U,S) = log(ratio_calc(value(D),U,S))
@pure ratio_calc(D::ExpGroup{ℯ},U,S) = exp(ratio_calc(value(D),U,S))
@pure ratio(::Dimension{D},U,S) where D = ratio_calc(UnitSystem(D),normal(U),normal(S))
@pure function ratio_calc(D::Group,U,S)
    Constant(
    UnitSystems.boltzmann(U,S)^D.v[1]*
    UnitSystems.planckreduced(U,S)^D.v[2]*
    UnitSystems.lightspeed(U,S)^D.v[3]*
    UnitSystems.permeability(U,S)^D.v[4]*
    UnitSystems.electronmass(U,S)^D.v[5]*
    UnitSystems.molarmass(U,S)^D.v[6]*
    UnitSystems.luminousefficacy(U,S)^D.v[7]*
    UnitSystems.angle(U,S)^D.v[8]*
    UnitSystems.rationalization(U,S)^D.v[9]*
    UnitSystems.lorentz(U,S)^D.v[10]*
    UnitSystems.gravity(U,S)^D.v[11])
end
export Ratio
const Ratio = ratio

import UnitSystems: normal, unitname
normal(x::Quantity) = quantity(x)
(s::UnitSystem)(q::Quantity{D,U}) where {D,U} = (S=normal(s); Quantity{D,S}(q.v*ratio(D,U,S)))
(q::Quantity{D,U})(s::UnitSystem) where {D,U} = (S=normal(s); Quantity{D,S}(q.v*ratio(D,U,S)))

function Quantity(u::UnitSystem)
    U = Constant(u); t = isone(normal(UnitSystems.boltzmann(U)))
    kB = Quantity{F*L/Θ,U}(UnitSystems.boltzmann(U))
    ħ = Quantity{F*L*T/A,U}(UnitSystems.planckreduced(U))
    c = Quantity{L/T,U}(UnitSystems.lightspeed(U))
    μ0 = Quantity{F*T*T*C*C/(Q*Q)/R,U}(UnitSystems.permeability(U))
    mₑ = Quantity{M,U}(UnitSystems.electronmass(U))
    Mᵤ = Quantity{M/N,U}(UnitSystems.molarmass(U))
    Kcd = Quantity{J*T/F/L,U}(UnitSystems.luminousefficacy(U))
    θ = Quantity{A,U}(UnitSystems.angle(U))
    λ = Quantity{R,U}(UnitSystems.rationalization(U))
    αL = Quantity{inv(C),U}(UnitSystems.lorentz(U))
    g₀ = Quantity{M*L/(F*T*T),U}(UnitSystems.gravity(U))
    τ = UnitSystems.tau(U)
    x = UnitSystems.two(U)
    y = UnitSystems.three(U)
    z = UnitSystems.five(U)
    u = UnitSystems.seven(U)
    v = UnitSystems.eleven(U)
    w = UnitSystems.nineteen(U)
    q = UnitSystems.fourtythree(U)
    UnitSystem(kB,ħ,c,μ0,mₑ,Mᵤ,Kcd,θ,λ,αL,g₀,Universe,τ,x,y,z,u,v,w,q)
end

const LD,JD = Constant(384399)*(𝟐*𝟓)^3,Constant(778479)*(𝟐*𝟓)^6
const μE☾ = Constant(UnitSystems.μE☾)

import UnitSystems: GaussSystem, EntropySystem, ElectricSystem, AstronomicalSystem
include("$dir/initdata.jl")

const Unified = Quantity(UnitSystem(F*L/Θ,F*L*T/A,L/T,F*T*T*C*C/(Q*Q)/R,M,M/N,J*T/F/L,A,R,inv(C),M*L/(F*T*T),Universe,τ,𝟐,𝟑,𝟓,𝟕,𝟏𝟏,𝟏𝟗,𝟒𝟑))
dimtext(::typeof(normal(Unified))) = Values("kB","ħ","𝘤","μ₀","mₑ","Mᵤ","Kcd","ϕ","λ","αL","g₀")

export Unified
unitname(::typeof(normal(Unified))) = "Unified"
(u::typeof(normal(Unified)))(d::Group) = UnitSystem(d)
(u::UnitSystem)(d::Group) = normal(Metric)(d)

for unit ∈ Convert
    if unit ∉ (:length,:time,:angle,:molarmass,:luminousefficacy)
        @eval const  $unit = dimensions(UnitSystems.$unit(UnitSystems.Natural,Natural))
    end
end
const gravityforce = acceleration/specificforce

@doc """
    (D::Dimension)(U::UnitSystem,S::UnitSystem) = ConvertUnit{D,U,S}()

Constant for unit conversion for `D::Dimension` from `U::UnitSystem` to `S::UnitSystem`.
```Julia
julia> energy(Metric,CGS)
$(energy(Metric,CGS))

julia> energy(Metric,English)
$(energy(Metric,English))
```
There still exists further opportunity to expand on the implementation of `ConvertUnit`.
""" ConvertUnit

@doc """
    (U::UnitSystem)(v::Number, D::Dimension) ↦ Quantity(D,U,v) = Quantity{D,U}(v)

Numerical `Quantity` having value `v` with `D::Dimension` specified in `U::UnitSystem`.
```Julia
julia> Metric(1,energy)
$(Metric(1,energy))

julia> English(1,energy)
$(English(1,energy))
```
An alternate syntax `Quantity(D::Dimension, U::UnitSystem, v::Number)` is also available as standard syntax.
When `using UnitSystems` instead of `using Similitude`, this same syntax can be written so that code doesn't need to be changed while the output is generated.
""" Quantity

include("derived.jl")

(u::typeof(normal(UnitSystems.QCD)))(d::Group) = normal(Planck)(d)

for dim ∈ (:angle, :length, :time)
    @eval begin
        Base.sqrt(x::typeof($dim)) = sqrt(evaldim(x))
        Base.cbrt(x::typeof($dim)) = cbrt(evaldim(x))
        Base.inv(x::typeof($dim)) = inv(evaldim(x))
        Base.log(x::typeof($dim)) = log(evaldim(x))
        Base.log2(x::typeof($dim)) = log2(evaldim(x))
        Base.log10(x::typeof($dim)) = log10(evaldim(x))
        Base.log(b::Number,x::typeof($dim)) = log(b,evaldim(x))
        Base.exp(x::typeof($dim)) = exp(evaldim(x))
        Base.exp2(x::typeof($dim)) = exp2(evaldim(x))
        Base.exp10(x::typeof($dim)) = exp10(evaldim(x))
        Base.:^(b::T,x::typeof($dim)) where T<:Number = b^evaldim(x)
        Base.:^(a::typeof($dim),b::Number) = evaldim(a)^b
        Base.:^(a::typeof($dim),b::Integer) = evaldim(a)^b
        Base.:^(a::typeof($dim),b::Rational) = evaldim(a)^b
        Base.:*(a::typeof($dim),b::Dimension) = evaldim(a)*b
        Base.:/(a::typeof($dim),b::Dimension) = evaldim(a)/b
        Base.:*(a::Dimension,b::typeof($dim)) = a*evaldim(b)
        Base.:/(a::Dimension,b::typeof($dim)) = a/evaldim(b)
        Base.:*(a::Number,b::typeof($dim)) = a*evaldim(b)
        Base.:*(a::typeof($dim),b::Number) = evaldim(a)*b
        Base.:/(a::Number,b::typeof($dim)) = a*inv(evaldim(b))
        Base.:/(a::typeof($dim),b::Number) = evaldim(a)*inv(b)
    end
    for dim2 ∈ (:angle, :length, :time)
        @eval begin
            Base.:*(a::typeof($dim),b::typeof($dim2)) = evaldim(a)*evaldim(b)
            Base.:/(a::typeof($dim),b::typeof($dim2)) = evaldim(a)*evaldim(b)
        end
    end
end

include("$dir/kinematicdocs.jl")
include("$dir/electromagneticdocs.jl")
include("$dir/thermodynamicdocs.jl")
include("$dir/physicsdocs.jl")
include("$dir/systems.jl")

@doc """
    Dimension{D} <: AbstractModule

Physical `Dimension` represented by `Group` element `D`.
```Julia
F, M, L, T, Q, Θ, N, J, A, R, C
```
Operations on `Dimension` are closed (`*`, `/`, `+`, `-`, `^`).
```Julia
julia> force(Unified)
$(force(Unified))

julia> mass(Unified)
$(mass(Unified))

julia> length(Unified)
$(length(Unified))

julia> time(Unified)
$(time(Unified))

julia> charge(Unified)
$(charge(Unified))

julia> temperature(Unified)
$(temperature(Unified))

julia> molaramount(Unified)
$(molaramount(Unified))

julia> luminousflux(Unified)
$(luminousflux(Unified))

julia> angle(Unified)
$(angle(Unified))

julia> rationalization(Unified)
$(rationalization(Unified))

julia> lorentz(Unified)
$(lorentz(Unified))
```
Derived `Dimension` can be obtained from multiplicative base of 11 fundamental `Dimension` symbols corresponding to `force`, `mass`, `length`, `time`, `charge`, `temperature`, `molaramount`, `luminousflux`, `angle`, `demagnetizingfactor`, and a `nonstandard` dimension.

Mechanics: `angle`, `$(listext(Kinematic))`, `$(listext(Mechanical))`;
Electromagnetics: `$(listext(Electromagnetic))`;
Thermodynamics: `$(listext(Thermodynamic))`,
`$(listext(Molar))`, `$(listext(Photometric))`.
""" Dimension

@doc """
    Unified = UnitSystem(...) # Unified System of Quantities (USQ)
    F, M, L, T, Q, Θ, N, J, A, R, C # fundamental base dimensions

Standard `Unified` system of `Quantities` (USQ) in terms of `UnitSystem` basis, transformed from the basis of `force`, `mass`, `length`, `time`, `charge`, `temperature`, `molaramount`, `luminousflux`, `angle`, `rationalization`, and a `nonstandard` dimension.

```Julia
julia> boltzmann(Unified) # entropy
$(boltzmann(Unified))

julia> planckreduced(Unified) # angularmomentum
$(planckreduced(Unified))

julia> lightspeed(Unified) # speed
$(lightspeed(Unified))

julia> vacuumpermeability(Unified) # permeability
$(vacuumpermeability(Unified))

julia> electronmass(Unified) # mass
$(electronmass(Unified))

julia> molarmass(Unified) # molarmass
$(molarmass(Unified))

julia> luminousefficacy(Unified) # luminousefficacy
$(luminousefficacy(Unified))

julia> radian(Unified) # angle
$(radian(Unified))

julia> rationalization(Unified) # demagnetizingfactor
$(rationalization(Unified))

julia> lorentz(Unified) # nonstandard
$(lorentz(Unified))

julia> gravity(Unified) # gravityforce
$(gravity(Unified))
```
""" Unified

end # module


text:
@testset "Test eigvals" begin

    function f(a1::T,b1,c1) where T
        a = Array{T, 2}(undef, 1, 1)
        a[1,1] = a1
        b = Array{T, 2}(undef, 1, 1)
        b[1,1] = b1
        c = Array{T, 2}(undef, 1, 1)
        c[1,1] = c1
        @test all(.≈(DeconvOptim.eigvals_symmetric_tullio(a,b,c), DeconvOptim.eigvals_symmetric(a,b,c)))
    end

    f(10.0, 20.0, -10.0)
    f(0f0, -12f0, 13f0)
end

@testset "Schatten norm consistent" begin

    x = [1 2 3; 1 1 1; 0 0 -1f0]
    @test DeconvOptim.HSp(x, p = 1) ≈ 0.9999999900000001
    @test DeconvOptim.HSp(x, p = 2) ≈ 1.732050831647567
    @test abs.(DeconvOptim.HSp(x, p = 1)) ≈ DeconvOptim.HS1(x)
end


text:
function initialize(::Type{TransmissionRates},
                    states::Vector{DiseaseState},
                    pop::Population,
                    rf::RiskFunctions{T},
                    rp::RiskParameters{T}) where T <: EpidemicModel
  n_ids = length(states)
  tr = TransmissionRates(n_ids)
  for i in findall(states .== Ref(State_S))
    # External exposure
    #tr.external[i] = rf.susceptibility(rp.susceptibility, pop, i) * rf.sparks(rp.sparks, pop, i)
    tr.external[i] = rf.sparks(rp.sparks, pop, i)
    # Internal exposure
    susceptibility = rf.susceptibility(rp.susceptibility, pop, i)
    for k in findall(states .== Ref(State_I))
      tr.internal[k, i] = susceptibility *
                          rf.infectivity(rp.infectivity, pop, i, k) *
                          rf.transmissibility(rp.transmissibility, pop, k)
    end
  end
  @debug "Initialization of $T TransmissionRates complete" ∑external = sum(tr.external) ∑internal = sum(tr.internal)
  return tr
end

function initialize(::Type{EventRates},
                    tr::TransmissionRates,
                    states::Vector{DiseaseState},
                    pop::Population,
                    rf::RiskFunctions{T},
                    rp::RiskParameters{T}) where T <: EpidemicModel
  n_ids = length(states)
  rates = EventRates{T}(n_ids)
  for i = 1:n_ids
    if states[i] == State_S
      if T in [SEIR; SEI]
        rates.exposure[i] = sum(tr, i)
      elseif T in [SIR; SI]
        rates.infection[i] = sum(tr, i)
      end
    elseif states[i] == State_E
      rates.infection[i] = rf.latency(rp.latency, pop, i)
    elseif states[i] == State_I
      if T in [SEIR; SIR]
        rates.removal[i] = rf.removal(rp.removal, pop, i)
      end
    end
  end
  @debug "Initialization of $T EventRates complete" rates = [sum(rates[k]) for k in _state_progressions[T][2:end]]
  return rates
end


text:
@testset "SRS" begin
    @testset "Default" begin
        tester = ImputorTester(SRS)
        test_hashing(tester)
        test_equality(tester)
        test_vector(tester)
        test_matrix(tester)
        test_dataframe(tester)
        # Behaviour is inconsistent for testing because of `rand` calls
        # test_groupby(tester)
        # test_axisarray(tester)
        # test_nameddimsarray(tester)
        # test_keyedarray(tester)
        test_columntable(tester)
        test_rowtable(tester)
    end

    @testset "Floats" begin
        a = allowmissing(1.0:1.0:20.0)
        a[[2, 3, 7]] .= missing
        result = impute(a, SRS(; rng=SequentialRNG()))
        expected = copy(a)
        expected[2] = 4.0
        expected[3] = 5.0
        expected[7] = 6.0

        @test result == expected

        @test result == Impute.srs(a; rng=SequentialRNG())

        a2 = copy(a)

        Impute.srs!(a2; rng=SequentialRNG())
        @test a2 == result
    end

    @testset "Ints" begin
        a = allowmissing(1:1:20)
        a[[2, 3, 7]] .= missing
        result = impute(a, SRS(; rng=SequentialRNG()))
        expected = copy(a)
        expected[2] = 4
        expected[3] = 5
        expected[7] = 6

        @test result == expected
    end

    @testset "Strings" begin
        a = allowmissing([randstring(4) for i in 1:20])
        a[[2, 3, 7]] .= missing
        result = impute(a, SRS(; rng=SequentialRNG()))
        expected = copy(a)
        expected[2] = expected[4]
        expected[3] = expected[5]
        expected[7] = expected[6]

        @test result == expected
    end
end


text:
# The InteratomicPotentials.jl API specification.

export energy_and_force, potential_energy, force, virial, virial_stress
export get_rcutoff, get_species
export get_parameters, set_parameters, serialize_parameters, deserialize_parameters
export get_hyperparameters, set_hyperparameters, serialize_hyperparameters, deserialize_hyperparameters

"""
    energy_and_force(s::AbstractSystem, p::AbstractPotential)::NamedTuple{(:e, :f), Tuple{Unitful.Energy,Vector{SVector{3, Unitful.Force}}}}

Calculate the unit-annotated potential energy of a system and the force acting on each particle in a system using the provided interatomic potential.
This combined function is offered because it is usually more efficient to calculate both properties simultaneously.
"""
function energy_and_force end
"""
    potential_energy(s::AbstractSystem, p::AbstractPotential)::Unitful.Energy

Calculate the unit-annotated potential energy of a system using the provided interatomic potential.
The default implementation uses the `:e` property of `energy_and_force(s,p)`.
"""
function potential_energy end
"""
    force(s::AbstractSystem, p::AbstractPotential)::Vector{SVector{3, Unitful.Force}}

Calculate the unit-annotated force acting on each particle in a system using the provided interatomic potential.
The default implementation uses the `:f` property of `energy_and_force(s,p)`.
"""
function force end
"""
    virial_stress(s::AbstractSystem, p::AbstractPotential)::SVector{6,Unitful.Energy}

Calculate the unit-annotated virial stress tensor of a system, officially calculated as the sum of radial-force outerproducts: `` \\sum r_{ij} \\bigotimes F_{ij}``, only returns the unique lower-diagonal components.
"""
function virial_stress end
"""
    virial(s::AbstractSystem, p::AbstractPotential)::Unitful.Energy

Calculate the unit-annotated virial of a system, officially calculated as the trace contraction of the sum of radial-force outerproducts: ``tr\\left( \\sum r_{ij} \\bigotimes F_{ij} \\right)``
"""
function virial end
"""
    get_rcutoff(p::AbstractPotential)::AbstractFloat

Retrieve the radius cutoff for the provided potential.
This is the cutoff used for neighbor list calculations.
(i.e. Any pairs beyond this cutoff will be ignored.)
Defaults to `Inf` if a potential type does not implement a custom method.
"""
function get_rcutoff end
"""
    get_species(p::AbstractPotential)::Union{Tuple,Missing}

Retrieve the species to be included in an interaction
(pairs including a species not in the list are ignored).
A value of `missing` indicates that all species should be included,
which is the default behavior if a potential type does not implement a custom method.
"""
function get_species end

"""
    get_parameters(p::AbstractPotential)::NamedTuple

Retrieve a `NamedTuple` with the trainable parameters of a potential.

For a `TrainablePotential{P,HP}`, the function automatically extracts the fields with names in the type parameter `P`.
For a `NonTrainablePotential`, the function returns an empty `NamedTuple`.
"""
function get_parameters end
"""
    set_parameters(p::P<:AbstractPotential, parameters::NamedTuple)::P

Generate a new potential of the same type as an exising potential but with modified trainable parameters.

For a `TrainablePotential{P,HP}`, the function copies over all fields and replaces the fields from the type parameter `P`.
For a `NonTrainablePotential`, the function returns the original potential, unchanged.
"""
function set_parameters end
"""
    serialize_parameters(p::AbstractPotential)::Vector{<:AbstractFloat}

Generate a `Vector` form of the trainable parameters of a potential.

For a `TrainablePotential{P,HP}`, the function automatically extracts the fields with names in the type parameter `P` and collects them into a `Vector`.
Note that this behavior should be overloaded for a potential type that has any non-scalar trainable parameters.
For a `NonTrainablePotential`, the function returns an empty `Vector`.
"""
function serialize_parameters end
"""
    deserialize_parameters(p::P<:AbstractPotential, parameters::AbstractVector{<:AbstractFloat})::P

Generate a new potential of the same type as an exising potential but with modified trainable parameters.
This function takes in the trainable parameters as a `Vector` in the same form as produced by `serialize_parameters(p)`.

For a `TrainablePotential{P,HP}`, the function converts the `Vector` of parameters into the type parameter `P` then calls the `set_parameters` implementation.
Note that this behavior should be overloaded for a potential type that has any non-scalar trainable parameters.
For a `NonTrainablePotential`, the function returns the original potential, unchanged.
"""
function deserialize_parameters end

"""
    get_hyperparameters(p::AbstractPotential)::NamedTuple

Retrieve a `NamedTuple` with the trainable hyperparameters of a potential.

For a `TrainablePotential{P,HP}`, the function automatically extracts the fields with names in the type parameter `HP`.
For a `NonTrainablePotential`, the function returns an empty `NamedTuple`.
"""
function get_hyperparameters end
"""
    set_hyperparameters(p::P<:AbstractPotential, hyperparameters::NamedTuple)::P

Generate a new potential of the same type as an exising potential but with modified trainable hyperparameters.

For a `TrainablePotential{P,HP}`, the function copies over all fields and replaces the fields from the type parameter `HP`.
For a `NonTrainablePotential`, the function returns the original potential, unchanged.
"""
function set_hyperparameters end
"""
    serialize_hyperparameters(p::AbstractPotential)::Vector{<:AbstractFloat}

Generate a `Vector` form of the trainable hyperparameters of a potential.

For a `TrainablePotential{P,HP}`, the function automatically extracts the fields with names in the type parameter `HP` and collects them into a `Vector`.
Note that this behavior should be overloaded for a potential type that has any non-scalar trainable hyperparameters.
For a `NonTrainablePotential`, the function returns an empty `Vector`.
"""
function serialize_hyperparameters end
"""
    deserialize_hyperparameters(p::P<:AbstractPotential, hyperparameters::AbstractVector{<:AbstractFloat})::P

Generate a new potential of the same type as an exising potential but with modified trainable hyperparameters.
This function takes in the trainable hyperparameters as a `Vector` in the same form as produced by `serialize_hyperparameters(p)`.

For a `TrainablePotential{P,HP}`, the function converts the `Vector` of hyperparameters into the type parameter `HP` then calls the `set_hyperparameters` implementation.
Note that this behavior should be overloaded for a potential type that has any non-scalar trainable hyperparameters.
For a `NonTrainablePotential`, the function returns the original potential, unchanged.
"""
function deserialize_hyperparameters end


text:
for (T, funcpairs) in (
    (Mag, (
        (:isfinite, :_is_finite),
        (:isinf, :_is_inf),
        (:isspecial, :_is_special),
        (:iszero, :_is_zero),)),
    (Arf, (
        (:isfinite, :_is_finite),
        (:isinf, :_is_inf),
        (:isinteger, :_is_int),
        (:isnan, :_is_nan),
        (:isneginf, :_is_neg_inf),
        (:isnormal, :_is_normal),
        (:isone, :_is_one),
        (:isposinf, :_is_pos_inf),
        (:isspecial, :_is_special),
        (:iszero, :_is_zero),
        )),
    (Arb, (
        (:isexact, :_is_exact),
        (:isfinite, :_is_finite),
        (:isinteger, :_is_int),
        (:isnegative, :_is_negative),
        (:isnonnegative, :_is_nonnegative),
        (:isnonpositive, :_is_nonpositive),
        (:isnonzero, :_is_nonzero),
        (:isone, :_is_one),
        (:ispositive, :_is_positive),
        (:iszero, :_is_zero),
    )),
    (Acb, (
        (:isexact, :_is_exact),
        (:isfinite, :_is_finite),
        (:isinteger, :_is_int),
        (:isone, :_is_one),
        (:isreal, :_is_real),
        (:iszero, :_is_zero),
    ))
    )
    for (jlf, arbsuffix) in funcpairs
        arbf = Symbol(cprefix(T), arbsuffix)
        @eval begin
            function $jlf(x::$T)
                return !iszero(ccall(@libarb($arbf), Cint, (Ref{$T},), x))
            end
        end
    end
end

for ArbT in (Arf, Arb, Acb, Mag)
    arbf = Symbol(cprefix(ArbT), :_equal)
    @eval begin
        function Base.isequal(x::$ArbT, y::$ArbT)
            res = ccall(@libarb($arbf), Cint, (Ref{$ArbT}, Ref{$ArbT}), x, y)
            return !iszero(res)
        end
    end

    ArbT == Mag && continue

    arbf = Symbol(cprefix(ArbT), :_equal_si)
    @eval begin
        function Base.isequal(x::$ArbT, y::Int)
            res = ccall(@libarb($arbf), Cint, (Ref{$ArbT}, Clong), x, y)
            return !iszero(res)
        end
        Base.isequal(y::Int, x::$ArbT) = isequal(x,y)
    end
end

Base.:(==)(x::Arf, y::Arf) = isequal(x, y)
Base.:(==)(x::Arf, y::Int) = isequal(x, y)
Base.:(==)(y::Int, x::Arf) = isequal(x, y)

for (suffix, jltype, ctype) in (
        ("",   Arf, Ref{Arf}),
        (:_si, Int, Clong),
        (:_ui, UInt, Culong),
        (:_d,  Float64, Cdouble)
        )
    arbf = Symbol(:arf_cmp, suffix)
    @eval begin
        function Base.isless(x::Arf, y::$jltype)
            res = ccall(@libarb($arbf), Cint, (Ref{Arf}, $ctype), x, y)
            return res < 0
        end
        Base.:(<)(x::Arf, y::$jltype) = isless(x, y)
        Base.:(<=)(x::Arf, y::$jltype) = isequal(x,y) | isless(x,y)
    end
end

for (ArbT, args) in (
    (Arb, (
        (:(==), :_eq), (:(!=), :_ne),
        (:(<),  :_lt), (:(<=), :_le),
        (:(>),  :_gt), (:(>=), :_ge),
        )),
    (Acb, (
        (:(==), :_eq), (:(!=), :_ne),
        )),
    )
    for (jlf, suffix) in args
        arbf = Symbol(cprefix(ArbT), suffix)
        @eval begin
            function Base.$jlf(x::$ArbT, y::$ArbT)
                res = ccall(@libarb($arbf), Cint, (Ref{$ArbT}, Ref{$ArbT}), x, y)
                return !iszero(res)
            end
        end
    end
end


text:
using Test
using FinancialToolbox,DualNumbers

#Test Parameters
testToll=1e-14;

#print_colored("Starting Standard Test\n",:green)
spot=10;K=10;r=0.02;T=2.0;sigma=0.2;d=0.01;
spotDual=dual(spot,1.0);

toll=1e-8

#EuropeanCall Option
PriceCall=blsprice(spot,K,r,T,sigma,d);
PriceCallBlack=blkprice(spot,K,r,T,sigma);

df(f,x,dx)=(f(x+dx)-f(x))/dx
dx=1e-8;
assert_(value,toll)=@test abs(value)<toll

sigma1=blsimpv(spot, K, r, T, PriceCall, d);
ResDual=blsimpv(spotDual, K, r, T, PriceCall, d);
DerDF_=df(spot->blsimpv(spot, K, r, T, PriceCall, d),spot,dx);

assert_(sigma1-ResDual,toll)
assert_(DerDF_-ResDual.epsilon,toll)

sigma1=blkimpv(spot, K, r, T, PriceCallBlack);
ResDual=blkimpv(spotDual, K, r, T, PriceCallBlack);
DerDF_=df(spot->blkimpv(spot, K, r, T, PriceCallBlack),spot,dx);

assert_(sigma1-ResDual,toll)
assert_(DerDF_-ResDual.epsilon,toll)

print_colored("Starting DualNumbers Test\n",:green)

#Test Parameters
testToll=1e-14;
spot=10;K=10;r=0.02;T=2.0;sigma=0.2;d=0.01;

#EuropeanCall Option
PriceCall=blsprice(spot,K,r,T,sigma,d);
DeltaCall=blsdelta(spot,K,r,T,sigma,d);
ThetaCall=blstheta(spot,K,r,T,sigma,d);
RhoCall=blsrho(spot,K,r,T,sigma,d);
LambdaCall=blslambda(spot,K,r,T,sigma,d);
VannaCall=blsvanna(spot,K,r,T,sigma,d);

#EuropeanPut Option
PricePut=blsprice(spot,K,r,T,sigma,d,false);
DeltaPut=blsdelta(spot,K,r,T,sigma,d,false);
ThetaPut=blstheta(spot,K,r,T,sigma,d,false);
RhoPut=blsrho(spot,K,r,T,sigma,d,false);
LambdaPut=blslambda(spot,K,r,T,sigma,d,false);
VannaPut=blsvanna(spot,K,r,T,sigma,d,false);

#Equals for both Options
Gamma=blsgamma(spot,K,r,T,sigma,d);
Vega=blsvega(spot,K,r,T,sigma,d);

########DUAL NUMBERS
DerToll=1e-13;
#Function definition
#Call
FcallDual(spot)=blsprice(spot,K,r,T,sigma,d);
GcallDual(r)=blsprice(spot,K,r,T,sigma,d);
HcallDual(T)=blsprice(spot,K,r,T,sigma,d);
LcallDual(sigma)=blsprice(spot,K,r,T,sigma,d);
PcallDual(spot)=blsprice(spot,K,r,T,sigma,d)*spot.value/blsprice(spot.value,K,r,T,sigma,d);
VcallDual(sigma)=blsdelta(spot,K,r,T,sigma,d);
#Put
FputDual(spot)=blsprice(spot,K,r,T,sigma,d,false);
GputDual(r)=blsprice(spot,K,r,T,sigma,d,false);
HputDual(T)=blsprice(spot,K,r,T,sigma,d,false);
PputDual(spot)=blsprice(spot,K,r,T,sigma,d,false)*spot.value/blsprice(spot.value,K,r,T,sigma,d,false);
VPutDual(sigma)=blsdelta(spot,K,r,T,sigma,d,false);
#Input
SpotDual=Dual(spot,1.0);
rDual=Dual(r,1.0);
TDual=Dual(T,1.0);
SigmaDual=Dual(sigma,1.0);

#Automatic Differentiation Test
#TEST
print_colored("--- European Call Sensitivities: DualNumbers\n",:yellow)
print_colored("-----Testing Delta\n",:blue);
assert_(FcallDual(SpotDual).epsilon-DeltaCall,DerToll)
print_colored("-----Testing Rho\n",:blue);
assert_(GcallDual(rDual).epsilon-RhoCall,DerToll)
print_colored("-----Testing Theta\n",:blue);
assert_(-HcallDual(TDual).epsilon-ThetaCall,DerToll)
print_colored("-----Testing Lambda\n",:blue);
assert_(PcallDual(SpotDual).epsilon-LambdaCall,DerToll)
print_colored("-----Testing Vanna\n",:blue);
assert_(VcallDual(SigmaDual).epsilon-VannaCall,DerToll)
#TEST
print_colored("--- European Put Sensitivities: DualNumbers\n",:yellow)
print_colored("-----Testing Delta\n",:blue);
assert_(FputDual(SpotDual).epsilon-DeltaPut,DerToll)
print_colored("-----Testing Rho\n",:blue);
assert_(GputDual(rDual).epsilon-RhoPut,DerToll)
print_colored("-----Testing Theta\n",:blue);
assert_(-HputDual(TDual).epsilon-ThetaPut,DerToll)
print_colored("-----Testing Lambda\n",:blue);
assert_(PputDual(SpotDual).epsilon-LambdaPut,DerToll)
print_colored("-----Testing Vanna\n",:blue);
assert_(VPutDual(SigmaDual).epsilon-VannaPut,DerToll)

print_colored("-----Testing Vega\n",:blue);
assert_(LcallDual(SigmaDual).epsilon-Vega,DerToll)
print_colored("Dual Numbers Test Passed\n",:green)
println("")

#TEST OF INPUT VALIDATION
print_colored("Starting Input Validation Test Dual\n",:magenta)
print_colored("----Testing Negative  Spot Price \n",:cyan)
@test_throws(ErrorException, blsprice(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blkprice(-SpotDual,K,r,T,sigma))
@test_throws(ErrorException, blsdelta(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blsgamma(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blstheta(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blsrho(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blspsi(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blslambda(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blsvanna(-SpotDual,K,r,T,sigma,d))
@test_throws(ErrorException, blsvega(-SpotDual,K,r,T,sigma,d))

print_colored("----Testing Negative  Strike Price \n",:cyan)
KK=Dual(K,0);
@test_throws(ErrorException, blsprice(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blkprice(spot,-KK,r,T,sigma))
@test_throws(ErrorException, blsdelta(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blsgamma(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blstheta(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blsrho(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blspsi(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blslambda(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blsvanna(spot,-KK,r,T,sigma,d))
@test_throws(ErrorException, blsvega(spot,-KK,r,T,sigma,d))

print_colored("----Testing Negative  Time to Maturity \n",:cyan)
@test_throws(ErrorException, blsprice(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blkprice(spot,K,r,-TDual,sigma))
@test_throws(ErrorException, blsdelta(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blsgamma(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blstheta(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blsrho(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blspsi(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blslambda(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blsvanna(spot,K,r,-TDual,sigma,d))
@test_throws(ErrorException, blsvega(spot,K,r,-TDual,sigma,d))

print_colored("----Testing Negative  Volatility \n",:cyan)
@test_throws(ErrorException, blsprice(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blkprice(spot,K,r,T,-SigmaDual))
@test_throws(ErrorException, blsdelta(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsgamma(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blstheta(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsrho(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blspsi(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blslambda(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsvanna(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsvega(spot,K,r,T,-SigmaDual,d))
@test_throws(ErrorException, blsimpv(spotDual,K,r,T,-0.2,d))

print_colored("Dual Input Validation Test Passed\n",:magenta)


text:
using ValueOrientedRiskManagementInsurance
using DataFrames
using LinearAlgebra
using Test
import LinearAlgebra.⋅

include("S2Life.jl")
VORMI = ValueOrientedRiskManagementInsurance

println("Start S2Life test")

ins_sum = df_portfolio[1, :ins_sum]

for 𝑡 ∈ 1: T
  @test prob_price[1, :sx] == prob_price[𝑡, :sx]
  @test rfr_price[1] == rfr_price[𝑡]
  @test λ_price[1, :infl] == λ_price[𝑡, :infl]
  @test λ_price[1, :eoy] == λ_price[𝑡, :eoy]
  if 𝑡 > 1
    @test λ_price[𝑡, :boy] == 0
  end
  @test  convert(Array, β[:sx]) == cumsum(fill(β[1, :sx], 5))
  ## following is used in presentation of pricing calculation:
  @test prob_price[𝑡, :qx] ≈ (10 + 𝑡 -1)/10000
end

for 𝑖 ∈ 1:nrow(df_portfolio)
  @test df_portfolio[𝑖, :ins_sum] == ins_sum
end

## Premium ------------------------------------------------------
prob_price[:px] = 1 .- prob_price[:, :qx] .- prob_price[1, :sx]
lx_price_boy =
  convert(Array, cumprod(prob_price[:px]) ./ prob_price[:px])

v_price_eoy = cumprod(1 ./ (1 .+ rfr_price))
v_price_boy = v_price_eoy .* (1 .+ rfr_price)

infl_price_eoy = convert(Array, cumprod(1 .+ λ_price[:infl]))
infl_price_boy =
  infl_price_eoy ./ convert(Array, 1 .+ λ_price[:infl])

prem_price_ratio =
  sum(lx_price_boy .* v_price_boy .*
      λ_price[:, :boy] .* infl_price_boy +
        lx_price_boy .* v_price_eoy .*
      (λ_price[:, :eoy] .* infl_price_eoy +
         prob_price[:, :qx] .* β[:, :qx] +
         prob_price[:, :px] .* β[:, :px])) /
  sum(lx_price_boy .* v_price_boy -
        prob_price[:, :sx] .*
      β[:, :sx] .* lx_price_boy .* v_price_eoy)

for 𝑖 ∈ 1:nrow(df_portfolio)
  for 𝑡 ∈ 1:𝑖
    @test prem_price_ratio ≈ liab_ins.mps[𝑖].β[𝑡, :prem] ./
                             liab_ins.mps[𝑖].β[𝑡, :qx]
  end
end

prem_price = prem_price_ratio * ins_sum

## techn. prov. (pricing) calc ----------------------------------
tp_price = zeros(Float64, T)
for 𝑡 ∈ (T-1):-1:1
  tp_price[𝑡] =
    - prem_price +
    1 / (1 + rfr_price[𝑡 + 1]) *
    (λ_price[𝑡 + 1, :eoy] * infl_price_eoy[𝑡 + 1] *ins_sum +
       prob_price[𝑡+1, :qx] * ins_sum +
       prob_price[𝑡 + 1, :sx] * β[𝑡 + 1, :sx] * prem_price +
       prob_price[𝑡 + 1, :px] * (β[𝑡 + 1, :px] * ins_sum +
                                   tp_price[𝑡 + 1]))
end

for 𝑖 ∈ 1:nrow(df_portfolio)
  t_contract = T + df_portfolio[𝑖, :t_start]
    @test tp_price[T-t_contract+1:T] ≈ liab_ins.mps[𝑖].tpg_price/
                                       df_portfolio[𝑖, :n]
end

## best estimate assumptions ------------------------------------

y_stock =
  vcat(proc_stock.x[1]/proc_stock.x_0 - 1,
       Float64[proc_stock.x[t] / proc_stock.x[t-1] - 1 for t in 2:T])
delta_qx = prob_price[1, :qx] - prob_be[1, :qx]
for 𝑡 ∈ 1:T
  @test y_stock[𝑡] ≈ proc_stock.x[1]-1
end
for 𝑡 ∈ 1:T
  @test prob_be[𝑡, :qx] + delta_qx == prob_price[𝑡, :qx]
end

## Costs ========================================================

## In our example the cost inflation is constant in time
for 𝑑 ∈ 1:nrow(df_portfolio)
  for 𝑥 ∈ cost_infl_be[𝑑]
    @test 𝑥 == cost_infl_be[1][1]
  end
end
## In our example all costs are constant in time
for 𝑡 ∈ 1:T
  if 𝑡 > 1
    @test  λ_be[𝑡, :boy] == 0
  end
  @test λ_be[1, :eoy] == λ_be[𝑡, :eoy]
  @test λ_invest[:IGCash][1, :rel] == λ_invest[:IGCash][𝑡, :rel]
  @test λ_invest[:IGCash][1, :abs] == λ_invest[:IGCash][𝑡, :abs]
  @test λ_invest[:IGStock][1, :rel] == λ_invest[:IGStock][𝑡, :rel]
  @test λ_invest[:IGStock][1, :abs] == λ_invest[:IGStock][𝑡, :abs]
end

## State of the economy =========================================

## notice that indices of tmp_stock, tmp_state are shifted by one
tmp_stock = [cap_mkt.stock.x_0; cap_mkt.stock.x]
state = Float64[(tmp_stock[𝑡+1]/tmp_stock[𝑡]-1-rfr[𝑡])/rfr[𝑡]
                for 𝑡 ∈ 1:T]
tmp_state =
  [cap_mkt.stock.yield_0 / cap_mkt.rfr.yield_0 - 1; state]
state_avg =
  Float64[(tmp_state[𝑡 + 1] + tmp_state[𝑡]) / 2 for 𝑡 ∈ 1:T]

allocation = 0.5 * (1 .- exp.(-max.(0, state_avg)))
## restore initial allocation
allocation[1] = invs.igs[:IGStock].alloc.total[1]

state_orig =
  Float64[VORMI.dynstate(𝑡,cap_mkt) for 𝑡 ∈ 1:T]
state_avg_orig =
  Float64[VORMI.dynstateavg(𝑡,cap_mkt) for 𝑡 ∈ 1:T]
@test state ≈ state_orig
@test state_avg ≈ state_avg_orig
@test allocation ≈ invs.igs[:IGStock].alloc.total

y_invest =
  Float64[allocation[𝑡] * y_stock[𝑡] + (1-allocation[𝑡]) * rfr[𝑡] for 𝑡 ∈ 1:T]

t_bonus_quota = dyn.bonus_factor * (y_invest - rfr_price)

sx_basis = Array{Vector{Float64}}(undef, nrow(df_portfolio))
for 𝑖 ∈ 1:nrow(df_portfolio)
  sx_basis[𝑖] = convert(Array, liab_ins.mps[𝑖].prob[:sx])
end

v_eoy = cumprod(1 ./ (1 .+ rfr))
v_boy = v_eoy .* (1 .+ rfr)

infl_eoy = cumprod(1 .+ cost_infl_be[end])
infl_boy = infl_eoy ./ (1 .+ cost_infl_be[end])

prob = deepcopy(prob_be)
prob[:px] = 1 .- prob[:qx] - prob[:sx]

rfr_cost = rfr - λ_invest[:IGCash][:, :rel]

function tpberec(tp_next, t, τ, prob_sx)
  prob_px = 1 .- prob[:qx] - prob_sx

  - prem_price +
    λ_be[t + τ - t_0 + 1, :boy] * infl_boy[t + 1] * ins_sum +
    1 / (1 + rfr[t + 1]) *
      (λ_be[t + τ - t_0 + 1, :eoy] * infl_eoy[t + 1] * ins_sum +
      prob[t + τ - t_0 + 1, :qx] * ins_sum +
      prob_sx[t + τ - t_0 + 1] * β[t + τ - t_0 + 1, :sx] * prem_price +
      prob_px[t + τ - t_0 + 1] * (β[t + τ - t_0 + 1, :px] * ins_sum + tp_next))
end

τ=1
d=4
t=3
prem_price
mp = deepcopy(liab_ins.mps[d])
fn = df_portfolio[d, :n]
prob_sx = convert(Array, prob[:sx]) * df_portfolio[d, :sx_be_fac]
prob_px = 1 .- prob[:qx] - prob_sx
@test fn * prem_price ≈ mp.β[t+1, :prem]
@test mp.λ[t + 1, :boy] *  mp.λ[t + 1, :cum_infl] /
      (1 + mp.λ[t + 1, :infl]) ≈
      λ_be[t + τ - t_0 + 1, :boy] * infl_eoy[t + 1] * ins_sum
disc = 1/(1 + rfr[t + 1])
@test 1 / (1 + rfr_cost[t + 1]) ≈
      disc/(1 - disc * invs.igs[:IGCash].cost.rel[τ + 1])
@test fn * λ_be[t + τ - t_0 + 1, :eoy] *
      infl_eoy[t + 1] * ins_sum ≈
      mp.λ[t + 1, :eoy] * mp.λ[t + 1, :cum_infl]
@test fn * λ_be[t + τ - t_0 + 1, :eoy] * ins_sum ≈
      mp.λ[t + 1, :eoy]

@test fn * prob[t + τ - t_0 + 1, :qx] * ins_sum  ≈
      mp.prob[t + 1, :qx] * mp.β[t + 1, :qx]
@test fn * prob_sx[t + τ - t_0 + 1] *
      β[t + τ - t_0 + 1, :sx] * prem_price ≈
      mp.prob[t + 1, :sx] * mp.β[t + 1, :sx]
@test fn * prob_px[t + τ - t_0 + 1] *
      (β[t + τ - t_0 + 1, :px] * ins_sum) ≈
      mp.prob[t + 1, :px] * (mp.β[t + 1, :px] + 0)

## Going concern ================================================
tmp_gc = zeros(Float64, T)
for 𝑑 ∈ 1:nrow(df_portfolio)
  lx_boy = df_portfolio[𝑑, :n]
  tmp_gc[1] += lx_boy
  for 𝑡 = 2:nrow(liab_ins.mps[𝑑].prob)
    lx_boy *=
      (1 - liab_ins.mps[𝑑].prob[𝑡 - 1,:qx] -
         liab_ins.mps[𝑑].prob[𝑡 - 1,:sx])
    tmp_gc[𝑡] += lx_boy

  end
end
tmp_gc /= sum(df_portfolio[:n])

@test tmp_gc ≈ liab_ins.gc

tmp_gc_extension = [liab_ins.gc; 0]
tmp_Δgc = zeros(Float64,T)
for 𝑡 ∈ 1:T
  tmp_Δgc[𝑡] = tmp_gc_extension[𝑡+1] - tmp_gc_extension[𝑡]
end
@test liab_ins.Δgc ≈ tmp_Δgc

## Going concern absolute costs ---------------------------------
cost_abs =
  Float64[liab_ins.gc[𝑡] *
            (λ_invest[:IGCash][𝑡, :abs] *
               prod(1 .+ (λ_invest[:IGCash][1:𝑡, :infl_abs])) +
               λ_invest[:IGStock][𝑡, :abs] *
               prod(1 .+ (λ_invest[:IGStock][1:𝑡, :infl_abs])))
          for 𝑡 ∈ 1:T]

@test proj.fixed_cost_gc ≈ cost_abs

## Subordinated debt --------------------------------------------
# cf_liab_other_unscaled =
#   fill(-df_sub_debt[1, :coupon], df_sub_debt[1, :t_mat] - t_0)
# cf_liab_other_unscaled[df_sub_debt[1, :t_mat] - t_0] -=
#   df_sub_debt[1, :nominal]

l_other = deepcopy(liab_other)
VORMI.goingconcern!(l_other, liab_ins.Δgc)
cf_l_other = Array{Vector{Float64}}(undef, T)
for 𝑡 ∈ 1:T
  cf_l_other[𝑡] = zeros(Float64, l_other.subord[𝑡].τ_mat)
  fill!(cf_l_other[𝑡], -l_other.subord[𝑡].coupon)
  cf_l_other[𝑡][l_other.subord[𝑡].τ_mat] -=
    l_other.subord[𝑡].nominal
end

cf_l_other_total = zeros(Float64, T)
for 𝑡 ∈ 1:T
  cf_l_other_total[1:𝑡] += cf_l_other[𝑡]
end

@test sum([cf_l_other[𝑡][end] for 𝑡 ∈ 1:T]) ≈
      -df_sub_debt[1, :coupon] - df_sub_debt[1, :nominal]
for 𝑡 ∈ 2:T
  @test cf_l_other[𝑡][1] / cf_l_other[𝑡][end]  ≈
        df_sub_debt[1, :coupon] /
        (df_sub_debt[1, :coupon] + df_sub_debt[1, :nominal])
end
@test sum(-[cumprod(1 ./ (1 .+ rfr))[1:𝑡] ⋅ cf_l_other[𝑡]
            for 𝑡 ∈ 1:T]) ≈
      proj.val_0[1, :l_other]

for 𝑑 = 1:(T-1)
  @test (-cumprod(1 ./ (1 .+ rfr[𝑑+1:T])) ⋅
        cf_l_other_total[𝑑+1:T]) ≈
        proj.val[𝑑, :l_other]
end

## gc surplus adjustment ----------------------------------------
@test proj.cf[:gc] ≈
      (proj.val_0[1, :invest] -
        proj.val_0[1, :tpg]-
        proj.val_0[1, :l_other]) *
       liab_ins.Δgc

## Technical provisions for guaranteed benefits, t = 0 ==========

tp = Array{Vector{Float64}}(undef, nrow(df_portfolio))
tp_0 = zeros(Float64, nrow(df_portfolio))

for 𝑑 ∈ 1:nrow(df_portfolio)
  # d = 4
  tp[𝑑] = zeros(Float64, T)
  local prob_sx =
    convert(Array, prob[:sx]) * df_portfolio[𝑑, :sx_be_fac]
  local τ = t_0 - df_portfolio[𝑑, :t_start]
  for 𝑡 ∈ (T-1-τ):-1:(1)
    tp[𝑑][𝑡] = tpberec(tp[𝑑][𝑡 + 1], 𝑡, τ, prob_sx)
  end
  tp_0[𝑑] = tpberec(tp[𝑑][1], 0, τ, prob_sx)
end

for 𝑑 ∈ 1:nrow(df_portfolio)
  for 𝑡 ∈ 1:(T + df_portfolio[𝑑, :t_start])
    @test VORMI.tpg(𝑡, cap_mkt.rfr.x, liab_ins.mps[𝑑]) ≈
          df_portfolio[𝑑, :n] * tp[𝑑][𝑡]
  end
end

lx = zeros(Float64, T+1)
tp_all_0 = 0.0
for 𝑑 ∈ 1:nrow(df_portfolio)
  local τ = t_0 - df_portfolio[d, :t_start]
  lx[1] = 1.0
  global tp_all_0 += lx[1] * tp_0[𝑑] * df_portfolio[𝑑, :n]
  for 𝑡 ∈ 1:T
    if 𝑡 + τ - t_0 <= T
      lx[𝑡 + 1] = lx[𝑡] * prob[𝑡 + τ - t_0, :px]
    end
  end
end

@test tp_all_0 ≈ proj.val_0[1, :tpg]

## Cashflows year 1 =============================================

## bi-quotient --------------------------------------------------
@test VORMI.bonusrate(1, y_invest[1], liab_ins.mps[1], dyn) ≈
      t_bonus_quota[1]
@test VORMI.getyield(1, cap_mkt.rfr) ≈ rfr[1]
@test VORMI.getyield(1, cap_mkt.stock) ≈ y_stock[1]

ind_bonus_1 =
  y_stock[1] / (t_bonus_quota[1] + liab_ins.mps[1].rfr_price[1])

@test VORMI.getyield(0, cap_mkt.stock) ≈ cap_mkt.stock.yield_0
@test VORMI.getyield(0, cap_mkt.rfr) ≈ cap_mkt.rfr.yield_0

ind_bonus_hypo =
  cap_mkt.stock.yield_0 /
  (liab_ins.mps[1].bonus_rate_hypo + liab_ins.mps[1].rfr_price_0)

bi_quot_1 =ind_bonus_1 / ind_bonus_hypo

@test bi_quot_1 ≈ VORMI.biquotient(1,
                                   y_invest[1],
                                   cap_mkt, invs,
                                   liab_ins.mps[1],
                                   dyn)

## In the text we assume that b^C,hypo does not depend on C
for 𝑑 ∈ 2:nrow(df_portfolio)
  @test df_portfolio[1,:bonus_rate_hypo] ≈
        df_portfolio[𝑑,:bonus_rate_hypo]
end

## δ_sx_one, qx_one, sx_one, px_one are vectors
## over all model points for time t == 1
δ_sx_one =
  Float64[VORMI.δsx(1, cap_mkt, invs, liab_ins.mps[𝑑], dyn)
   for 𝑑 ∈ 1:nrow(df_portfolio)]

t = 1
for 𝑑 ∈ 1:nrow(df_portfolio)
  if length(sx_basis[𝑑]) >= t
    @test sx_basis[𝑑][t] ≈ liab_ins.mps[𝑑].prob[t, :sx]
  end
end

sx_one = δ_sx_one .* Float64[sx_basis[𝑑][1] for 𝑑 ∈ 1:nrow(df_portfolio)]

qx_one =
  Float64[liab_ins.mps[𝑑].prob[t, :qx]
          for 𝑑 ∈ 1:nrow(df_portfolio)]
px_one = 1 .- qx_one .- sx_one

## Cashflows ----------------------------------------------------

cf_prem_one =
  sum([liab_ins.mps[𝑑].β[1,:prem] for 𝑑 ∈ 1:nrow(df_portfolio)])
cf_λ_boy_one =
  -sum([liab_ins.mps[𝑑].λ[1,:boy] for 𝑑 ∈ 1:nrow(df_portfolio)])
cf_λ_eoy_one =
  -sum([liab_ins.mps[𝑑].λ[1,:eoy] *
          (1 + liab_ins.mps[𝑑].λ[1,:infl])
        for 𝑑 ∈ 1:nrow(df_portfolio)]) -
  invs.igs[:IGCash].mv_0 *
  (1 + (cf_prem_one + cf_λ_boy_one) / invs.mv_0) *
  invs.igs[:IGCash].cost.rel[1] *
  (1 + invs.igs[:IGCash].cost.infl_rel[1]) -
  invs.igs[:IGStock].mv_0 *
  (1 + (cf_prem_one + cf_λ_boy_one) / invs.mv_0) *
  invs.igs[:IGStock].cost.rel[1] *
  (1 + invs.igs[:IGStock].cost.infl_rel[1]) -
  invs.igs[:IGCash].cost.abs[1] *
  (1 + invs.igs[:IGCash].cost.infl_abs[1]) -
  invs.igs[:IGStock].cost.abs[1] *
  (1 + invs.igs[:IGStock].cost.infl_abs[1])
cf_qx_one = -sum([qx_one[𝑑] *liab_ins.mps[𝑑].β[1,:qx]
                  for 𝑑 ∈ 1:nrow(df_portfolio)])
cf_sx_one = -sum([sx_one[𝑑] *liab_ins.mps[𝑑].β[1,:sx]
                  for 𝑑 ∈ 1:nrow(df_portfolio)])
cf_px_one = -sum([px_one[𝑑] *liab_ins.mps[𝑑].β[1,:px]
                  for 𝑑 ∈ 1:nrow(df_portfolio)])
cf_bonus_one =
  -sum([t_bonus_quota[1] * liab_ins.mps[𝑑].tpg_price_0
        for 𝑑 ∈ 1:5])
cf_invest_one =
  (invs.mv_0+cf_prem_one + cf_λ_boy_one) * y_invest[1]

@test cf_prem_one ≈ ins_sum * sum(df_portfolio[:,:n]) *
                    prem_price_ratio
@test cf_qx_one ≈ -sum([qx_one[𝑑] * df_portfolio[𝑑,:n] * ins_sum
                        for 𝑑 ∈ 1:nrow(df_portfolio)])
@test cf_sx_one ≈ -sum([sx_one[𝑑] * sx_fac * df_portfolio[𝑑,:n] *
                        (T - 𝑑 +1) *  ins_sum * prem_price_ratio
                        for 𝑑 ∈ 1:nrow(df_portfolio)])
@test cf_px_one ≈ -px_one[1] * df_portfolio[1,:n] *  ins_sum
@test cf_bonus_one ≈ proj.cf[1, :bonus]
@test cf_prem_one ≈ proj.cf[t, :prem]
@test cf_λ_boy_one ≈ proj.cf[t, :λ_boy]
@test cf_λ_eoy_one ≈ proj.cf[t, :λ_eoy]
@test cf_qx_one ≈ proj.cf[t, :qx]
@test cf_sx_one ≈ proj.cf[t, :sx]
@test cf_px_one ≈ proj.cf[t, :px]
@test cf_invest_one ≈ proj.cf[t, :invest]
@test cf_l_other_total[1] ≈ proj.cf[t, :l_other]

## technical provisions for guaranteed benefits -----------------

## the following is used in the text
@test(abs(δ_sx_one[1]-1) > eps(1.) ? true : false)

probabs = Array{DataFrame}(undef, nrow(df_portfolio))
for 𝑑 ∈ 1:nrow(df_portfolio)
  probabs[𝑑] = deepcopy(prob_be)
  probabs[𝑑][:sx] =
    δ_sx_one[𝑑] * convert(Array, probabs[𝑑][:sx]) *
    df_portfolio[𝑑, :sx_be_fac]
  probabs[𝑑][:px] = 1 .- probabs[𝑑][:qx] - probabs[𝑑][:sx]
end

## We recalculate technical provisions with updated sx
for 𝑑 ∈ 1:nrow(df_portfolio)
  tp[𝑑] = zeros(Float64, T)
  local τ = t_0 - df_portfolio[𝑑, :t_start]
  for 𝑡 ∈ (T-1-τ):-1:(1)
    tp[𝑑][𝑡] = tpberec(tp[𝑑][𝑡 + 1], 𝑡, τ, probabs[𝑑][:sx])
  end
end

tpg_0_1_1 =
  sum([probabs[𝑑][T-𝑑+1, :px] * tp[𝑑][1] *
         df_portfolio[𝑑, :n]
       for 𝑑 ∈ 1:nrow(df_portfolio)])

@test tpg_0_1_1 ≈ proj.val[1,:tpg]
@test proj.val[1,:tpg]-proj.val_0[1,:tpg] ≈ -proj.cf[1,:Δtpg]

# Profit and tax ------------------------------------------------
@test sum(convert(Array, proj.cf[1, [:prem, :λ_boy, :λ_eoy, :qx,
                                     :sx, :px, :invest, :Δtpg, :l_other, :bonus]])) ≈
      proj.cf[1, :profit]

tax_pre = zeros(Float64, T)
tax= zeros(Float64, T)
tax_credit = tax_credit_0
tax_credit_vec = zeros(Float64, T)

for 𝑡 ∈ 1:T
  tax_pre[𝑡] = tax_rate * proj.cf[𝑡,:profit]
  if tax_credit > tax_pre[𝑡]
    tax[𝑡] = 0
    global tax_credit -= tax_pre[𝑡]
  else
    tax[𝑡] = tax_pre[𝑡] - tax_credit
    global tax_credit = 0
  end
  tax_credit_vec[𝑡] = tax_credit
end

@test tax ≈ -proj.cf[:,:tax]

## Dividends  ----------------------------------------------------
surp_quota = zeros(Float64,T)
for 𝑡 ∈ 1:T
  surp_quota[𝑡] =
    proj.val[𝑡, :invest] /
    (proj.val[𝑡, :tpg] + proj.val[𝑡, :l_other]) - 1
end

invest_eoy_prev =
  Float64[𝑡 == 1 ?
            proj.val_0[𝑡, :invest] :
            proj.val[𝑡-1, :invest] for 𝑡 ∈ 1:T]
invest_boy =
  convert(Array,
          invest_eoy_prev + proj.cf[:prem] + proj.cf[:λ_boy])
invest_eoy_pre_divid =
  invest_eoy_prev + proj.cf[:profit] + proj.cf[:tax] -
  proj.cf[:Δtpg] + proj.cf[:gc]

@test invest_boy ≈ invs.mv_boy
@test invest_eoy_pre_divid ≈
      Float64[VORMI.investpredivid(𝑡, invs, proj) for 𝑡 ∈ 1:T]

val_liab = convert(Array, proj.val[:tpg] .- proj.val[:l_other])
q_surp =
  ((invest_eoy_pre_divid .- proj.val[:tpg] .-
    proj.val[:l_other]) ./
  (proj.val[:tpg] + proj.val[:l_other]))

## The following was used in the text:
@test q_surp[1] > dyn.quota_surp
@test  min.(0, convert(Array,
              (1+dyn.quota_surp) * ( proj.val[:tpg] +
                                    proj.val[:l_other])-
                                    invest_eoy_pre_divid)) ≈
      proj.cf[:divid]

## Dividend mechanism works:
for 𝑡 ∈ 1:T
  if (proj.cf[𝑡, :divid] < 0) &
      (proj.val[𝑡, :tpg] + proj.val[𝑡, :l_other] > 0)
    @test dyn.quota_surp ≈
          (proj.val[𝑡, :invest] - proj.val[𝑡, :tpg] -
            proj.val[𝑡, :l_other]) /
          (proj.val[𝑡, :tpg] + proj.val[𝑡, :l_other])
  end
  if proj.cf[t, :divid] ≥ 0
    @test proj.cf[t, :divid] ≈ 0
  end
end

@test proj.val[1, :invest] ≈
      invest_eoy_pre_divid[1] + proj.cf[1, :divid]

## Balance sheet for ≥ 0 ========================================

fdb = zeros(Float64, T)
for 𝑡 ∈ (T-1):-1:1
  fdb[𝑡] =
    (fdb[𝑡 + 1] -  proj.cf[𝑡 + 1, :bonus]) /
    (1 .+ rfr[𝑡 + 1])
end
fdb_0 =
  (fdb[1] -  proj.cf[1, :bonus]) /  (1 .+ rfr[1])

@test -cumprod(1 ./ (1 .+ rfr)) ⋅ proj.cf[:bonus] ≈
      proj.val_0[1,:bonus]
@test fdb_0 ≈ proj.val_0[1,:bonus]
@test fdb ≈ proj.val[:bonus]

balance = vcat(proj.val_0, proj.val)
## here reserves for boni are considered part of capital.
@test balance[:surplus] ≈
      balance[:invest] - balance[:tpg] - balance[:l_other]

## recall that balance[t, :] = proj.val[t-1, :] for t>1
for 𝑡 ∈ 2:(T+1)
  for 𝑤 ∈ [:invest, :tpg, :l_other, :surplus, :bonus, :cost_prov]
     @test balance[𝑡, 𝑤] ≈ proj.val[𝑡-1, 𝑤]
  end
end

## technical provisions for absolute costs & investment costs ---

# provisions absolute costs & investment costs are correct:
for 𝜏 ∈ 1:5
  x = balance[𝜏, :cost_prov]
  for 𝑡 ∈ 𝜏:T
    x *= (1 + rfr[𝑡])
    x -= sum(convert(Array,
                     balance[𝑡,
                             [:tpg, :bonus,
                              :l_other, :cost_prov]])) *
      invs.igs[:IGCash].cost.cum_infl_rel[𝑡] *
      invs.igs[:IGCash].cost.rel[𝑡]
    x -= proj.fixed_cost_gc[𝑡]
  end
  @test x ≈ 0 atol=1.0e-14
end

tpgprev(tp_next, t) =
  (tp_next + cost_abs[t + 1])/(1 + rfr[t + 1])

tp_cost_abs = zeros(Float64, T)
for 𝑡 ∈ (T-1):-1:1
  tp_cost_abs[𝑡] = tpgprev(tp_cost_abs[𝑡 + 1], 𝑡)
end
tp_cost_abs_0 = tpgprev(tp_cost_abs[1], 0)

@test [tp_cost_abs_0; tp_cost_abs] ≈
      Float64[VORMI.tpgfixed(𝑡, cap_mkt.rfr.x[1:liab_ins.dur],
                             proj.fixed_cost_gc)
              for 𝑡 ∈ 0:T]

## S2 Example ###################################################

liabs_mod_0 = balance[1,:tpg] +balance[1,:bonus] + balance[1, :cost_prov]
assets_mod_0 = balance[1,:invest] + proj.tax_credit_0
bof_0 = assets_mod_0 - liabs_mod_0
symb_bal = [:invest, :tpg, :l_other, :surplus, :bonus]
@test convert(Array, balance[1,symb_bal]) ≈
      convert(Array, s2.balance[1, symb_bal])
@test VORMI.bof(s2, :be) ≈ bof_0

## S2 Example Interest ------------------------------------------

ind_mkt = findfirst( (in)([:S2Mkt]), ds2[:mdl])
ind_mkt_int = findfirst( (in)([:S2MktInt]), ds2_mkt[:mdl])
s2_mkt_int = s2.mds[ind_mkt].mds[ind_mkt_int]
rfr_up = VORMI.rfrshock(cap_mkt.rfr.x, s2_mkt_int, :spot_up)
rfr_down = VORMI.rfrshock(cap_mkt.rfr.x, s2_mkt_int, :spot_down)

@test rfr ≈ VORMI.spot2forw(VORMI.forw2spot(rfr))
@test rfr_down ≈ VORMI.spot2forw(VORMI.forw2spot(rfr) .*
                 (1 .+ ds2_mkt_int[:shock][:spot_down][1:T]))
@test rfr_up ≈ VORMI.spot2forw(
                   VORMI.forw2spot(rfr) .+
                  max.(ds2_mkt_int[:spot_up_abs_min],
                      VORMI.forw2spot(rfr) .*
                      ds2_mkt_int[:shock][:spot_up][1:T]))

## S2 Example Equity --------------------------------------------
ind_mkt = findfirst( (in)([:S2Mkt]), ds2[:mdl])
ind_mkt_eq = findfirst( (in)([:S2MktEq]), ds2_mkt[:mdl])
s2_mkt_eq = s2.mds[ind_mkt].mds[ind_mkt_eq]
bal = s2_mkt_eq.balance

@test bal[bal[:scen] .== :type_1,:invest][1] ≈
      (1 + eq_shock[:type_1]) * sum(df_stock[:mv_0]) + sum(df_cash[:mv_0])

## S2 Example Market risk----------------------------------------
ind_mkt = findfirst( (in)([:S2Mkt]), ds2[:mdl])
s2_mkt = s2.mds[ind_mkt]
@test s2_mkt_int.scen_up == false
corr_mkt = s2_mkt.corr_down[1:2,1:2]
scr_mkt_net = [s2_mkt_int.scr[NET], s2_mkt_eq.scr[NET]]
scr_mkt_gross = [s2_mkt_int.scr[GROSS], s2_mkt_eq.scr[GROSS]]
@test sqrt(scr_mkt_net ⋅ (corr_mkt * scr_mkt_net)) ≈
      s2_mkt.scr[NET]
@test sqrt(scr_mkt_gross ⋅ (corr_mkt * scr_mkt_gross)) ≈
      s2_mkt.scr[GROSS]

## Default Risk type 1 -----------------------------------------
ind_def = findfirst( (in)([:S2Def]), ds2[:mdl])
s2_def = s2.mds[ind_def]
inv_len = length(invs.igs[:IGCash].investments)
@test inv_len == 2

accs_mv_0 =
  Float64[invs.igs[:IGCash].investments[𝑖].mv_0
          for 𝑖 ∈ 1:inv_len]
accs_cqs =
  Int[parse(Int, string(string(invs.igs[:IGCash].
                                 investments[𝑖].cqs)[end]))
      for 𝑖 ∈ 1:inv_len]
accs_tlgd =
  Float64[s2_def.mds[1].tlgd[accs_cqs[𝑖]+1] for 𝑖 ∈ 1:inv_len]
accs_slgd =
  Float64[s2_def.mds[1].slgd[accs_cqs[𝑖]+1] for 𝑖 ∈ 1:inv_len]
accs_defu =
  Float64[s2_def.mds[1].u[accs_cqs[𝑖]+1,accs_cqs[𝑗]+1]
          for 𝑖 ∈ 1:inv_len, 𝑗 ∈ 1:inv_len]
accs_defv =
  Float64[s2_def.mds[1].v[accs_cqs[𝑖]+1]
          for 𝑖 ∈ 1:inv_len]

accs_var_t = accs_tlgd ⋅ (accs_defu * accs_tlgd)
accs_var_s = accs_slgd ⋅ accs_defv
accs_var = accs_var_t + accs_var_s
accs_sigma_norm = sqrt(accs_var) / sum(accs_tlgd)
accs_scr_low = s2.mds[ind_def].mds[1].scr_par[:low][1]
accs_scr_low_fac = s2.mds[ind_def].mds[1].scr_par[:low][2]
accs_scr_def = accs_scr_low_fac * sqrt(accs_var)

@test s2_def.scr[NET] ≈ accs_scr_def
@test s2_def.scr[GROSS] ≈ accs_scr_def

## Life mortality -----------------------------------------------
ind_life = findfirst( (in)([:S2Life]), ds2[:mdl])
ind_life_qx = 1
s2_life_qx = s2.mds[ind_life].mds[ind_life_qx]

@test collect(keys(s2_life_qx.shock))[1] == :qx
@test s2_life_qx.shock[:qx] > 0
@test ! s2_life_qx.mp_select[:qx][1]
for 𝑖 ∈ 2:length(s2_life_qx.mp_select[:qx])
  @test s2_life_qx.mp_select[:qx][𝑖]
end

## Life longevity -----------------------------------------------

ind_life_px = 2
s2_life_px = s2.mds[ind_life].mds[ind_life_px]
@test collect(keys(s2_life_px.shock))[1] == :px
@test s2_life_px.shock[:px] < 0
for 𝑖 ∈ 2:length(s2_life_px.mp_select[:px])
  @test !s2_life_px.mp_select[:px][𝑖]
end

## Life surrender -----------------------------------------------

ind_life_sx = 4
s2_life_sx = s2.mds[ind_life].mds[ind_life_sx]

@test (sort(collect(keys(s2_life_sx.shock))) ==
         sort([:sx_down, :sx_up,
               :sx_mass_other, :sx_mass_pension]))

for 𝑚 ∈ 2: length(s2_life_sx.mp_select[:sx_down])
  @test ! s2_life_sx.mp_select[:sx_down][𝑚]
end
for 𝑚 ∈ 2:length(s2_life_sx.mp_select[:sx_up])
  @test s2_life_sx.mp_select[:sx_up][𝑚]
end
for 𝑚 ∈ 2:length(s2_life_sx.mp_select[:sx_mass_other])
  @test s2_life_sx.mp_select[:sx_mass_other][𝑚]
end

## Life cost ----------------------------------------------------

ind_life_cost = 5
s2_life_cost = s2.mds[ind_life].mds[ind_life_cost]
@test collect(keys(s2_life_cost.shock)) == [:cost]

## Life cat -----------------------------------------------------
s2.mds[ind_life].mds
ind_life_cat = 7
s2_life_cat = s2.mds[ind_life].mds[ind_life_cat]
@test collect(keys(s2_life_cat.shock)) == [:cat]
for 𝑚 ∈ 2:length(s2_life_cat.mp_select[:cat])
  @test s2_life_cat.mp_select[:cat][𝑚]
end

## Life aggregation ---------------------------------------------

s2_life = s2.mds[ind_life]

ind_life_risks =
  [ind_life_qx, ind_life_sx, ind_life_cost, ind_life_cat]
life_corr = s2_life.corr[ind_life_risks, ind_life_risks]

scrs_life_net = [s2_life_qx.scr[NET],
                 s2_life_sx.scr[NET],
                 s2_life_cost.scr[NET],
                 s2_life_cat.scr[NET]]
scrs_life = [s2_life_qx.scr[GROSS],
             s2_life_sx.scr[GROSS],
             s2_life_cost.scr[GROSS],
             s2_life_cat.scr[GROSS]]

@test sqrt(scrs_life ⋅ (life_corr * scrs_life)) ≈
      s2_life.scr[GROSS]
@test sqrt(scrs_life_net ⋅ (life_corr * scrs_life_net)) ≈
      s2_life.scr[NET]

## BSCR =========================================================

ind = [ind_mkt, ind_def, ind_life]

bscr_corr = ds2[:corr][ind, ind]
bscrs_gross =
  [s2_mkt.scr[GROSS], s2_def.scr[GROSS], s2_life.scr[GROSS]]

bscrs_net = [s2_mkt.scr[NET], s2_def.scr[NET], s2_life.scr[NET]]

@test sqrt(bscrs_net ⋅ (bscr_corr * bscrs_net)) ≈ s2.bscr[NET]
@test sqrt(bscrs_gross ⋅ (bscr_corr * bscrs_gross)) ≈
      s2.bscr[GROSS]

## operational risk =============================================
@test s2.op.comp_tp ≈
      (proj.val_0[1, :tpg] + proj.val_0[1, :bonus]) *
      s2.op.fac[:tp]

@test s2.op.comp_prem ≈
      s2.op.fac[:prem] * s2.op.prem_earned + s2.op.fac[:prem] *
        max(0, s2.op.prem_earned -
               s2.op.fac[:prem_py] * s2.op.prem_earned_prev)

## Adjustments ==================================================
## Adj technical provisions -------------------------------------
@test -max(0.0, min(s2.bscr[GROSS] - s2.bscr[NET],
                    VORMI.fdb(s2, :be))) ≈
      s2.adj_tp
@test s2.adj_dt == 0
@test s2.liabs_mod ≈
       s2.balance[1,:tpg] + s2.balance[1,:bonus] +  s2.balance[1,:cost_prov]
@test s2.invest_mod ≈ s2.balance[1,:invest]

ds2[:coc]

coc = 0.06
balance = vcat(proj.val_0, proj.val)
length(balance)
tpbe =
  convert(Array,
          balance[:tpg] + balance[:bonus] + balance[:cost_prov])
src_future = (tpbe * s2.scr / tpbe[1])[1:T]
discount = 1 ./ cumprod(1 .+ rfr)
risk_margin = coc * src_future ⋅ discount

@test s2.risk_margin ≈ risk_margin

@test s2.scr_ratio ≈
      (s2.invest_mod - s2.liabs_mod - s2.risk_margin) / s2.scr
@test s2.scr_ratio > 1

#################################################################
println("End S2Life test")

sqrt


text:
# using Base.Test
# using MIToS.SIFTS

print("""

Tests for SIFTS Mappings
========================
""")

print("""

Parse SIFTS
-----------
""")

const sifts_file = joinpath(pwd(), "data", "2vqc.xml.gz")

print("""
all against all
""")

let dbs = [(dbUniProt, "P20220"), (dbPfam, "PF09645"), (dbNCBI, "244589"), (dbPDB, "2vqc"), (dbSCOP, "153426"), (dbPDBe, "2vqc")]
  for (to,toid) in dbs, (from, fromid) in dbs
    map = siftsmapping(sifts_file, from, fromid, to, toid)
    for (k,v) in map
      to == from && @test k == v
    end
  end
end

print("""

siftsmapping
------------
""")

print("""
missings
""")

let map = siftsmapping(sifts_file, dbPDBe, "2vqc", dbPDB, "2vqc", chain="A", missings=false)
  @test_throws KeyError map[9]  # Missing
  @test_throws KeyError map[80] # Missing
  @test_throws KeyError map[1]  # Missing
  @test map[10] == "4"
  @test map[79] == "73"
end

let map = siftsmapping(sifts_file, dbPDBe, "2vqc", dbPDB, "2vqc", chain="A")
  @test map[9] == "3"   # Missing
  @test map[80] == "74" # Missing
  @test map[1] == "-5"  # Missing # Negative Resnum
  @test map[10] == "4"
  @test map[79] == "73"
end

print("""
1SSX => Residues with insert codes: 15A 15B
""")

let map = siftsmapping(joinpath(pwd(), "data", "1ssx.xml.gz"), dbPDBe, "1ssx", dbPDB, "1ssx", chain="A")
  residue_A = map[1]
  residue_B = map[2]
  residue_C = map[3]
  @test residue_A == "15A"
  @test residue_B == "15B"
  @test residue_C == "16"
end

let map = siftsmapping(joinpath(pwd(), "data", "1ssx.xml.gz"), dbPDB, "1ssx", dbUniProt, "P00778", chain="A")
  @test map["15A"] == 200
  @test map["15B"] == 201
  @test map["16"] == 202
end

print("""
1CBN => Multiple InterProt annotations, the last is used.
Identical PDBe ResNum for Residue 22:

        <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="22" dbResName="SER">
          <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="1cbn" dbResNum="22" dbResName="SER" dbChainId="A"/>
          <crossRefDb dbSource="UniProt" dbCoordSys="UniProt" dbAccessionId="P01542" dbResNum="22" dbResName="P"/>
          ....
        <residue dbSource="PDBe" dbCoordSys="PDBe" dbResNum="22" dbResName="PRO">
          <crossRefDb dbSource="PDB" dbCoordSys="PDBresnum" dbAccessionId="1cbn" dbResNum="22" dbResName="PRO" dbChainId="A"/>
          <crossRefDb dbSource="UniProt" dbCoordSys="UniProt" dbAccessionId="P01542" dbResNum="22" dbResName="P"/>
          ...

""")

let map = siftsmapping(joinpath(pwd(), "data", "1cbn.xml.gz"), dbPDBe, "1cbn", dbInterPro, "IPR001010", chain="A")
  @test_throws KeyError map[1] # Without InterPro
  @test map[2] == "2" # Same ResNum for different InterPros
end

let map = read(joinpath(pwd(), "data", "1cbn.xml.gz"), SIFTSXML, chain="A")
  @test length(map[1].InterPro) == 0 # Without InterPro
  @test length(map[2].InterPro) == 4 # Same ResNum for different InterPros
end

print("""
1AS5 => NMR
""")

let map = siftsmapping(joinpath(pwd(), "data", "1as5.xml.gz"), dbPDBe,"1as5", dbUniProt, "P56529", chain="A") # missings=true : NMR there are not missing residues
  @test map[23] == 73
  @test_throws KeyError map[24] # Without UniProt
end

let map = siftsmapping(joinpath(pwd(), "data", "1as5.xml.gz"), dbPDBe,"1as5", dbPDB, "1as5", chain="A")
  @test map[24] == "24"
end

print("""
1DPO => Inserted residues lack insertion letters
Single unnamed chain in 1DPO contains insertions at postions 184 (Gly, Phe), 188 (Gly, Lys), and 221 (Ala, Leu) but no insertion letters.
""")

let map = siftsmapping(joinpath(pwd(), "data", "1dpo.xml.gz"), dbPDBe, "1dpo", dbPDB, "1dpo", chain="A") # Unnamed chain is "A" in SIFTS
  @test map[164] == "184"
  @test map[165] == "184A" # Has insertion code in SIFTS
  @test map[169] == "188"
  @test map[170] == "188A" # Has insertion code in SIFTS
  @test map[198] == "221"
  @test map[199] == "221A" # Has insertion code in SIFTS
end

print("""
1IGY => Insertions have more than one copy of the same amino acid in a single insertion block.
For example, chain B in 1IGY contains a block of four residues inserted at sequence position 82. The block contains Leu-Ser-Ser-Leu.
""")

let map = siftsmapping(joinpath(pwd(), "data", "1igy.xml.gz"), dbPDBe, "1igy", dbCATH, "2.60.40.10", chain="B")
  @test map[82] == "82"
  @test map[83] == "82A"
  @test map[84] == "82B"
  @test map[85] == "82C"
end

print("""
1HAG => Chain E begins with 1H, 1G, 1F, ... 1A, then 1 (in reverse alphabetic order)
""")

let map = siftsmapping(joinpath(pwd(), "data", "1hag.xml.gz"), dbPDBe, "1hag", dbPDB, "1hag", chain="E")
  @test map[1] == "1H"
  @test map[2] == "1G"
  @test map[3] == "1F"
  @test map[4] == "1E"
  @test map[5] == "1D"
  @test map[6] == "1C"
  @test map[7] == "1B"
  @test map[8] == "1A"
  @test map[9] == "1"
end

print("""

Functions for SIFTSResidues
===========================
""")

print("""
1NSA => Contains a single (unnamed) protein chain with sequence 7A-95A that continues 4-308.
""")

print("""
Test findobjects & read
""")

let mapp = read(joinpath(pwd(), "data", "1nsa.xml.gz"), SIFTSXML)
  four = findobjects(mapp, dbPDB, Is(:id, "1nsa"), Is(:number, "4"))[1]
  @test findobjects(mapp, dbPDB, Is(:id, "1nsa"), Is(:number, "95A"))[1] + 1 == four
  @test get(mapp[ findobjects(mapp, dbPDB, Is(:id, "1nsa"), Is(:number, "95A")) ][1].PDB).number == "95A"
end

print("""
Test collectcaptures
""")

let cap = collectcaptures(read(joinpath(pwd(), "data", "1nsa.xml.gz"), SIFTSXML), dbPDB, :number, dbPDB, Is(:id, "1nsa"), Is(:number, "95A"))
  @test get(cap[ [!isnull(x) for x in cap] ][1]) == "95A"
end

print("""
1IAO => Contains in chain B (in this order) 1S, 323P-334P, 6-94, 94A, 95-188, 1T, 2T
""")

print("""
Test collectobjects and isobject
""")

let mapp = read(joinpath(pwd(), "data", "1iao.xml.gz"), SIFTSXML)
  @test collectobjects(mapp, dbPDB, Is(:id, "1iao"), Is(:number, "1S"), Is(:chain, "B"))[1].PDBe.number == 1
  i = findobjects(mapp, dbPDB, Is(:id, "1iao"), Is(:number, "1S"), Is(:chain, "B"))[1]
  @test isobject(mapp[i+2], dbPDB, Is(:id, "1iao"), Is(:chain, "B"), Is(:number, "323P"))
end

print("""

Test download
=============
""")

let pdb = "2vqc"
  filename = downloadsifts(pdb)
  try
    @test_throws ErrorException downloadsifts("2vqc_A")
    @test length(read(filename, SIFTSXML)) == length(read(joinpath(pwd(), "data", "$(pdb).xml.gz"), SIFTSXML))
  finally
    rm(filename)
  end
end


text:
include("prelude.jl")

# Ensure the output buffer can be drained after closing the input and/or after
# the program exists. For this, we use a large buffer in order to delay any
# pending read long enough to trigger potential issues
BUFSIZE = 1024*128

# TODO: both tests have issues

#=
@testset "pty" begin
    interact(`cat`, 1) do proc
        buf = " " ^ BUFSIZE
        print(proc, buf)
        sendeof(proc)
        wait(proc)

        # we intentionally ignore anything after the newline
        tmp = readline(proc)
        @test length(tmp) == BUFSIZE
        @test success(proc)
    end
end

@testset "nopty" begin
    interact(`cat`, 1; pty=false) do proc
        buf = " " ^ BUFSIZE
        print(proc, buf)
        sendeof(proc)
        wait(proc)

        # we intentionally ignore anything after the newline
        tmp = readline(proc)
        @test length(tmp) == BUFSIZE
        @test success(proc)
    end
end
=#


text:
struct Singleton{T}
  value::T
end

Base.show(io::IO, x::Singleton) = print(io, "Singleton(", x.value, ")")

collapse(x::Bernoulli) =
  mean(x) == 1.0 ? Singleton(true) :
  mean(x) == 0.0 ? Singleton(false) :
  x

# TODO make this a distribution, approximate logpdf etc

struct Empirical{T}
  samples::Vector{T}
end

@forward Empirical.samples Statistics.mean, Statistics.var, Statistics.std, Base.length

rounded(x) = @sprintf("%.2g", x)

function Base.show(io::IO, d::Empirical{Bool})
  print(io, "Empirical Bernoulli(p=$(rounded(mean(d)))), N=$(length(d))")
end

function Base.show(io::IO, d::Empirical{<:Real})
  print(io, "Empirical(mean = $(rounded(mean(d))), std = $(rounded(std(d))), N=$(length(d)))")
end


text:
## ROC COMPUTATION

"""
    _idx_unique_sorted(v)

Internal function to return the index of unique elements in `v` under the
assumption that the vector `v` is sorted in decreasing order.
"""
function _idx_unique_sorted(v::Vec{<:Real})
    n    = length(v)
    idx  = ones(Int, n)
    p, h = 1, 1
    cur  = v[1]
    @inbounds while h < n
        h     += 1                  # head position
        cand   = v[h]               # candidate value
        cand   < cur || continue    # is it new? otherwise skip
        p     += 1                  # if new store it
        idx[p] = h
        cur    = cand               # and update the last seen value
    end
    p < n && deleteat!(idx, p+1:n)
    return idx
end

"""
    fprs, tprs, ts = roc_curve(ŷ, y) = roc(ŷ, y)

Return the ROC curve for a two-class probabilistic prediction `ŷ` given the
ground  truth `y`. The true positive rates, false positive rates over a range
of thresholds `ts` are returned. Note that if there are `k` unique scores,
there are correspondingly  `k` thresholds and `k+1` "bins" over which the FPR
and TPR are constant:

* `[0.0 - thresh[1]]`
* `[thresh[1] - thresh[2]]`
* ...
* `[thresh[k] - 1]`

consequently, `tprs` and `fprs` are of length `k+1` if `ts` is of length `k`.

To draw the curve using your favorite plotting backend, do `plot(fprs, tprs)`.
"""
function roc_curve(ŷm, ym)
    ŷ, y    = skipinvalid(ŷm, ym)
    n       = length(y)
    lab_pos = levels(y)[2]
    scores  = pdf.(ŷ, lab_pos)
    ranking = sortperm(scores, rev=true)

    scores_sort = scores[ranking]
    y_sort_bin  = (y[ranking] .== lab_pos)

    idx_unique = _idx_unique_sorted(scores_sort)
    thresholds = scores_sort[idx_unique]

    # detailed computations with example:
    # y = [  1   0   0   1   0   0   1]
    # s = [0.5 0.5 0.2 0.2 0.1 0.1 0.1] thresh are 0.5 0.2, 0.1 // idx [1, 3, 5]
    # ŷ = [  0   0   0   0   0   0   0] (0.5 - 1.0] # no pos pred
    # ŷ = [  1   1   0   0   0   0   0] (0.2 - 0.5] # 2 pos pred
    # ŷ = [  1   1   1   1   0   0   0] (0.1 - 0.2] # 4 pos pred
    # ŷ = [  1   1   1   1   1   1   1] [0.0 - 0.1] # all pos pre

    idx_unique_2 = idx_unique[2:end]   # [3, 5]
    n_ŷ_pos      = idx_unique_2 .- 1   # [2, 4] implicit [0, 2, 4, 7]

    cs   = cumsum(y_sort_bin)          # [1, 1, 1, 2, 2, 2, 3]
    n_tp = cs[n_ŷ_pos]                 # [1, 2] implicit [0, 1, 2, 3]
    n_fp = n_ŷ_pos .- n_tp             # [1, 2] implicit [0, 1, 2, 4]

    # add end points
    P = sum(y_sort_bin) # total number of true positives
    N = n - P           # total number of true negatives

    n_tp = [0, n_tp..., P] # [0, 1, 2, 3]
    n_fp = [0, n_fp..., N] # [0, 1, 2, 4]

    tprs = n_tp ./ P  # [0/3, 1/3, 2/3, 1]
    fprs = n_fp ./ N  # [0/4, 1/4, 2/4, 1]

    return fprs, tprs, thresholds
end

const roc = roc_curve


text:
using DynamicBounds, Plots
pyplot()

x0(p) = [34.0; 20.0; 0.0; 0.0; 16.0; 0.0]
function f!(du, u, p, t)
    du[1] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[6]*u[6]
    du[2] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[3]*u[3]
    du[3] =  p[1]*u[1]*u[2] - p[2]*u[3] - p[3]*u[3]
    du[4] =  p[3]*u[3] - p[4]*u[4]*u[5] + p[5]*u[6]
    du[5] = -p[4]*u[4]*u[5] + p[5]*u[6] + p[6]*u[6]
    du[6] =  p[4]*u[4]*u[5] - p[5]*u[6] - p[6]*u[6]
    return
end

function add_trajectory(plt, integrator, pL, pU)
    ratio = rand(6)
    pstar = pL.*ratio .+ pU.*(1.0 .- ratio)
    setall!(integrator, ParameterValue(), pstar)
    integrate!(integrator)
    plot!(plt, integrator.local_problem_storage.integrator_t,
    integrator.local_problem_storage.pode_x[1,:], label="", linecolor = :green,
    markershape = :+, markercolor = :green, linestyle = :dash, markersize = 2, lw=0.75)
    return plt
end

# Example from Improved relaxations for the parametric solutionsof ODEs using differential inequalities
function generate_plot()
    tspan = (0.0,18.0e-5*250)
    pL = [0.1; 0.033; 16.0; 5.0; 0.5; 0.3]
    pU = 10.0*pL

    prob = ODERelaxProb(f!, tspan, x0, pL, pU)

    A = [0.0 -1.0 -1.0  0.0  0.0  0.0;
         0.0  0.0  0.0  0.0 -1.0 -1.0;
         1.0 -1.0 0.0 1.0 -1.0 0.0]
    b = [-20.0; -16.0; -2.0]
    set!(prob, PolyhedralConstraint(A, b))

    xL = zeros(6)
    xU = [34.0; 20.0; 20.0; 34.0; 16.0; 16.0]
    set!(prob, ConstantStateBounds(xL,xU))

    Scott2013(prob, calculate_relax = false, calculate_subgradient = false)
    integrator = Scott2013(prob, calculate_relax = false, calculate_subgradient = false)

    relax!(integrator)

    plt = plot(integrator.relax_t, integrator.relax_lo[1,:], label="Interval Bounds", linecolor = :darkblue, linestyle = :dash, lw=1.5, legend=:bottomleft)
    plot!(plt, integrator.relax_t, integrator.relax_hi[1,:], label="", linecolor = :darkblue, linestyle = :dash, lw=1.5)
    ylabel!("x[1] (M)")
    xlabel!("Time (seconds)")

    ratio = rand(6)
    pstar = pL.*ratio .+ pU.*(1.0 .- ratio)
    setall!(integrator, ParameterValue(), pstar)
    integrate!(integrator)
    plot!(plt, integrator.local_problem_storage.integrator_t,
    integrator.local_problem_storage.pode_x[1,:], label="Trajectories", linecolor = :green,
    markershape = :+, markercolor = :green, linestyle = :dash, markersize = 2, lw=0.75)

    for i in 1:39
        add_trajectory(plt, integrator, pL, pU)
    end
    plt
end

plt = generate_plot()
Plots.eps(plt, "C:\\Users\\wilhe\\Desktop\\Package Development Work\\sample_plot")

#=
function add_trajectory(plt, integrator, pL, pU)
    ratio = rand(6)
    pstar = pL.*ratio .+ pU.*(1.0 .- ratio)
    setall!(integrator, ParameterValue(), pstar)
    integrate!(integrator)
    plot!(plt, integrator.local_problem_storage.integrator_t,
    integrator.local_problem_storage.pode_x[1,:], label="", lw=1)
    return plt
end

plt = add_trajectory(plt, integrator, pL, pU)
plt = add_trajectory(plt, integrator, pL, pU)
plot(plt)
show(plt)
=#

#invariant = DenseLinearInvariant(A, b)
#xnew = Interval.(xL,xU)
#xprior = Interval.(xL,xU)

#=
using Plots
plot([integrator.local_problem_storage.pode_x[1,:],
      integrator.local_problem_storage.pode_x[2,:],
      integrator.local_problem_storage.pode_x[3,:],
      integrator.local_problem_storage.pode_x[4,:],
      integrator.local_problem_storage.pode_x[5,:],
      integrator.local_problem_storage.pode_x[6,:]])
=#

#using Profile, ProfileView
#test_func(x::Int) = relax!(integrator)
#@profview test_func(1)
#@profview test_func(2)
#=
t = integrator.relax_t
println("length(t): $(length(t))")

#cv1 = integrator.relax_cv[1,:]
#cc1 = integrator.relax_cc[1,:]
lo1 = integrator.relax_lo[1,:]
hi1 = integrator.relax_hi[1,:]
val1 = integrator.local_problem_storage.pode_x[1,:]

#cv2 = integrator.relax_cv[2,:]
#cc2 = integrator.relax_cc[2,:]
lo2 = integrator.relax_lo[2,:]
hi2 = integrator.relax_hi[2,:]
val2 = integrator.local_problem_storage.pode_x[2,:]

p1 = Plots.plot(t, lo1, label="lo[1,:]", title="Dim1")
#Plots.plot!(t, cc1, label="cc[1,:]")
#Plots.plot!(t, cv1, label="cv[1,:]")
Plots.plot!(t, hi1, label="hi[1,:]")
Plots.plot!(integrator.local_problem_storage.integrator_t, val1, label="val[1,:]")

prob = integrator.relax_ode_prob
#@code_warntype prob.f.f(ScottA2013a_pL, ScottA2013a_pL, ScottA2013a_pL, ScottA2013a_tf)
=#


text:
E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force

problems = Any[PointLoadCantilever(Val{:Linear}, (80, 20, 20), (1.0, 1.0, 1.0), E, v, f),
            PointLoadCantilever(Val{:Linear}, (160, 40), (1.0, 1.0), E, v, f),
            HalfMBB(Val{:Linear}, (60, 20), (1.0, 1.0), E, v, f),
            LBeam(Val{:Linear}, Float64),
            TieBeam(Val{:Quadratic}, Float64)]
problem_names = ["cantilever beam", "half MBB beam", "L-beam", "tie-beam"]
approx_objvals = [330.0, 175.0, 65.0, 1413.0]

@testset "Continuation SIMP - $(problem_names[i])" for i in 1:4
    # Define the problem
    problem = problems[i]
    # Parameter settings
    V = 0.5 # volume fraction
    xmin = 0.001 # minimum density
    steps = 40 # maximum number of penalty steps, delta_p0 = 0.1
    reuse = true # adaptive penalty flag
    convcriteria = MMA.KKTCriteria()
    #penalty = TopOpt.PowerPenalty(1.0)
    penalty = TopOpt.PowerPenalty(1.0)
    pcont = Continuation(penalty, steps = steps, xmin = xmin, pmax = 5.0)

    mma_options = options = MMA.Options(maxiter=1000)
    if convcriteria isa MMA.KKTCriteria
        maxtol = 0.1 # maximum tolerance
        mintol = 0.001 # minimum tolerance

        b = log(mintol / maxtol) / steps
        a = maxtol / exp(b)
        mma_options_gen = TopOpt.MMAOptionsGen(steps = steps, initial_options = mma_options, kkttol_gen = ExponentialContinuation(a, b, 0.0, steps + 1, mintol))
    else
        maxtol = 0.01 # maximum tolerance
        mintol = 0.0001 # minimum tolerance

        b = log(mintol / maxtol) / steps
        a = maxtol / exp(b)
        mma_options_gen = TopOpt.MMAOptionsGen(steps = steps, initial_options = mma_options, ftol_gen = ExponentialContinuation(a, b, 0.0, steps + 1, mintol))
    end
    csimp_options = TopOpt.CSIMPOptions(steps = steps,
                                        options_gen = mma_options_gen,
                                        p_gen = pcont,
                                        reuse = reuse
                                        )

    # Define a finite element solver
    solver = FEASolver(Displacement, Direct, problem, xmin = xmin, penalty = penalty)
    # Define compliance objective
    filtering = problem isa TopOptProblems.TieBeam ? false : true
    obj = Objective(ComplianceFunction(problem, solver, filtering = filtering,
        rmin = 3.0, tracing = true, logarithm = false))
    #cu_obj = TopOpt.cu(obj)
    # Define volume constraint
    constr = Constraint(VolumeFunction(problem, solver), V)
    # Define subproblem optimizer
    #optimizer = MMAOptimizer{CPU}(cu_obj, constr, MMA.MMA87(),
    #    ConjugateGradient(), maxiter=1000); optimizer.obj.fevals = 0
    optimizer = MMAOptimizer{CPU}(obj, constr, MMA.MMA87(),
        ConjugateGradient(), options = mma_options, convcriteria = convcriteria); optimizer.obj.f.fevals = 0

    # Define continuation SIMP optimizer
    simp = SIMP(optimizer, penalty.p)

    cont_simp = ContinuationSIMP(simp, steps, csimp_options)

    # Solve
    x0 = fill(1.0, length(solver.vars))
    result = cont_simp(x0)

    @test round(result.objval, digits=0) == approx_objvals[i]
end

@testset "Continuation SIMP 2 - $(problem_names[i])" for i in 1:4
    # Define the problem
    problem = problems[i]
    # Parameter settings
    V = 0.5 # volume fraction
    xmin = 0.001 # minimum density
    maxtol = 0.01 # maximum tolerance
    mintol = 0.0001 # minimum tolerance
    steps = 40 # maximum number of penalty steps, delta_p0 = 0.1
    reuse = true # adaptive penalty flag

    # Define a finite element solver
    solver = FEASolver(Displacement, CG, MatrixFree, problem, xmin = xmin,
        penalty = TopOpt.PowerPenalty(1.0))
    # Define volume constraint
    obj = Objective(VolumeFunction(problem, solver))
    # Define compliance objective
    filtering = problem isa TopOptProblems.TieBeam ? false : true
    constr = Constraint(ComplianceFunction(problem, solver, filtering = filtering,
        rmin = 3.0, tracing = true, logarithm = false), approx_objvals[i])
    cu_constr = TopOpt.cu(constr)
    # Define subproblem optimizer

    mma_options = options = MMA.Options(maxiter=1000)
    optimizer = MMAOptimizer{CPU}(obj, constr, MMA.MMA87(),
        ConjugateGradient(), options = mma_options); optimizer.obj.f.fevals = 0

    # Define continuation SIMP optimizer
    simp = SIMP(optimizer, 1.0)

    b = log(mintol / maxtol) / steps
    a = maxtol / exp(b)
    mma_options_gen = TopOpt.MMAOptionsGen(steps = steps, initial_options = mma_options, ftol_gen = ExponentialContinuation(a, b, 0.0, steps + 1, mintol))
    csimp_options = TopOpt.CSIMPOptions(steps = steps, options_gen = mma_options_gen, pstart = 1.0, pfinish = 5.0)

    cont_simp = ContinuationSIMP(simp, steps, csimp_options)

    # Solve
    x0 = fill(1.0, length(solver.vars))
    result = cont_simp(x0)

    @test round(result.objval, digits=0) == approx_objvals[i]
end


text:
# Submodule for imaging functions

#######################################
# Function to resize an array representation
# of an image to a certain number of pixeles
function resize_image(f, N)
  if size(f,1) > size(f,2)
		f = reshape(f,size(f,2),size(f,1),size(f,3));
  end
	P = size(f,1);
	# add 1 pixel boundary
	g = f;
	g = cat(2, g, reshape(g[:,1,:],size(g,1),1,size(g,3)));
    g = cat(1, g, reshape(g[1,:,:],1,size(g,2),size(g,3)));
	# interpolate
	t = linspace(1,P,N);
	ti = round(Int64,floor(t)) ; tj = round(Int64,ceil(t));
	fi = round(Int64,t-floor(t)); fj = 1-fi;
	h = zeros(N,N,size(f,3));
	for s in 1:size(f,3)
	    h[:,:,s] = g[ti,ti,s] .* (fj*fj') + g[tj,tj,s] .* (fi*fi') + g[ti,tj,s] .* (fj*fi') + g[tj,ti,s] .* (fi*fj');
	end
	return h;
end

#######################################
# Function that load an image in a local path
# and resize it to a certain number of pixeles
"""
...
load_image(path::string,pixels::int) load an image in a local path
with N=nxn pixeles
...
"""
function load_image(name,N)
	g = PyPlot.imread(name);
	g = resize_image(g, N);
	return g;
end #load_image

#######################################
# Function that show the image associated
# to a nxn array
"""
...
imageplot(g::Array,title::str="", a=-1, b=-1, c=-1) show the image
associated to a nxn array, with title coordinates in multiplo (a,b,c)
...

"""
function imageplot(g, str="", a=-1, b=-1, c=-1)
	if size(g,3)==1
		# special for b&w images.
		g = g[:,:,1];
	end
	if c>0
		subplot(a,b,c);
	end
	imshow(g, interpolation="nearest")
	set_cmap("gray")
	axis("off")
	if str!=""
		title(str);
	end
end #imageplot

#####################################
# Function that rescales date in a,b
"""
...
rescale(x::Array,a=0,b=1) rescales data in
[a,b]
...
"""
function rescale(x,a=0,b=1)
	m = minimum(x[:]);
	M = maximum(x[:]);

	if M-m<1e-10
	    y = x;
	else
	    y = (b-a) * (x-m)/(M-m) + a;
	end
	return y;

end ## rescale

#####################################################
# Function that rescales a wavelet transform in an array A
function rescaleWav(A)
    v = maximum( abs(A[:]) );
    B = copy(A)
    if v > 0
        B = .5 + .5 * A / v;
    end
    return B;
end

#####################################################
# Function that makes a plot of wavelet coefficients
"""
...
plot_wavelet(fW::Array, Jmin=0) plot of the wavelet
transform with the transformed image data fW an the
minimum scale
...
"""
function plot_wavelet(fW, Jmin=0)
    n = size(fW,1);
    Jmax = round(Int64,log2(n)) - 1;
    U = copy(fW);
    for j = Jmax:-1:Jmin
    	L = 1:2^j; H = (2^j)+1:2^(j+1);
        U[L,H] = rescaleWav( U[L,H] );
        U[H,L] = rescaleWav( U[H,L] );
        U[H,H] = rescaleWav( U[H,H] );
    end
    # coarse scale
    L = 1:2^Jmin;
    U[L,L] = rescale(U[L,L]);
    # plot underlying image
    imageplot(U);
    # display crosses
    for j=Jmax:-1:Jmin
        plot([0, 2^(j+1)], [2^j, 2^j], "r");
        plot([2^j, 2^j], [0, 2^(j + 1)], "r");
    end
    # display box
    plot([0, n], [0, 0], "r");
    plot([0, n], [n, n], "r");
    plot([0, 0], [0, n], "r");
    plot([n, n], [0, n], "r");
    return U;

end # plot_wavelet