text:
using PaddedMatrices: DynamicPaddedVector, DynamicPaddedMatrix, AbstractPaddedMatrix, AbstractDynamicPaddedMatrix
using PaddedMatrices: StackPointer, AbstractMutableFixedSizeArray

abstract type AbstractFixedSizeCovarianceMatrix{M,T,R,L} <: PaddedMatrices.AbstractMutableFixedSizeMatrix{M,M,T,R,L} end
mutable struct FixedSizeCovarianceMatrix{M,T,R,L} <: AbstractFixedSizeCovarianceMatrix{M,T,R,L}
    data::NTuple{L,T}
    FixedSizeCovarianceMatrix{M,T,R,L}(::UndefInitializer) where {M,T,R,L} = new{M,T,R,L}()
end
struct PtrFixedSizeCovarianceMatrix{M,T,R,L} <: AbstractFixedSizeCovarianceMatrix{M,T,R,L}
    ptr::Ptr{T}
end
function padl(M, T = Float64)
    Wm1 = VectorizationBase.pick_vector_width(M,T) - 1
    R = (M + Wm1) & ~Wm1
    R, R*M
end
@generated function PtrFixedSizeCovarianceMatrix{M,T}(sp::StackPointer, ::UndefInitializer = undef) where {M,T}
    R = M
    L = R * M
    :(sp + $(sizeof(T)*L), PtrFixedSizeCovarianceMatrix{$M,$T,$R,$L}(pointer(sp, $T)))
end
@generated function FixedSizeCovarianceMatrix{M,T}(::UndefInitializer) where {M,T}
    R = M
    L = R * M
    :(FixedSizeCovarianceMatrix{$M,$T,$R,$L}(undef))
end
@inline FixedSizeCovarianceMatrix(::Val{M}, ::Type{T}) where {M,T} = FixedSizeCovarianceMatrix{M,T}(undef)
@inline FixedSizeCovarianceMatrix(sp::StackPointer, ::Val{M}, ::Type{T}) where {M,T} = PtrFixedSizeCovarianceMatrix{M,T}(sp)

#DynamicCovarianceMatrix{T}(sp::StackPointer, ::UndefInitializer, N) where {T} = DynamicCovarianceMatrix{T}(sp, N)

Base.pointer(A::PtrFixedSizeCovarianceMatrix) = A.ptr
Base.unsafe_convert(::Type{Ptr{T}}, A::PtrFixedSizeCovarianceMatrix{M,T}) where {M,T} = A.ptr
Base.pointer(A::FixedSizeCovarianceMatrix{M,T}) where {M,T} = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
PtrFixedSizeCovarianceMatrix(S::FixedSizeCovarianceMatrix{M,T,R,L}) where {M,T,R,L} = PtrFixedSizeCovarianceMatrix{M,T,R,L}(pointer(S))
Base.unsafe_convert(::Type{Ptr{T}}, A::FixedSizeCovarianceMatrix{M,T}) where {M,T} = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
LinearAlgebra.checksquare(A::AbstractFixedSizeCovarianceMatrix{M}) where {M} = M
@inbounds Base.stride(A::AbstractFixedSizeCovarianceMatrix{M,T,R}, i::Integer) where {M,T,R} = (1, R)[i]
function Base.zeros(::Type{<:FixedSizeCovarianceMatrix{M,T}}) where {M,T}
    A = FixedSizeCovarianceMatrix{M,T}(undef)
    @inbounds for i ∈ 1:PaddedMatrices.full_length(A)
        A[i] = zero(T)
    end
    A
end

struct DynamicCovarianceMatrix{T,ADT <: AbstractDynamicPaddedMatrix{T}} <: AbstractDynamicPaddedMatrix{T} #L
    data::ADT
end
function DynamicCovarianceMatrix{T}(::UndefInitializer, N) where {T}
    DynamicCovarianceMatrix(
        DynamicPaddedMatrix{T}(undef, (N,N))
    )
end
function DynamicCovarianceMatrix{T}(sp::StackPointer, N) where {T}
    sp, data = DynamicPaddedMatrix{T}(sp, (N,N))
    sp, DynamicCovarianceMatrix( data )
end
DynamicCovarianceMatrix{T}(sp::StackPointer, ::UndefInitializer, N) where {T} = DynamicCovarianceMatrix{T}(sp, N)
Base.pointer(A::DynamicCovarianceMatrix) = pointer(A.data)
Base.unsafe_convert(::Type{Ptr{T}}, A::DynamicCovarianceMatrix{T}) where {T} = pointer(A.data)
LinearAlgebra.checksquare(A::DynamicCovarianceMatrix) = size(A.data,2)
Base.stride(A::DynamicCovarianceMatrix, i) = stide(A.data, i)
Base.size(A::DynamicCovarianceMatrix) = (m = size(A.data,2); (m,m))
Base.getindex(A::DynamicCovarianceMatrix, I...) = A.data[I...]
const AbstractCovarianceMatrix{T} = Union{AbstractFixedSizeCovarianceMatrix{M,T} where M,DynamicCovarianceMatrix{T}}

struct MissingDataVector{T,S,M,R,L}
    indices::FixedSizeVector{S,Int,S} # note, indices start from 0
    ∂Σ::FixedSizeCovarianceMatrix{M,T,R,L} # ∂Σ, rather than stack pointer, so that we can guarantee elements are 0
end
struct MissingDataVectorAdjoint{T,S,M,R,L}#,VT}
    # mdv::MissingDataVector{T,S,M,R,L}#,VT}
    indices::PtrVector{S,Int,S,false}
    ∂Σ::PtrFixedSizeCovarianceMatrix{M,T,R,L}
end
function MissingDataVectorAdjoint(mdv::MissingDataVector{T,S,M,R,L}) where {T,S,M,R,L}
    MissingDataVectorAdjoint{T,S,M,R,L}(
        PtrVector{S,Int,S}(pointer(mdv.indices)),
        PtrFixedSizeCovarianceMatrix(mdv.∂Σ),
    )
end
function MissingDataVector{T}(bitmask) where {T}
    MissingDataVector{T}(PaddedMatrices.FixedSizeVector, bitmask)
end
## Not Type Stable
function MissingDataVector{T}(::Type{<:PaddedMatrices.AbstractFixedSizeVector}, bitmask) where {T}
    Σ = zeros(FixedSizeCovarianceMatrix{length(bitmask),T})
    inds = findall(i -> i == 1, bitmask)
    indices = PaddedMatrices.FixedSizeVector{length(inds),Int}(inds)#::PaddedMatrices.FixedSizeVector
    MissingDataVector(
        indices, Σ
    )
end
Base.size(::MissingDataVectorAdjoint{T,S,M}) where {T,S,M} = (S,M)

#They assume lower triangular symetric
function mask!(
    B::Union{DynamicCovarianceMatrix{T,ADT1},Symmetric{T,ADT1}},
    A::Union{DynamicCovarianceMatrix{T,ADT2},Symmetric{T,ADT2}},
    inds::AbstractVector{<:Integer}
) where {T, ADT1 <: AbstractDynamicPaddedMatrix{T}, ADT2 <: AbstractDynamicPaddedMatrix{T}}
    N = length(inds)
    @inbounds for (ca,cb) ∈ enumerate(inds), ra ∈ ca:N
        rb = inds[ra]
        B[rb,cb] = A[ra,ca]
    end
end
function mask!(
    B::AbstractFixedSizeCovarianceMatrix,
    A::AbstractFixedSizeCovarianceMatrix,
    inds::AbstractVector{<:Integer}
)
    N = length(inds)
    @inbounds for (ca,cb) ∈ enumerate(inds), ra ∈ ca:N
        rb = inds[ra]
        B[rb,cb] = A[ra,ca]
    end
end

function subset!(B::AbstractPaddedMatrix, A::AbstractPaddedMatrix, mdr::MissingDataVector, mdc::MissingDataVector)
    indsr = mdr.indices
    indsc = mdc.indices
    @boundscheck begin
        size(B) == (length(indsr), length(indsc)) || PaddedMatrices.ThrowBoundsError()
    end
    @inbounds for ics ∈ eachindex(indsc), irs ∈ eachindex(indsr)
        B[irs,ics] = A[indsr[irs],indsc[ics]]
    end
end
function subset!(
    B::Union{DynamicCovarianceMatrix{T,ADT1},Symmetric{T,ADT1}},
    A::Union{DynamicCovarianceMatrix{T,ADT2},Symmetric{T,ADT2}},
    mdv::MissingDataVector
    #mrow::BitVector, mcol::Vector{<:Integer}
) where {T, ADT1 <: AbstractDynamicPaddedMatrix{T}, ADT2 <: AbstractDynamicPaddedMatrix{T}}
    inds = mdv.indices
    N = length(inds)
    @boundscheck begin
        size(B) == (N, N) || PaddedMatrices.ThrowBoundsError()
    end
    @inbounds for ics ∈ eachindex(inds), irs ∈ ics:N #eachindex(indsr)
        B[irs,ics] = A[inds[irs],inds[ics]]
    end
end
function subset!(
    B::AbstractFixedSizeCovarianceMatrix{S,T},
    A::AbstractFixedSizeCovarianceMatrix{L,T},
    mdv::MissingDataVector{T,S,L}
) where {T, S, L}
    inds = mdv.indices
    @inbounds for ics ∈ 1:S, irs ∈ ics:S
        B[irs,ics] = A[inds[irs],inds[ics]]
    end
end

function Base.getindex( A::DynamicCovarianceMatrix{T}, mdv::MissingDataVector ) where {T}
    B = DynamicCovarianceMatrix{T}( undef, length(mdv.indices) )
    @inbounds subset!( B, A, mdv )
    B
end
function Base.getindex( sp::StackPointer, A::DynamicCovarianceMatrix{T}, mdv::MissingDataVector ) where {T}
    sp, B = DynamicCovarianceMatrix{T}( sp, length(mdv.indices) )
    @inbounds subset!( B, A, mdv )
    sp, B
end
function PaddedMatrices.∂getindex( A::DynamicCovarianceMatrix{T}, mdv::MissingDataVector ) where {T}
    B = DynamicCovarianceMatrix{T}( undef, length(mdv.indices) )
    @inbounds subset!( B, A, mdv )
    B, MissingDataVectorAdjoint(mdv)
end
function PaddedMatrices.∂getindex( sp::StackPointer, A::DynamicCovarianceMatrix{T}, mdv::MissingDataVector ) where {T}
    sp, B = DynamicCovarianceMatrix{T}( sp, length(mdv.indices) )
    @inbounds subset!( B, A, mdv )
    sp, (B, MissingDataVectorAdjoint(mdv))
end

function Base.getindex( A::AbstractFixedSizeCovarianceMatrix{L,T}, mdv::MissingDataVector{T,S,L} ) where {T,S,L}
    B = FixedSizeCovarianceMatrix{S,T}(undef)
    @inbounds subset!( B, A, mdv )
    B
end
function Base.getindex( sp::StackPointer, A::AbstractFixedSizeCovarianceMatrix{L,T}, mdv::MissingDataVector{T,S,L} ) where {T,S,L}
    sp, B = PtrFixedSizeCovarianceMatrix{S,T}( sp )
    @inbounds subset!( B, A, mdv )
    sp, B
end
function PaddedMatrices.∂getindex( A::AbstractFixedSizeCovarianceMatrix{L,T}, mdv::MissingDataVector{T,S,L} ) where {T,S,L}
    B = FixedSizeCovarianceMatrix{S,T}(undef)
    @inbounds subset!( B, A, mdv )
    B, MissingDataVectorAdjoint(mdv)
end
function PaddedMatrices.∂getindex( sp::StackPointer, A::AbstractFixedSizeCovarianceMatrix{L,T}, mdv::MissingDataVector{T,S,L} ) where {T,S,L}
    sp, B = PtrFixedSizeCovarianceMatrix{S,T}( sp )
    @inbounds subset!( B, A, mdv )
    sp, (B, MissingDataVectorAdjoint(mdv))
end

function Base.getindex(
    a::AbstractMutableFixedSizeVector{M,T},
    mdv::MissingDataVector{T,S,M}
) where {T,S,M}
    inds = mdv.indices
    b = FixedSizeVector{S,T}(undef)
    @inbounds for s ∈ 1:S
        b[s] = a[inds[s]]
    end
    b
end
function PaddedMatrices.∂getindex(
    a::AbstractFixedSizeVector{M,T},
    mdv::MissingDataVector{T,S,M}
) where {T,S,M}
    inds = mdv.indices
    b = FixedSizeVector{S,T}(undef)
    @inbounds for s ∈ 1:S
        b[s] = a[inds[s]]
    end
    b, MissingDataVectorAdjoint(mdv)
end
function Base.getindex(
    sp::StackPointer,
    a::AbstractFixedSizeVector{M,T},
    mdv::MissingDataVector{T,S,M}
) where {T,S,M}
    inds = mdv.indices
    sp, b = PaddedMatrices.PtrVector{S,T}(sp)
    @inbounds for s ∈ 1:S
        b[s] = a[inds[s]]
    end
    sp, b
end
function PaddedMatrices.∂getindex(
    sp::StackPointer,
    a::AbstractFixedSizeVector{M,T},
    mdv::MissingDataVector{T,S,M}
) where {T,S,M}
    inds = mdv.indices
    sp, b = PaddedMatrices.PtrVector{S,T}(sp)
    @inbounds for s ∈ 1:S
        b[s] = a[inds[s]]
    end
    sp, (b, MissingDataVectorAdjoint(mdv))
end
@generated function Base.getindex(
    sp::StackPointer,
    a::NTuple{K,V},
    mdv::MissingDataVector{T,S,M}
) where {T,K,S,M,V<:PaddedMatrices.AbstractMutableFixedSizeVector{M,T}}
    quote
        inds = mdv.indices
        Base.Cartesian.@nexprs $K k -> (a_k = @inbounds a[k]; (sp, b_k) = PaddedMatrices.PtrVector{$S,$T}(sp))
        for s ∈ 1:$S
            i = inds[s]
            Base.Cartesian.@nexprs $K k-> (b_k[s] = a_k[i])
        end
        sp, (Base.Cartesian.@ntuple $K b)
    end
end

@generated function PaddedMatrices.∂getindex(
    sp::StackPointer,
    a::NTuple{K,V},
    mdv::MissingDataVector{T,S,M}
) where {T,K,S,M,V<:PaddedMatrices.AbstractFixedSizeVector{M,T}}
    quote
        inds = mdv.indices
        Base.Cartesian.@nexprs $K k -> (a_k = @inbounds a[k]; (sp, b_k) = PaddedMatrices.PtrVector{$S,$T}(sp))
        @inbounds for s ∈ 1:$S
            i = inds[s]
            Base.Cartesian.@nexprs $K k-> (b_k[s] = a_k[i])
        end
        sp, ((Base.Cartesian.@ntuple $K b), MissingDataVectorAdjoint(mdv))
    end
end

function Base.:*(A::AbstractCovarianceMatrix, mdv::MissingDataVectorAdjoint)
    ∂Σ = mdv.∂Σ
    mask!(∂Σ, A, mdv.indices)
    ∂Σ
end
@generated function Base.:*(
    a::NTuple{K,V},
    mdv::MissingDataVectorAdjoint{T,S,M}
) where {K,T,S,M,V<:PaddedMatrices.AbstractFixedSizeVector{S,T}}
    quote
        Base.Cartesian.@nexprs $K k -> (b_k = zeros(FixedSizeVector{$M,$T}); a_k = a[k])
        inds = mdv.indices
        @inbounds for s ∈ 1:$S
            i = inds[s]
            Base.Cartesian.@nexprs $K k -> b_k[i] = a_k[s]
        end
        Base.Cartesian.@ntuple $K b
    end
end
@generated function Base.:*(
    sp::StackPointer,
    a::NTuple{K,V},
    mdv::MissingDataVectorAdjoint{T,S,M}
) where {K,T,S,M,V<:PaddedMatrices.AbstractFixedSizeVector{S,T}}
    Wm1 = VectorizationBase.pick_vector_width(M,T) - 1
    total_length = (K*M + Wm1) & + ~Wm1
    quote
        # we zero the used data in a single loop.
        zero_out = PtrVector{$total_length,$T,$total_length}(pointer(sp, $T))
        @inbounds @simd ivdep for i ∈ 1:$total_length
            zero_out[i] = zero($T)
        end
        Base.Cartesian.@nexprs $K k -> begin
            b_k = PtrVector{$M,$T,$M}(pointer(sp,$T) + $(sizeof(T)*M) * (k-1))
            a_k = a[k]
        end
        inds = mdv.indices
        @inbounds for s ∈ 1:$S
            i = inds[s]
            Base.Cartesian.@nexprs $K k -> b_k[i] = a_k[s]
        end
        sp + $(VectorizationBase.align(sizeof(T)*M*K)), (Base.Cartesian.@ntuple $K b)
    end
end

@generated function DynamicCovarianceMatrix(
                rhos::PaddedMatrices.AbstractFixedSizeVector{K,T}, L::StructuredMatrices.AbstractLowerTriangularMatrix{K},
                times::ConstantFixedSizeVector{nT}, workspace
            ) where {K,T,nT}
    Wm1 = VectorizationBase.pick_vector_width(nT, T) - 1
    nTl = PaddedMatrices.calc_padding(nT, T) #(nT + Wm1) & ~Wm1
    quote
        ARs = workspace.ARs
        @fastmath @inbounds for k ∈ 1:$K
            ρ = rhos[k]
            logρ = $(T == Float64 ? :(ccall(:log,Float64,(Float64,),ρ))  :  :(Base.log(ρ)))
            for tc ∈ 1:nT
                for tr ∈ 1:tc-1
                    ARs[tr,tc,k] = ARs[tc,tr,k]
                end
                ARs[tc,tc,k] = one(T)
                for tr ∈ tc+1:nT
                    deltatimes = times[tr] - times[tc]
                    logrhot = deltatimes*logρ
                    ARs[tr,tc,k] = $(T == Float64 ? :(ccall(:exp,Float64,(Float64,),logrhot)) : :(Base.exp(logrhot)))
                end
            end
        end
        Sigfull = workspace.Sigfull
        @inbounds begin
        Base.Cartesian.@nexprs $K kc -> begin
            for kr ∈ kc-1:K-1
                Base.Cartesian.@nexprs kc j -> l_j = L[kr+1,j] * L[kc,j]
                sigrow, sigcol = $nT*kr, $nT*(kc-1)
                for tc ∈ 1:$nT
                    @simd ivdep for tr ∈ 1:$nTl
                        ari = l_1 * ARs[tr, tc, 1]
                        Base.Cartesian.@nexprs kc-1 j -> ari = muladd(l_{j+1}, ARs[tr, tc, j+1], ari)
                        Sigfull[ tr + sigrow, tc + sigcol] = ari
                    end
                end
            end
        end
        end
        DynamicCovarianceMatrix(Sigfull)
    end
end

struct Covariance_LAR_AR_Adjoint{K, nT, T, nTP, L, MFPD<:AbstractMutableFixedSizeArray{Tuple{nT,nT,K},T,3,nTP,L}, LKJ_T <: StructuredMatrices.AbstractLowerTriangularMatrix}
    ∂ARs::MFPD#FixedSizeArray{Tuple{nT,nT,K},T,3,nTP,L}
    LKJ::LKJ_T
end
struct Covariance_LAR_L_Adjoint{K, nT, T, nTP, L, MFPD<:AbstractMutableFixedSizeArray{Tuple{nT,nT,K},T,3,nTP,L}, LKJ_T <: StructuredMatrices.AbstractLowerTriangularMatrix}
    ARs::MFPD#FixedSizeArray{Tuple{nT,nT,K},T,3,nTP,L}
    LKJ::LKJ_T
end

@generated function ∂DynamicCovarianceMatrix(
    rhos::PaddedMatrices.AbstractFixedSizeVector{K,T},
    L::StructuredMatrices.AbstractLowerTriangularMatrix{K,T},
    times::AbstractFixedSizeVector{nT},
    workspace, ::Val{(true,true)}
) where {K,T,nT}
    # We will assume rho > 0
    Wm1 = VectorizationBase.pick_vector_width(nT, T) - 1
    nTl = PaddedMatrices.calc_padding(nT, T) #(nT + Wm1) & ~Wm1
    quote
        # nT = length(times) + 1
        ARs = workspace.ARs
        ∂ARs = workspace.∂ARs
        @inbounds @fastmath for k ∈ 1:$K
            ρ = rhos[k]
            # We want to use Base.log, and not fastmath log.
            logρ = $(T == Float64 ? :(ccall(:log,Float64,(Float64,),ρ))  :  :(Base.log(ρ)))
            for tc ∈ 1:nT
                for tr ∈ 1:tc-1
                    ARs[tr,tc,k] = ARs[tc,tr,k]
                    ∂ARs[tr,tc,k] = ∂ARs[tc,tr,k]
                end
                ARs[tc,tc,k] = one(T)
                ∂ARs[tc,tc,k] = zero(T)
                for tr ∈ tc+1:nT
                    deltatimes = times[tr] - times[tc]
                    # We want to specifically use Base.exp, and not fastmath exp
                    logrhot = (deltatimes - one($T))*logρ
                    rhot = $(T == Float64 ? :(ccall(:exp,Float64,(Float64,),logrhot)) : :(Base.exp(logrhot)))
                    ARs[tr,tc,k] = ρ*rhot
                    ∂ARs[tr,tc,k] = deltatimes*rhot
                end
            end
        end
        Sigfull = workspace.Sigfull
        @inbounds begin
        Base.Cartesian.@nexprs $K kc -> begin
            for kr ∈ kc-1:K-1
                Base.Cartesian.@nexprs kc j -> l_j = L[kr+1,j] * L[kc,j]
                sigrow, sigcol = $nT*kr, $nT*(kc-1)
                for tc ∈ 1:$nT
                    @simd ivdep for tr ∈ 1:$nTl
                        ari = l_1 * ARs[tr, tc, 1]
                        Base.Cartesian.@nexprs kc-1 j -> ari = muladd(l_{j+1}, ARs[tr, tc, j+1], ari)
                        Sigfull[ tr + sigrow, tc + sigcol] = ari
                    end
                end
            end
        end
        end
        DynamicCovarianceMatrix(Sigfull), Covariance_LAR_AR_Adjoint(∂ARs,L), Covariance_LAR_L_Adjoint(ARs,L)
    end
end

@generated function CovarianceMatrix(
    sp::StackPointer,
    rhos::PaddedMatrices.AbstractFixedSizeVector{K,T},
    L::StructuredMatrices.AbstractLowerTriangularMatrix{K,T},
    times::AbstractFixedSizeVector{nT},
    ::Val{SIMD} = Val{true}()
) where {K,T,nT,SIMD}
    # We will assume rho > 0
    W, Wshift = VectorizationBase.pick_vector_width_shift(nT, T)
    Wm1 = W - 1
    V = Vec{W,T}
    nTl = PaddedMatrices.calc_padding(nT, T) #(nT + Wm1) & ~Wm1
    KT = K*nT
    T_size = sizeof(T)
    quote
        # only Sigfull escapes, so we allocate it first
        # and return the stack pointer pointing to its end.
        sp, Sigfull = PtrFixedSizeCovarianceMatrix{$KT,$T}(sp)
        ARs = PaddedMatrices.PtrArray{Tuple{$nT,$nT,$K},$T,3}(pointer(sp, $T))

        ptr_time = pointer(times)
        ptr_ARs = pointer(ARs)
        @inbounds for k ∈ 1:$K
            ρ = rhos[k]
            # We want to use Base.log, and not fastmath log.
            logρ = $(T == Float64 ? :(ccall(:log,Float64,(Float64,),ρ))  :  :(Base.log(ρ)))
            $(SIMD ? quote
            vρ = SIMDPirates.vbroadcast($V, ρ)
            vlogρ = SIMDPirates.vbroadcast($V, logρ)
            for tcouter ∈ 0:$((nTl>>>Wshift)-1)
                for tcinner ∈ 1:$W
                    Wtcouter = $W*tcouter
                    tc = tcinner + Wtcouter
                    tc > $nT && break
                    for tr ∈ 1:Wtcouter
                        ARs[tr,tc,k] = ARs[tc,tr,k]
                    end
                    vtime_tc = SIMDPirates.vbroadcast($V, times[tc])
                    for i ∈ tcouter:$((nTl>>>Wshift)-1)
                        vtimes_tr = SIMDPirates.vload($V, ptr_time + i * $(T_size*W))
                        vδt = SIMDPirates.vabs(SIMDPirates.vsub(vtimes_tr, vtime_tc))
                        vδtm1 = SIMDPirates.vsub(vδt, SIMDPirates.vbroadcast($V,one($T)))
                        vρt = SLEEFPirates.exp(SIMDPirates.vmul(vδtm1, vlogρ))
                        AR_offset = (k-1)*$(T_size*nTl*nT) + (tc-1)*$(T_size*nTl) + i*$(W*T_size)

                        SIMDPirates.vstore!(ptr_ARs + AR_offset, SIMDPirates.vmul(vρ,vρt))
                    end
                end
              end
              end : quote
            for tc ∈ 1:$nT
                for tr ∈ 1:tc-1
                    ARs[tr,tc,k] = ARs[tc,tr,k]
                end
                ARs[tc,tc,k] = one(T)
                for tr ∈ tc+1:$nT
                    deltatimes = times[tr] - times[tc]
                    #                    rhot = rho^(deltatimes - one(T))
                    # We want to specifically use Base.exp, and not fastmath exp
                    logrhot = (deltatimes - one($T))*logρ
                    rhot = $(T == Float64 ? :(ccall(:exp,Float64,(Float64,),logrhot)) : :(Base.exp(logrhot)))
                    ARs[tr,tc,k] = ρ*rhot
                end
              end
              end )
        end
        @inbounds begin
            Base.Cartesian.@nexprs $K kc -> begin
                for kr ∈ kc-1:K-1
                    Base.Cartesian.@nexprs kc j -> l_j = L[kr+1,j] * L[kc,j]
                    sigrow, sigcol = $nT*kr, $nT*(kc-1)
                    for tc ∈ 0:$(nT-1)
                        @vectorize for tr ∈ 1:$nT
                            ari = l_1 * ARs[tr + tc * $nTl]
                            Base.Cartesian.@nexprs kc-1 j -> ari = muladd(l_{j+1}, ARs[tr + tc*$nTl + j*$(nTl*nT)], ari)
                            Sigfull[ tr + sigrow + (tc + sigcol)*$KT] = ari
                        end
                    end
                end
            end
        end
        sp, Sigfull
    end
end
@generated function ∂CovarianceMatrix(
    sp::StackPointer,
    rhos::PaddedMatrices.AbstractFixedSizeVector{K,T},
    L::StructuredMatrices.AbstractLowerTriangularMatrix{K,T},
    times::AbstractFixedSizeVector{nT},
#    missing::MissingDataMatrix,
    ::Val{(true,true)},::Val{SIMD} = Val{true}()
) where {K,T,nT,SIMD}
    # We will assume rho > 0
    W, Wshift = VectorizationBase.pick_vector_width_shift(nT, T)
    Wm1 = W - 1
    V = Vec{W,T}
    nTl = PaddedMatrices.calc_padding(nT, T)
    KT = K*nT
    T_size = sizeof(T)
    sym = true
    ARquote = if SIMD && sym
        quote
            vρ = SIMDPirates.vbroadcast($V, ρ)
            vlogρ = SIMDPirates.vbroadcast($V, logρ)
            for tcouter ∈ 0:$((nTl>>>Wshift)-1)
                for tcinner ∈ 1:$W
                    Wtcouter = $W*tcouter
                    tc = tcinner + Wtcouter
                    tc > $nT && break
                    for tr ∈ 1:Wtcouter
                        ARs[tr,tc,k] = ARs[tc,tr,k]
                        ∂ARs[tr,tc,k] = ∂ARs[tc,tr,k]
                    end
                    vtime_tc = SIMDPirates.vbroadcast($V, times[tc])
                    for i ∈ tcouter:$((nTl>>>Wshift)-1)
                        vtimes_tr = SIMDPirates.vload($V, ptr_time + i * $(T_size*W))
                        vδt = SIMDPirates.vabs(SIMDPirates.vsub(vtimes_tr, vtime_tc))
                        vδtm1 = SIMDPirates.vsub(vδt, SIMDPirates.vbroadcast($V,one($T)))
                        vρt = SLEEFPirates.exp(SIMDPirates.vmul(vδtm1, vlogρ))
                        AR_offset = (k-1)*$(T_size*nTl*nT) + (tc-1)*$(T_size*nTl) + i*$(W*T_size)

                        SIMDPirates.vstore!(ptr_ARs + AR_offset, SIMDPirates.vmul(vρ,vρt))
                        SIMDPirates.vstore!(ptr_∂ARs + AR_offset, SIMDPirates.vmul(vδt,vρt))
                    end
                end
              end
        end
    elseif SIMD
        quote
            vρ = SIMDPirates.vbroadcast($V, ρ)
            vlogρ = SIMDPirates.vbroadcast($V, logρ)
            @inbounds for c ∈ 0:$(nT-1)
                offset = (k-1)*$(T_size*nTl*nT) + c*$(T_size*nTl)
                vtime_tc = SIMDPirates.vbroadcast($V, times[c+1])
                for r ∈ 0:$((nTl >>> Wshift)-1)
                    vtimes_tr = SIMDPirates.vload($V, ptr_time + r * $(T_size*W))
                    vδt = SIMDPirates.vabs(SIMDPirates.vsub(vtimes_tr, vtime_tc))
                    vδtm1 = SIMDPirates.vsub(vδt, SIMDPirates.vbroadcast($V,one($T)))
                    vρt = SLEEFPirates.exp(SIMDPirates.vmul(vδtm1, vlogρ))
                    AR_offset = offset + r*$(W*T_size)

                    SIMDPirates.vstore!(ptr_ARs + AR_offset, SIMDPirates.vmul(vρ,vρt))
                    SIMDPirates.vstore!(ptr_∂ARs + AR_offset, SIMDPirates.vmul(vδt,vρt))
                end
            end
        end
    else
        quote
            for tc ∈ 1:$nT
                for tr ∈ 1:tc-1
                    ARs[tr,tc,k] = ARs[tc,tr,k]
                    ∂ARs[tr,tc,k] = ∂ARs[tc,tr,k]
                end
                ARs[tc,tc,k] = one(T)
                ∂ARs[tc,tc,k] = zero(T)
                for tr ∈ tc+1:$nT
                    deltatimes = times[tr] - times[tc]
                    #                    rhot = rho^(deltatimes - one(T))
                    # We want to specifically use Base.exp, and not fastmath exp
                    logrhot = (deltatimes - one($T))*logρ
                    rhot = $(T == Float64 ? :(ccall(:exp,Float64,(Float64,),logrhot)) : :(Base.exp(logrhot)))
                    ARs[tr,tc,k] = ρ*rhot
                    ∂ARs[tr,tc,k] = deltatimes*rhot
                end
            end
        end
    end
    quote
        sp,  ARs = PaddedMatrices.PtrArray{Tuple{$nT,$nT,$K},$T,3}(sp)
        sp, ∂ARs = PaddedMatrices.PtrArray{Tuple{$nT,$nT,$K},$T,3}(sp)
        ptr_time = pointer(times)
        ptr_ARs = pointer(ARs)
        ptr_∂ARs = pointer(∂ARs)
        @inbounds for k ∈ 1:$K
            ρ = rhos[k]
            # We want to use Base.log, and not fastmath log.
            logρ = $(T == Float64 ? :(ccall(:log,Float64,(Float64,),ρ))  :  :(Base.log(ρ)))
            $ARquote
        end
        sp, Sigfull = PtrFixedSizeCovarianceMatrix{$KT,$T}(sp)
        @inbounds begin
            Base.Cartesian.@nexprs $K kc -> begin
                for kr ∈ kc-1:K-1
                    Base.Cartesian.@nexprs kc j -> l_j = L[kr+1,j] * L[kc,j]
                    sigrow, sigcol = $nT*kr, $nT*(kc-1)
                    for tc ∈ 0:$(nT-1)
                        @vectorize for tr ∈ 1:$nT
                            ari = l_1 * ARs[tr + tc * $nTl]
                            Base.Cartesian.@nexprs kc-1 j -> ari = muladd(l_{j+1}, ARs[tr + tc*$nTl + j*$(nTl*nT)], ari)
                            Sigfull[ tr + sigrow + (tc + sigcol)*$KT] = ari
                        end
                    end
                end
            end
        end
        sp, (Sigfull, Covariance_LAR_AR_Adjoint(∂ARs,L), Covariance_LAR_L_Adjoint(ARs,L))
    end
end

@generated function Base.:*(
    C::AbstractMatrix{T},
    adj::Covariance_LAR_AR_Adjoint{K, nT, T, nTP}
) where {K, nT, T, nTP}
    # C is a DynamicCovarianceMatrix
    Wm1 = VectorizationBase.pick_vector_width(T)-1
    KL = PaddedMatrices.calc_padding(K, T) #(K+Wm1) & ~Wm1
    outtup = Expr(:tuple, [Symbol(:κ_,k) for k ∈ 1:K]..., [zero(T) for k ∈ K+1:KL]...)
    quote
        $(Expr(:meta,:inline))
        # Calculate all in 1 pass?
        L = adj.LKJ
        ∂ARs = adj.∂ARs
        Base.Cartesian.@nexprs $K j -> κ_j = zero($T)
        @inbounds begin
            Base.Cartesian.@nexprs $K kc -> begin
                # Diagonal block
                Base.Cartesian.@nexprs kc j -> l_j = L[kc,j]^2
                ccol = crow = $nT*(kc-1)
                for tc ∈ 1:$nT
                    tco = tc + ccol
                    for tr ∈ tc+1:$nT
                        cij = C[tr + crow, tco]
                        Base.Cartesian.@nexprs kc j -> begin
                            κ_j = Base.FastMath.add_fast(Base.FastMath.mul_fast(cij,l_j,∂ARs[tr,tc,j]), κ_j)
                        end
                    end
                end

                for kr ∈ kc:$(K-1)
                    Base.Cartesian.@nexprs kc j -> l_j = L[kc,j]*L[kr+1,j]
                    crow = $nT*kr; ccol = $nT*(kc-1)
                    for tc ∈ 1:$nT
                        tco = tc + ccol
                        for tr ∈ 1:$nT
                            cij = C[tr + crow, tco]
                            Base.Cartesian.@nexprs kc j -> begin
                                κ_j = Base.FastMath.add_fast(Base.FastMath.mul_fast(cij,l_j,∂ARs[tr,tc,j]), κ_j)
                            end
                        end
                    end
                end
            end
        end
        ConstantFixedSizeVector{$K,$T,$KL,$KL}( $outtup )'
    end
end
@generated function Base.:*(
    sp::StackPointer,
    C::AbstractMatrix{T},
    adj::Covariance_LAR_AR_Adjoint{K, nT, T, nTP}
) where {K, nT, T, nTP}
    # C is a DynamicCovarianceMatrix
    Wm1 = VectorizationBase.pick_vector_width(T)-1
    KL = (K+Wm1) & ~Wm1
    quote
        $(Expr(:meta,:inline))
        # Calculate all in 1 pass?
        L = adj.LKJ
        ∂ARs = adj.∂ARs
        Base.Cartesian.@nexprs $K j -> κ_j = zero($T)
        @inbounds begin
            Base.Cartesian.@nexprs $K kc -> begin
                # Diagonal block
                Base.Cartesian.@nexprs kc j -> l_j = L[kc,j]^2
                ccol = crow = $nT*(kc-1)
                for tc ∈ 1:$nT
                    tco = tc + ccol
                    for tr ∈ tc+1:$nT
                        cij = C[tr + crow, tco]
                        Base.Cartesian.@nexprs kc j -> begin
                            κ_j = Base.FastMath.add_fast(Base.FastMath.mul_fast(cij,l_j,∂ARs[tr,tc,j]), κ_j)
                        end
                    end
                end

                for kr ∈ kc:$(K-1)
                    Base.Cartesian.@nexprs kc j -> l_j = L[kc,j]*L[kr+1,j]
                    crow = $nT*kr; ccol = $nT*(kc-1)
                    for tc ∈ 1:$nT
                        tco = tc + ccol
                        for tr ∈ 1:$nT
                            cij = C[tr + crow, tco]
                            Base.Cartesian.@nexprs kc j -> begin
                                κ_j = Base.FastMath.add_fast(Base.FastMath.mul_fast(cij,l_j,∂ARs[tr,tc,j]), κ_j)
                            end
                        end
                    end
                end
            end
            b = PaddedMatrices.PtrVector{$K,$T,$KL}(pointer(sp,$T))
            Base.Cartesian.@nexprs $K k -> b[k] = κ_k
        end
        sp + $(VectorizationBase.align(KL*sizeof(T))), b'
    end
end

function Base.:*(
    C::AbstractMatrix{T},
    adj::Covariance_LAR_L_Adjoint{K, nT, T, nTP}
) where {K, nT, T, nTP}
    ∂LKJ = StructuredMatrices.MutableLowerTriangularMatrix{K,T}(undef)
    LKJ = adj.LKJ
    ARs = adj.ARs
    # Outer loops are for ∂lkj
    for kc ∈ 1:K
        @inbounds for kr ∈ kc:K
            ∂lkj = zero(T)
            # Move left
            for lr ∈ kc:kr-1
                lkj = LKJ[lr,kc]
                for tc ∈ 1:nT
                    # C: kr = row, lr = col
                    @fastmath for tr ∈ 1:nT
                        ∂lkj += lkj * ARs[tr,tc,kc] * C[tr + (kr-1)*nT,tc + (lr-1)*(nT)]
                    end
                end
            end
            # squared diagonal derivative
            lkj = 2LKJ[kr,kc]
            for tc ∈ 1:nT
                # C: kr = row, lr = col
                @fastmath ∂lkj += lkj * ARs[tc,tc,kc] * C[tc + (kr-1)*nT,tc + (kr-1)*(nT)]# * T(0.5)  #* (tc == tr ? T(0.5) : T(1.0))
                @fastmath for tr ∈ tc+1:nT
                    ∂lkj += lkj * ARs[tr,tc,kc] * C[tr + (kr-1)*nT,tc + (kr-1)*(nT)] #* (tc == tr ? T(0.5) : T(1.0))
                end
            end
            # move down
            for lr ∈ kr+1:K
                lkj = LKJ[lr,kc]
                for tc ∈ 1:nT
                    # C: lr = row, kr = col
                    @fastmath for tr ∈ 1:nT
                        ∂lkj += lkj * ARs[tr,tc,kc] * C[tr + (lr-1)*nT,tc + (kr-1)*(nT)]
                    end
                end
            end
            ∂LKJ[kr,kc] = ∂lkj
        end
    end
    ∂LKJ
end

function Base.:*(
    sp::StackPointer,
    C::AbstractMatrix{T},
    adj::Covariance_LAR_L_Adjoint{K, nT, T, nTP}
) where {K, nT, T, nTP}

    (sp,∂LKJ) = StructuredMatrices.PtrLowerTriangularMatrix{K,T}(sp)

    LKJ = adj.LKJ
    ARs = adj.ARs

    # Outer loops are for ∂lkj
    for kc ∈ 1:K
        # coffset = (kc-1)*stride(C,2)
        @inbounds for kr ∈ kc:K
            ∂lkj = zero(T)
            # Move left
            for lr ∈ kc:kr-1
                lkj = LKJ[lr,kc]
                for tc ∈ 1:nT
                    # C: kr = row, lr = col
                    @fastmath for tr ∈ 1:nT
                        ∂lkj += lkj * ARs[tr,tc,kc] * C[tr + (kr-1)*nT,tc + (lr-1)*(nT)]
                    end
                end
            end
            # squared diagonal derivative
            lkj = 2LKJ[kr,kc]
            for tc ∈ 1:nT
                # C: kr = row, lr = col
                @fastmath ∂lkj += lkj * ARs[tc,tc,kc] * C[tc + (kr-1)*nT,tc + (kr-1)*(nT)] #* (tc == tr ? T(0.5) : T(1.0))
                @fastmath for tr ∈ tc+1:nT
                    ∂lkj += lkj * ARs[tr,tc,kc] * C[tr + (kr-1)*nT,tc + (kr-1)*(nT)] #* (tc == tr ? T(0.5) : T(1.0))
                end
            end
            # move down
            for lr ∈ kr+1:K
                lkj = LKJ[lr,kc]
                for tc ∈ 1:nT
                    # C: lr = row, kr = col
                    @fastmath for tr ∈ 1:nT
                        ∂lkj += lkj * ARs[tr,tc,kc] * C[tr + (lr-1)*nT,tc + (kr-1)*(nT)]
                    end
                end
            end
            ∂LKJ[kr,kc] = ∂lkj
        end
    end
#    @show ∂LKJ
    sp, ∂LKJ
end


text:
#!/usr/bin/env julia
# This is a Julia Processing server to do non-linearity correction of HPF data

using NPZ
@everywhere using Dierckx
@everywhere using JLD
using Memento
include("/storage/home/jpn23/scratch/JuliaServerForMath.jl/src/JuliaProcessingServer.jl")
using JuliaProcessingServer

logger = Memento.config!("info",fmt="[{date} | {level}]: {msg}")
push!(logger, DefaultHandler("/storage/home/j/jpn23/scratch/logfile_JuliaServerForMath.log",DefaultFormatter("[{date} | {level}]: {msg}")))

info(logger,"Loading the NLcorrTCKdic from jdl file")

# #Load and merge the Nonlinearity dictionaires in parallel
# NLcorrTCKdic = @parallel merge for i in 0:128:1920
#     jldfilename = "/storage/home/j/jpn23/scratch/LongRamp/BSPLINE_NLC_20171117_$i-$(i+128).jld"
#     load(jldfilename, "NLcorrTCKdic")
# end

#Load the NLCdictionaries as a list without merging
ListOfNLcorrTCKdic = pmap((fname)->load(fname,"NLcorrTCKdic"), ["/storage/home/j/jpn23/scratch/LongRamp/BSPLINE_NLC_20171117_$i-$(i+128).jld" for i in 0:128:1920])

# NLcorrTCKdic = load("/storage/home/j/jpn23/scratch/LongRamp/BSPLINE_NLC_20171117_0-128.jld", "NLcorrTCKdic")
info(logger,"Finished Loading the NLcorrTCKdic from jdl file")
function ApplyNonLinearityCorrection(datacube)
    datacube = SharedArray(datacube)
    @parallel for NLcorrTCKdic in ListOfNLcorrTCKdic
        for ij in keys(NLcorrTCKdic)
            tck = NLcorrTCKdic[ij]
            # Since Julia has 1 indexing, we should add +1 to python coorindates to match julia's coordinates
    	    datacube[:,ij[1]+1,ij[2]+1] = Dierckx._evaluate(tck[1],tck[2],tck[3],datacube[:,ij[1]+1,ij[2]+1],0)
        end
    end
    return datacube
end

#################################
ListenSOCKET = 8006
ReaderFunction = NPZ.npzreadarray   # Function to read .npy file received via socket
ProcessFunction = ApplyNonLinearityCorrection # Function to apply on to the recived data
PackerFunction = NPZ.npzwritearray # Function to write .npy file to be send via socket
##################################

open("/storage/home/jpn23/scratch/NonLinCorrServerStarted.txt", "w") do f
    write(f, "NonLinCorr Julia Server Started\n")
end
StartProcessingServer(ListenSOCKET,ReaderFunction,ProcessFunction,PackerFunction;logger=logger)

ext: jl


text:
module JoinTest1

using HPAT

#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(3)

@acc hpat function test1(file_name, file_name2)
    t1 = DataSource(DataTable{:userid=Int64, :val2=Float64}, HDF5, file_name)
    t2 = DataSource(DataTable{:userid=Int64, :val3=Float64}, HDF5, file_name2)
    t3 = join(t1,t2, :userid==:userid, :userid)
    return t3[:userid], t3[:val3]
end

using HDF5
using MPI
using Base.Test

function main()

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        h5write("test1_1.hdf5", "/userid", [1,2,3,1,2])
        h5write("test1_1.hdf5", "/val2", [1.1,2.1,3.1,3.2,1.9])
        h5write("test1_2.hdf5", "/userid", [1,3])
        h5write("test1_2.hdf5", "/val3", [7.1,8.3])
    end
    a,b = test1("test1_1.hdf5","test1_2.hdf5")
    println(a)
    println(b)

    @test a==[1,1,3]
    @test_approx_eq b [7.1,7.1,8.3]

    if MPI.Comm_rank(MPI.COMM_WORLD)==0
        rm("test1_1.hdf5")
        rm("test1_2.hdf5")
    end

end

end

JoinTest1.main()


text:
module google
    using ..ProtoBuf
    import ..ProtoBuf: meta
    module protobuf
        using ..ProtoBuf
        import ..ProtoBuf: meta

        include("empty_pb.jl")
        include("any_pb.jl")
        include("timestamp_pb.jl")
        include("duration_pb.jl")
        include("wrappers_pb.jl")
        include("source_context_pb.jl")
        include("type_pb.jl")
        include("field_mask_pb.jl")
        include("api_pb.jl")
        include("struct_pb.jl")
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -16.424759279035918 <= q <= 20.23452089871623)

                     add_NL_constraint(m, :((0.8420310557117703 + 0.34198141401583015*(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.1887519896617813*(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.2882490310085508*(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))))/(1 + exp(-(0.8420310557117703 + 0.34198141401583015*(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.1887519896617813*(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.2882490310085508*(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))))))) + (-0.35188345522271236 + -0.11415509244934263*(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.8088728193024375*(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + 0.8271426964720727*(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))))/(1 + exp(-(-0.35188345522271236 + -0.11415509244934263*(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.8088728193024375*(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + 0.8271426964720727*(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))))))) + (0.4798469380896484 + 0.8134725364327324*(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.12585836968848696*(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + 0.613011911779838*(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))))/(1 + exp(-(0.4798469380896484 + 0.8134725364327324*(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.5389217244927917 + 0.5798025520034438*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + 0.4428632765655749*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + 0.9033224385695577*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + -0.12585836968848696*(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.23921485861490321 + 0.5030319708107833*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.2081481501163327*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.2939216839670169*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))))))) + 0.613011911779838*(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))/(1 + exp(-(-0.6545847926105979 + 0.9774628615021363*(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4]))/(1 + exp(-(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])))) + -0.21180887804098658*(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4]))/(1 + exp(-(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])))) + -0.7013494059895882*(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))/(1 + exp(-(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function expected_improvement(modl::KrigingModel, fmin::Real, PointToExamine::Array{<:Real,1})
    # This is equation 15 of Jones, Scholau, Welch
    output = modl(PointToExamine)
    fmin_m_y = fmin - output.value
    fmin_m_y_on_s = fmin_m_y / output.standard_error
    return fmin_m_y * Distributions.cdf(Normal(), fmin_m_y_on_s) + output.standard_error *  Distributions.pdf(Normal(), fmin_m_y_on_s)
end

function get_next_query_point_through_expected_improvement(modl::KrigingModel, lower =  repeat([0.0], size(modl.X)[2]), upper = repeat([1.0], size(modl.X)[2]) )
    fmin = minimum(modl.f)
    opt = optimize(x -> -expected_improvement(modl, fmin, x), lower, upper, (lower .+ upper) ./ 2, SAMIN(), Optim.Options(iterations = 10000))
    return opt.minimizer
end

function get_predicted_minimum(modl::KrigingModel, lower =  repeat([0.0], size(modl.X)[2]), upper = repeat([1.0], size(modl.X)[2]) )
    opt = optimize(x -> -modl(x).value, lower, upper, (lower .+ upper) ./ 2, SAMIN(), Optim.Options(iterations = 10000))
    return opt.minimizer
end


text:
"Recursive structure for representing goals."
mutable struct GoalTree
    term::Term # Term to be proven
    parent::Union{GoalTree,Nothing} # Parent goal
    children::Vector{Term} # List of subgoals
    active::Int # Index of active subgoal
    env::Subst # Dictionary of variable mappings for the current goal
    vmap::Subst # Variables inherited from parent
end

GoalTree(g::GoalTree) =
    GoalTree(g.term, g.parent, copy(g.children), g.active, copy(g.env), g.vmap)

"Built-in arithmetic operations."
const math_ops = Set([:+, :-, :*, :/, :mod])
"Built-in comparison operations."
const comp_ops = Set([:(==), :<=, :>=, :<, :>, :(!=)])
"Built-in operators."
const ops = union(math_ops, comp_ops)
"Built-in functions."
const default_funcs = Dict(op => eval(op) for op in ops)
"Built-in logical connectives."
const logicals = Set([true, false, :and, :or, :not, :!,
                      :exists, :forall, :imply, :(=>)])
"Built-in predicates with special handling during SLD resolution."
const builtins = union(comp_ops, logicals,
    Set([:is, :call, :unifies, :≐, :cut, :fail, :findall, :countall]))

"""
    eval_term(term, env[, funcs])

Given an environment, evaluate all variables in a Julog term to constants.
Returns a term that is as fully evaluated as possible.

# Arguments
- `term::Term`: A term to evaluate.
- `env::Dict{Var,Term}`: An environment mapping variables to terms.
- `funcs::Dict=Dict()`: Additional custom functions (e.g. custom math).
"""
eval_term(term::Term, env::Subst, funcs::Dict=Dict()) = error("Not implemented")
eval_term(term::Var, env::Subst, funcs::Dict=Dict()) =
    term in keys(env) ? eval_term(env[term], env, funcs) : term
function eval_term(term::Const, env::Subst, funcs::Dict=Dict())
    val = get(funcs, term.name, term)
    if val === term
        return val
    elseif isa(val, Function)
        return Const(val())
    else
        return Const(val)
    end
end
function eval_term(term::Compound, env::Subst, funcs::Dict=Dict())
    args = Term[eval_term(a, env, funcs) for a in term.args]
    func = get(funcs, term.name) do
        get(default_funcs, term.name, nothing)
    end
    if func !== nothing && all(isa(a, Const) for a in args)
        if isa(func, Function)
            # Evaluate function if all arguments are fully evulated
            return Const(func((a.name for a in args)...))
        elseif isa(func, Dict)
            # Lookup value if custom function is a lookup table
            key = Tuple(a.name for a in args)
            if key in keys(func) return Const(func[key]) end
            # Leave un-evaluated if lookup table has missing entries
            return Compound(term.name, args)
        else
            error("$(term.name) is neither a custom function or lookup table.")
        end
    else
        return Compound(term.name, args)
    end
end

"""
    unify(src, dst[, occurs_check])

Unifies src with dst and returns a dictionary of any substitutions needed.

# Arguments
- `src::Term`: A term to unify.
- `dst::Term`: A term to unify (src/dst order does not matter).
- `occurs_check::Bool=true`: Whether to perform the occurs check.
- `funcs::Dict=Dict()`: Custom functions to evaluate.
"""
function unify(src::Term, dst::Term,
               occurs_check::Bool=true, funcs::Dict=Dict())
    src_stack, dst_stack = Term[src], Term[dst]
    src_defer, dst_defer = Term[], Term[]
    subst = Subst()
    success = true
    while length(src_stack) > 0
        src, dst = pop!(src_stack), pop!(dst_stack)
        # Unify $src with $dst
        if isa(src, Const) && isa(dst, Const)
            success = _unify!(src::Const, dst, funcs)
            if !success break end
        elseif isa(src, Var)
            success = _unify!(src::Var, dst, src_stack, dst_stack,
                              subst, occurs_check)
            if !success break end
        elseif isa(dst, Var)
            success = _unify!(dst::Var, src, dst_stack, src_stack,
                              subst, occurs_check)
            if !success break end
        elseif isa(src, Compound) && isa(dst, Compound)
            success = _unify!(src::Compound, dst, src_stack, dst_stack)
            if !success break end
        elseif isa(src, Const) && isempty(dst.args)
            success = _unify!(src::Const, dst, funcs)
            if !success break end
        elseif isa(dst, Const) && isempty(src.args)
            success = _unify!(dst::Const, src, funcs)
            if !success break end
        else # Reaches here if one term is compound and the other is constant
            push!(src_defer, src)
            push!(dst_defer, dst)
        end
    end
    # Try evaluating deferred Const vs Compound comparisons
    for (src, dst) in zip(src_defer, dst_defer)
        src, dst = eval_term(src, subst, funcs), eval_term(dst, subst, funcs)
        if src != dst # "No: cannot unify $src and $dst after evaluation"
            success = false; break
        end
    end
    return success ? subst : nothing
end

function _unify!(src::Const, dst::Term, funcs::Dict=Dict())
     if src.name != dst.name
         src_val = get(funcs, src.name, src.name)
         dst_val = get(funcs, dst.name, dst.name)
         return src_val == dst_val
     end
     return true
end

function _unify!(src::Var, dst::Term, src_stack, dst_stack,
                 subst, occurs_check::Bool=true)
    if isa(dst, Var) && src.name == dst.name # "Yes: same variable"
        return true
    elseif occurs_check && occurs_in(src, dst) # "No: src occurs in dst"
        return false
    end
    # Replace src with dst in stack
    for i in 1:length(src_stack)
        s, d = src_stack[i], dst_stack[i]
        ss = substitute(s, src, dst)
        if s !== ss src_stack[i] = ss end
        dd = substitute(d, src, dst)
        if d !== dd dst_stack[i] = dd end
    end
    # Replace src with dst in substitution values
    map!(v -> substitute(v, src, dst), values(subst))
    # Add substitution of src to dst
    subst[src] = dst
    return true
end

function _unify!(src::Compound, dst::Term, src_stack, dst_stack)
    if src.name != dst.name # "No: diff functors"
        return false
    elseif length(src.args) != length(dst.args) # "No: diff arity"
        return false
    end
    # "Yes: pushing args onto stack"
    append!(src_stack, src.args)
    append!(dst_stack, dst.args)
    return true
end

"Handle built-in predicates"
function handle_builtins!(queue, clauses, goal, term; options...)
    funcs = get(options, :funcs, Dict())
    occurs_check = get(options, :occurs_check, false)
    vcount = get(options, :vcount, Ref(UInt(0)))
    if term.name == true
        return true
    elseif term.name == false
        return false
    elseif term.name == :call
        # Handle meta-call predicate
        pred, args = term.args[1], term.args[2:end]
        if isa(pred, Var) pred = get(goal.env, pred, nothing) end
        if isnothing(pred) error("$term not sufficiently instantiated.") end
        term = Compound(pred.name, [pred.args; args]) # Rewrite term
        goal.children[goal.active] = term # Replace term
        goal.active -= 1
        return true
    elseif term.name == :is
        # Handle is/2 predicate
        qn, ans = term.args[1], eval_term(term.args[2], goal.env, funcs)
        # Failure if RHS is insufficiently instantiated
        if !isa(ans, Const) return false end
        # LHS can either be a variable or evaluate to a constant
        if isa(qn, Var) && !(qn in keys(goal.env))
            # If LHS is a free variable, bind to RHS
            for (k, v) in goal.env
                goal.env[k] = substitute(v, qn, ans)
            end
            goal.env[qn] = ans
            return true
        else
            # If LHS evaluates to a constant, check if it is equal to RHS
            qn = eval_term(qn, goal.env, funcs)
            return isa(qn, Const) ? qn == ans : false
        end
        return false
    elseif term.name in [:unifies, :≐]
        # Check if LHS and RHS unify
        term = substitute(term, goal.env)
        lhs, rhs = term.args[1], term.args[2]
        unifier = unify(lhs, rhs, occurs_check, funcs)
        if isnothing(unifier) return false end
        compose!(goal.env, unifier) # Update variable bindings if satisfied
        return true
    elseif term.name == :and
        # Remove self and add all arguments as children to the goal, succeed
        splice!(goal.children, goal.active, term.args)
        goal.active -= 1
        return true
    elseif term.name == :or
        # Create new goals for each term in disjunct, add to queue
        for arg in term.args
            g = GoalTree(goal)
            g.children[g.active] = arg # Replace disjunct with term
            push!(queue, g)
        end
        # No longer work on this goal
        return false
    elseif term.name in [:not, :!]
        # Try to resolve negated predicate, return true upon failure
        neg_goal = term.args[1]
        sat, _ = resolve(Term[neg_goal], clauses; options...,
                         vcount=vcount, env=copy(goal.env), mode=:any)
        return !sat # Success if no proof is found
    elseif term.name == :exists
        # exists(Cond, Body) holds if Body holds for at least 1 binding of Cond
        cond, body = term.args
        sat, subst = resolve(Term[cond, body], clauses; options...,
                             vcount=vcount, env=copy(goal.env), mode=:any)
        # Update variable bindings if satisfied
        goal.env = sat ? compose!(goal.env, subst[1]) : goal.env
        return sat
    elseif term.name == :forall
        # forall(Cond, Body) holds if Body holds for all bindings of Cond
        cond, body = term.args
        term = @julog(not(and(:cond, not(:body)))) # Rewrite term
        goal.children[goal.active] = term # Replace term
        goal.active -= 1
        return true
    elseif term.name in [:imply, :(=>)]
        # imply(Cond, Body) holds if or(not(Cond), Body) holds
        cond, body = term.args
        sat, _ = resolve(Term[cond], clauses; options...,
                         vcount=vcount, env=copy(goal.env), mode=:any)
        # Return true if Cond does not hold
        if !sat return true end
        # Otherwise replace original term with [Cond, Body], return true
        splice!(goal.children, goal.active, [cond, body])
        goal.active -= 1
        return true
    elseif term.name == :cut
        # Remove all other goals and succeed
        empty!(queue)
        return true
    elseif term.name == :fail
        # Fail and skip goal
        return false
    elseif term.name == :findall
        # Find list of all cond matches, substituted into template
        template, cond, list = term.args
        _, subst = resolve(Term[cond], clauses; options...,
                           vcount=vcount, env=copy(goal.env), mode=:all)
        matches = to_term_list([substitute(template, s) for s in subst])
        unifier = unify(list, matches, occurs_check, funcs)
        if isnothing(unifier) return false end
        compose!(goal.env, unifier) # Update variable bindings if satisfied
        return true
    elseif term.name == :countall
        # Count number of ways to prove the condition
        cond, count = term.args
        _, subst = resolve(Term[cond], clauses; options...,
                           vcount=vcount, env=copy(goal.env), mode=:all)
        unifier = unify(count, Const(length(subst)), occurs_check, funcs)
        if isnothing(unifier) return false end
        compose!(goal.env, unifier) # Update variable bindings if satisfied
        return true
    elseif term.name in comp_ops || term.name in keys(funcs)
        defer_eval = get(options, :defer_eval, false)
        result = eval_term(term, goal.env, funcs)
        if isa(result, Const)
            # Return true if term evaluates to true
            return result.name == true
        elseif defer_eval && goal.active < length(goal.children)
            # Defer evaluation
            goal.children[goal.active] = goal.children[goal.active + 1]
            goal.children[goal.active + 1] = term
            goal.active -= 1
            return true
        end
    end
    return false
end

"""
    resolve(goals, clauses; <keyword arguments>)
    bwd_chain(goals, clauses; <keyword arguments>)

SLD-resolution of goals with additional Prolog-like control flow.

# Arguments
- `goals::Vector{<:Term}`: A list of Julog terms to be prove or query.
- `clauses::Vector{Clause}`: A list of Julog clauses.
- `env::Subst=Subst()`: An initial environment mapping variables to terms.
- `funcs::Dict=Dict()`: Custom functions for evaluating terms.
  A function `f` should be stored as `funcs[:f] = f`
- `mode::Symbol=:all`: How results should be returned.
  `:all` returns all possible substitiutions. `:any` returns the first
  satisfying substitution found. `:interactive` prompts for continuation
  after each satisfying substitution is found.
- `search::Symbol=:bfs`: search either breadth (`:bfs`) or depth-first (`:dfs`)
- `occurs_check::Bool=false`: Flag for occurs check during unification
- `defer_eval::Bool=false`: Flag to defer evaluation of (custom) operators.
"""
function resolve(goals::Vector{<:Term}, clauses::Vector{<:AbstractClause}; options...)
    return resolve(goals, index_clauses(clauses); options...)
end

function resolve(goals::Vector{<:Term}, clauses::ClauseTable; options...)
    # Unpack options
    env = get(options, :env, Subst())::Subst
    funcs = get(options, :funcs, Dict())::Dict
    mode = get(options, :mode, :all)::Symbol
    search = get(options, :search, :bfs)::Symbol
    occurs_check = get(options, :occurs_check, false)::Bool
    vcount = get(options, :vcount, Ref(UInt(0)))::Ref{UInt}
    # Construct top level goal and put it on the queue
    queue = [GoalTree(Const(false), nothing, Vector{Term}(goals), 1, env, Subst())]
    subst = Subst[]
    # Iterate across queue of goals
    while length(queue) > 0
        goal = (search == :dfs) ? pop!(queue) : popfirst!(queue)
        @debug string("Goal: ", Clause(goal.term, goal.children), " ",
                      "Env: ", goal.env)
        if goal.active > length(goal.children)
            # All subgoals are done
            if isnothing(goal.parent)
                # If goal has no parent, we are successful
                @debug string("Success: ", goal.env)
                if !(goal.env in subst)
                    push!(subst, goal.env)
                end
                if (mode == :all) continue
                elseif (mode == :any || length(queue) == 0) break
                elseif mode == :interactive
                    # Interactively prompt for continuation
                    println(goal.env)
                    print("Continue? [y/n]: ")
                    s = readline()
                    if s in ["y", "Y", ""] continue
                    else break
                    end
                end
            end
            @debug string("Done, returning to parent.")
            # Copy construct parent to create fresh copy of environment
            parent = GoalTree(goal.parent)
            # Remap variables from child to parent
            vmap = Subst()
            for (pvar, cvar) in goal.vmap
                if !(cvar in keys(goal.env)) continue end
                if pvar == goal.env[cvar] continue end
                vmap[pvar] = goal.env[cvar]
            end
            parent.env = compose!(parent.env, vmap)
            # Advance parent to next subgoal and put it back on the queue
            parent.active += 1
            push!(queue, parent)
            continue
        end
        # Examine current subgoal
        term = goal.children[goal.active]
        @debug string("Subgoal: ", term)
        # Handle built-in special terms
        if term.name in builtins || term.name in keys(funcs)
            success = handle_builtins!(queue, clauses, goal, term;
                                       vcount=vcount, options...)
            # If successful, go to next subgoal, otherwise skip to next goal
            if success
                @debug string("Done, returning to parent.")
                goal.active += 1
                push!(queue, goal)
            end
            continue
        end
        # Substitute and freshen variables in term
        vmap = Subst()
        term = freshen!(substitute(term, goal.env), vmap, vcount)
        # Iterate across clause set with matching heads
        matched_clauses = retrieve_clauses(clauses, term, funcs)
        matched = false
        for c in matched_clauses
            # Freshen variables in clause
            c = freshen!(c, Subst(), vcount)
            # If term unifies with head of a clause, add it as a subgoal
            unifier = unify(term, c.head, occurs_check, funcs)
            if isnothing(unifier) continue end
            child = GoalTree(c.head, goal, copy(c.body), 1, unifier, vmap)
            push!(queue, child)
            matched = true
        end
        if !matched @debug string("Failed, no matching clauses.") end
    end
    # Goals were satisfied if we found valid substitutions
    return (length(subst) > 0), subst
end

resolve(goal::Term, clauses::Vector{<:AbstractClause}; options...) =
    resolve(Term[goal], clauses; options...)

resolve(goals::Vector{Clause}, clauses::Vector{<:AbstractClause}; options...) =
    resolve([convert(Term, g) for g in goals], clauses; options...)

const bwd_chain = resolve

"Return all facts derivable in `n` steps from the initial set of clauses."
function derivations(clauses::Vector{<:AbstractClause}, n::Real=1; options...)
    rules = filter(c -> length(c.body) > 0, clauses)
    facts = filter(c -> length(c.body) == 0, clauses)
    return derivations(rules, facts, n; options...)
end

function derivations(rules::Vector{T}, facts::Vector{T}, n::Real=1;
                     options...) where {T <: AbstractClause}
    return derivations(rules, index_clauses(facts), n; options...)
end

function derivations(rules::Vector{T}, facts::ClauseTable{T}, n::Real=1;
                     as_indexed::Bool=false, options...) where {T <: AbstractClause}
    step = 0
    n_facts = num_clauses(facts)
    while step < n
        # Iteratively add facts derivable from each rule
        derived = derive_step(rules, facts; options...)
        insert_clauses!(facts, derived)
        # Terminate early if we reach a fixed point
        n_facts_new = num_clauses(facts)
        if (n_facts < n_facts_new) n_facts = n_facts_new else break end
        step += 1
    end
    return as_indexed ? facts : deindex_clauses(facts)
end

"Iteratively add facts derivable from each rule."
function derive_step(rules::Vector{T}, facts::ClauseTable{T};
                     options...) where {T <: AbstractClause}
    # Iteratively add facts derivable from each rule
    derived = Term[]
    for r in rules
        # Find all valid substitutions of each rule's body
        _, subst = resolve(r.body, facts; options...)
        append!(derived, Term[substitute(r.head, s) for s in subst])
    end
    return [Clause(d, []) for d in derived]
end

"""
    derive(goals, clauses; <keyword arguments>)
    fwd_chain(goals, clauses; <keyword arguments>)

Derive goals via forward-chaining from the initial set of clauses.

# Arguments
- `goals::Vector{<:Term}`: A list of Julog terms to be prove or query.
- `clauses::Vector{Clause}`: A list of Julog clauses.
- `max-steps::Real=100`: Maximum steps before terminating with failure.

See `resolve` for other supported arguments.
"""
function derive(goals::Vector{<:Term}, clauses::Vector{T};
                max_steps::Real=100, options...) where {T <: AbstractClause}
    grounded = all(is_ground.(goals))
    rules = filter(c -> length(c.body) > 0, clauses)
    facts = index_clauses(filter(c -> length(c.body) == 0, clauses))
    n_facts = num_clauses(facts)
    step = 0
    while step < max_steps
        # Return all goals are grounded and satisfied
        if grounded
            sat, subst = resolve(goals, facts; options...)
            if sat return sat, subst end
        end
        # Iteratively add facts derivable from each rule
        derived = derive_step(rules, facts; options...)
        insert_clauses!(facts, derived)
        # Terminate early if we reach a fixed point
        n_facts_new = num_clauses(facts)
        if (n_facts < n_facts_new) n_facts = n_facts_new else break end
        step += 1
    end
    return resolve(goals, facts; options...)
end

derive(goal::Term, clauses::Vector{<:AbstractClause}; options...) =
    derive(Term[goal], clauses; options...)

derive(goals::Vector{Clause}, clauses::Vector{<:AbstractClause}; options...) =
    derive([convert(Term, g) for g in goals], clauses; options...)

const fwd_chain = derive


text:
using Primes

function smallest_multiple(n)
    prime_powers = Dict()

    for i in 2:n
        prime_factors = factor(i)
        for (prime, power) in prime_factors
            prime_powers[prime] = max(power, get(prime_powers, prime, 0))
        end
    end

    println(prod([prime^power for (prime, power) in prime_powers]))
end

smallest_multiple(20)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
@testset "Search" begin
    t = SegmentMap{Int,Int}()
    @test !haskey(t, (1, 2))

    n = 10000
    global maxend = 1000000
    segments = [randsegment(1, maxend) for i in 1:n]

    @testset "true positives" begin
        for (i, segment) in enumerate(segments)
            t[segment] = i
        end

        shuffle!(segments)
        @test all(Bool[haskey(t, segment) for segment in segments])
        @test all(Bool[findfirst(t, segment) !== nothing for segment in segments])

        results = Bool[]
        for segment in segments
            x = findfirst(t, segment)
            push!(results, x.first == segment[1] && x.last == segment[2])
        end
        @test all(results)

        @test all(Bool[(get!(t, segment, -1) != SegmentValue{Int,Int}(segment[1], segment[2], -1))
                       for segment in segments])
        @test all(Bool[get(t, segment, -1) != -1 for segment in segments])
    end

    @testset "true negatives" begin
        @test all(Bool[!haskey(t, segment)
                       for segment in [randsegment(maxend + 1, 2 * maxend)
                                        for i in 1:n]])
        @test all(Bool[(findfirst(t, segment) === nothing)
                       for segment in [randsegment(maxend + 1, 2 * maxend)
                                        for i in 1:n]])

        @test all(Bool[(findfirst(t, segment, (a, b) -> false) === nothing) for segment in segments])

        @test all(Bool[get!(t, segment, -1) == SegmentValue{Int,Int}(segment[1], segment[2], -1)
                       for segment in [randsegment(maxend + 1, 2 * maxend)
                                        for i in 1:n]])
        @test all(Bool[get(t, segment, -1) == -1
                       for segment in [randsegment(maxend + 1, 2 * maxend)
                                        for i in 1:n]])
    end
end


text:
function Skewness()
    return Skewness(())
end

function Skewness(arg0::Skewness)
    return Skewness((Skewness,), arg0)
end

function Skewness(arg0::ThirdMoment)
    return Skewness((ThirdMoment,), arg0)
end

function clear(obj::Skewness)
    return jcall(obj, "clear", void, ())
end

function copy(obj::Skewness)
    return jcall(obj, "copy", Skewness, ())
end

function evaluate(obj::Skewness, arg0::Vector{jdouble}, arg1::jint, arg2::jint)
    return jcall(obj, "evaluate", jdouble, (Vector{jdouble}, jint, jint), arg0, arg1, arg2)
end

function get_n(obj::Skewness)
    return jcall(obj, "getN", jlong, ())
end

function get_result(obj::Skewness)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::Skewness, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/MarketDamagesBurke.jl")

    marketdamages = addmarketdamagesburke(m)

    set_param!(m, :MarketDamagesBurke, :rtl_realizedtemperature, readpagedata(m, "test/validationdata/$valdir/rtl_realizedtemperature.csv"))
    set_param!(m, :MarketDamagesBurke, :rcons_per_cap_SLRRemainConsumption, readpagedata(m,"test/validationdata/$valdir/rcons_per_cap_SLRRemainConsumption.csv"))
    set_param!(m, :MarketDamagesBurke, :rgdp_per_cap_SLRRemainGDP, readpagedata(m,"test/validationdata/$valdir/rgdp_per_cap_SLRRemainGDP.csv"))
    set_param!(m, :MarketDamagesBurke, :isatg_impactfxnsaturation, 28.333333333333336)

    p = load_parameters(m)
    p["y_year_0"] = 2015.
    p["y_year"] = Mimi.dim_keys(m.md, :time)
    set_leftover_params!(m, p)

    run(m)

    rcons_per_cap = m[:MarketDamagesBurke, :rcons_per_cap_MarketRemainConsumption]
    rcons_per_cap_compare = readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_MarketRemainConsumption.csv")
    @test rcons_per_cap ≈ rcons_per_cap_compare rtol=1e-1

    rgdp_per_cap = m[:MarketDamagesBurke, :rgdp_per_cap_MarketRemainGDP]
    rgdp_per_cap_compare = readpagedata(m, "test/validationdata/$valdir/rgdp_per_cap_MarketRemainGDP.csv")
    @test rgdp_per_cap ≈ rgdp_per_cap_compare rtol=1e-2
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE