text:
abstract type StochasticDataSeries{T,N} <: AbstractArray{T,N} end

#Type for holding a stochastic data series:
# nd   - dimension of the space of data
# nt+1 - number of time steps in the DataSeries including the initial condition
# ns   - number of sample paths of the stochastic process
# ni   - number of the initial conditions
# d    - the array holding the data
# T    - type of the held data
# N    - the number of dimensions of the array d
struct SStochasticDataSeries{T,N} <: StochasticDataSeries{T,N}
    nd::Int
    nt::Int
    ns::Int
    ni::Int
    d::Array{T,N}

    function SStochasticDataSeries{T,N}(nd, nt, ns, ni) where {T,N}
        @assert T <: Number
        @assert nd > 0
        @assert nt ≥ 0
        @assert ns ≥ 0
        @assert ni ≥ 0

        @assert N ∈ (1,2,3,4)

        if N == 1
            d = zeros(T, nt+1)
        elseif N == 2
            d = zeros(T, nd, nt+1)
        elseif N == 3
            ni==1 ? d = zeros(T, nd, nt+1, ns) : d = zeros(T, nd, nt+1, ni)
        elseif N==4
            d = zeros(T, nd, nt+1, ns, ni)
        end

        new(nd, nt, ns, ni, d)
    end

    function SStochasticDataSeries{T,N}(nd, nt, ns, ni, d) where {T,N}
        new(nd, nt, ns, ni, d)
    end
end

function SStochasticDataSeries(T, nd::Int, nt::Int, ns::Int, ni::Int)

    if nd==ns==ni==1
        N=1
    elseif ns==ni==1
        N=2
    elseif ns!=1!=ni
        N=4
    else
        N=3
    end

    return SStochasticDataSeries{T,N}(nd, nt, ns, ni)
end

function SStochasticDataSeries(d::Array{T,1}) where {T}
    nd = 1
    nt = size(d,1)-1
    ns = 1
    ni = 1
    return SStochasticDataSeries{T,1}(nd, nt, ns, ni, d)
end

function SStochasticDataSeries(d::Array{T,2}) where {T}
    nd = size(d,1)
    nt = size(d,2)-1
    ns = 1
    ni = 1
    return SStochasticDataSeries{T,2}(nd, nt, ns, ni, d)
end

# A 3-dimensional matrix d can represent a data series with ns=1 or ni=1.
# The second argument specifies whether there are multiple initial conditions
# (IC=true, so ni>1) or a single one (default IC=false, so ni=1)
function SStochasticDataSeries(d::Array{T,3}; IC=false) where {T}
    nd = size(d,1)
    nt = size(d,2)-1

    if IC==false
        ns = size(d,3)
        ni = 1
    else
        ni = size(d,3)
        ns = 1
    end

    return SStochasticDataSeries{T,3}(nd, nt, ns, ni, d)
end

function SStochasticDataSeries(d::Array{T,4}) where {T}
    nd = size(d,1)
    nt = size(d,2)-1
    ns = size(d,3)
    ni = size(d,4)
    return SStochasticDataSeries{T,4}(nd, nt, ns, ni, d)
end

#TODO: PStochasticDataSeries

function Base.show(io::IO, ds::StochasticDataSeries{T,N}) where {T,N}
    print(io, "StochasticDataSeries with data type ", T, " and ", N, " dimensions:\n")
    print(io, ds.d)
end

function similar(ds::StochasticDataSeries{T,N}) where {T,N}
    typeof(ds){T,N}(ds.nd, ds.nt, ds.ns, ds.ni)
end

Base.parent(ds::StochasticDataSeries) = ds.d
Base.eltype(ds::StochasticDataSeries{T,N}) where {T,N} = T
Base.ndims(ds::StochasticDataSeries{T,N}) where {T,N} = N

errmsg_size(ds::StochasticDataSeries) = error("size not supported for data series with axes $(axes(ds))")
Base.size(ds::StochasticDataSeries) = errmsg_size(ds)
Base.size(ds::StochasticDataSeries, d) = errmsg_size(ds)

Base.eachindex(::IndexCartesian, ds::StochasticDataSeries) = CartesianIndices(axes(ds))
Base.eachindex(::IndexLinear, ds::StochasticDataSeries) = axes(ds, 1)

Base.lastindex(ds::StochasticDataSeries{T,1}) where {T} = (ds.nt)
Base.lastindex(ds::StochasticDataSeries{T,2}) where {T} = (ds.nd, ds.nt)
Base.lastindex(ds::StochasticDataSeries{T,3}) where {T} = ds.ni==1 ? (ds.nd, ds.nt, ds.ns) : (ds.nd, ds.nt, ds.ni)
Base.lastindex(ds::StochasticDataSeries{T,4}) where {T} = (ds.nd, ds.nt, ds.ns, ds.ni)

@inline Base.axes(ds::StochasticDataSeries{T,1}) where {T} = (0:ds.nt)
@inline Base.axes(ds::StochasticDataSeries{T,2}) where {T} = (1:ds.nd, 0:ds.nt)
@inline Base.axes(ds::StochasticDataSeries{T,3}) where {T} = ds.ni==1 ? (1:ds.nd, 0:ds.nt, 1:ds.ns) : (1:ds.nd, 0:ds.nt, 1:ds.ni)
@inline Base.axes(ds::StochasticDataSeries{T,4}) where {T} = (1:ds.nd, 0:ds.nt, 1:ds.ns, 1:ds.ni)

Base.strides(ds::StochasticDataSeries) = (strides(ds.d))

function get_data!(ds::StochasticDataSeries{T,1}, n) where {T}
    j = n+1
    @assert j ≤ size(ds.d, 1)
    @inbounds return ds.d[j]
end

function get_data!(ds::StochasticDataSeries{T,2}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n) where {T}
    j = n+1
    @assert j ≤ size(ds.d, 2)
    @inbounds for i in eachindex(x)
        x[i] = ds.d[i,j]
    end
end

function get_data!(ds::StochasticDataSeries{T,3}, x::Union{Array{T,2}, Array{TwicePrecision{T},2}}, n) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @assert j ≤ size(ds.d, 2)
    @inbounds for k in axes(x, 2)
        for i in axes(x, 1)
            x[i,k] = ds.d[i,j,k]
        end
    end
end

function get_data!(ds::StochasticDataSeries{T,3}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n, k) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert j ≤ size(ds.d, 2)
    @assert k ≤ size(ds.d, 3)
    @inbounds for i in eachindex(x)
        x[i] = ds.d[i,j,k]
    end
end

function get_data!(ds::StochasticDataSeries{T,4}, x::Union{Array{T,3}, Array{TwicePrecision{T},3}}, n) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @assert size(x,3) == size(ds.d, 4)
    @assert j ≤ size(ds.d, 2)
    @inbounds for m in axes(x,3)
        for k in axes(x, 2)
            for i in axes(x, 1)
                x[i,k,m] = ds.d[i,j,k,m]
            end
        end
    end
end

# Copies the data for the time step n and initial condition m from ds.d to x
function get_data!(ds::StochasticDataSeries{T,4}, x::Union{Array{T,2}, Array{TwicePrecision{T},2}}, n, m) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @assert j ≤ size(ds.d, 2)
    @assert m ≤ size(ds.d, 4)
    @inbounds for k in axes(x, 2)
        for i in axes(x, 1)
            x[i,k] = ds.d[i,j,k,m]
        end
    end
end

function get_data!(ds::StochasticDataSeries{T,4}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n, k, m) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert j ≤ size(ds.d, 2)
    @assert k ≤ size(ds.d, 3)
    @assert m ≤ size(ds.d, 4)
    @inbounds for i in eachindex(x)
        x[i] = ds.d[i,j,k,m]
    end
end

# Assigns a number to the nth time slot in a one-dimensional data series
function set_data!(ds::StochasticDataSeries{T,1}, x::Union{T, TwicePrecision{T}}, n) where {T}
    j = n+1
    @assert j ≤ size(ds.d, 1)
    @inbounds ds.d[j] = x
end

# Assigns a number to the nth time slot in a one-dimensional data series,
# when the number is passed as a single-entry array
function set_data!(ds::StochasticDataSeries{T,1}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n) where {T}
    @assert length(x)==1
    set_data!(ds,x[1],n)
end

# Assigns a vector to the nth time slot in a two-dimensional data series
function set_data!(ds::StochasticDataSeries{T,2}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n) where {T}
    j = n+1
    @assert length(x) == size(ds.d, 1)
    @assert j ≤ size(ds.d, 2)
    @inbounds for i in axes(ds.d, 1)
        ds.d[i,j] = x[i]
    end
end

# Assigns a two dimensional array (nd x ns or nd x ni) to the nth time slot in a 3-dimensional data series
function set_data!(ds::StochasticDataSeries{T,3}, x::Union{Array{T,2}, Array{TwicePrecision{T},2}}, n) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @assert j ≤ size(ds.d, 2)
    @inbounds for k in axes(ds.d, 3)
        for i in axes(ds.d, 1)
            ds.d[i,j,k] = x[i,k]
        end
    end
end

# Assigns a one dimensional array (nd elements) to the nth time slot and the kth sample path in a 3-dimensional data series
function set_data!(ds::StochasticDataSeries{T,3}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n, k) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert j ≤ size(ds.d, 2)
    @assert k ≤ size(ds.d, 3)
    @inbounds for i in axes(ds.d, 1)
        ds.d[i,j,k] = x[i]
    end
end

# Assigns a one dimensional array (nd elements) to the nth time slot and to all sample paths in a 3-dimensional data series
function set_data!(ds::StochasticDataSeries{T,3}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert j ≤ size(ds.d, 2)
    @inbounds for k in axes(ds.d, 3)
        @inbounds for i in axes(ds.d, 1)
            ds.d[i,j,k] = x[i]
        end
    end
end

# Assign an nd x ns x ni array to the nth time slot in a 4-dimensional data series
function set_data!(ds::StochasticDataSeries{T,4}, x::Union{Array{T,3}, Array{TwicePrecision{T},3}}, n) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @assert size(x,3) == size(ds.d, 4)
    @assert j ≤ size(ds.d, 2)
    @inbounds for m in axes(x,3)
        for k in axes(x, 2)
            for i in axes(x, 1)
                ds.d[i,j,k,m] = x[i,k,m]
            end
        end
    end
end

# Assings the data stored in x to the data series ds as the time step n and for
# the initial condition m
function set_data!(ds::StochasticDataSeries{T,4}, x::Union{Array{T,2}, Array{TwicePrecision{T},2}}, n, m) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @assert j ≤ size(ds.d, 2)
    @assert m ≤ size(ds.d, 4)
    @inbounds for k in axes(x, 2)
        for i in axes(x, 1)
            ds.d[i,j,k,m] = x[i,k]
        end
    end
end

# Assings an nd x ni matrix to the data series ds as the time step n and for
# all sample paths
function set_data!(ds::StochasticDataSeries{T,4}, x::Union{Array{T,2}, Array{TwicePrecision{T},2}}, n) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 4)
    @assert j ≤ size(ds.d, 2)
    @inbounds for m in axes(ds.d, 3)
        @inbounds for k in axes(x, 2)
            for i in axes(x, 1)
                ds.d[i,j,m,k] = x[i,k]
            end
        end
    end
end

function set_data!(ds::StochasticDataSeries{T,4}, x::Union{Array{T,1}, Array{TwicePrecision{T},1}}, n, k, m) where {T}
    j = n+1
    @assert size(x,1) == size(ds.d, 1)
    @assert j ≤ size(ds.d, 2)
    @assert k ≤ size(ds.d, 3)
    @assert m ≤ size(ds.d, 4)
    @inbounds for i in eachindex(x)
        ds.d[i,j,k,m] = x[i]
    end
end

function reset!(ds::StochasticDataSeries{T,1}) where {T}
    # ds[0]=ds[end-1] was wrong for some reason
    # ds[0]=ds[end] would be ok, but to be safe below using ds.d
    @inbounds ds.d[1] = ds.d[end]
end

function reset!(ds::StochasticDataSeries{T,2}) where {T}
    # when 'end' is passed to getindex, because of the shift it becomes end+1
    # in ds.d, therefore end-1 has to be passed to getindex instead
	# -> should be fixed now (but check!)
    @inbounds for i in axes(ds,1)
        ds[i,0] = ds[i,end]
    end
end

function reset!(ds::StochasticDataSeries{T,3}) where {T}
    @inbounds for k in axes(ds,3)
        for i in axes(ds,1)
            ds[i,0,k] = ds[i,end,k]
        end
    end
end

function reset!(ds::StochasticDataSeries{T,4}) where {T}
    @inbounds for m in axes(ds,4)
        for k in axes(ds,3)
            for i in axes(ds,1)
                ds[i,0,k,m] = ds[i,end,k,m]
            end
        end
    end
end

@inline function Base.getindex(ds::StochasticDataSeries{T,1}, j::Int) where {T}
    @boundscheck checkbounds(ds.d, j+1)
    @inbounds r = getindex(ds.d, j+1)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,2}, i::Int, j::Int) where {T}
    @boundscheck checkbounds(ds.d, i, j+1)
    @inbounds r = getindex(ds.d, i, j+1)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,2}, j::Int) where {T}
    @boundscheck checkbounds(ds.d, :, j+1)
    @inbounds r = getindex(ds.d, :, j+1)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,3}, i::Int, j::Int, k::Int) where {T}
    @boundscheck checkbounds(ds.d, i, j+1, k)
    @inbounds r = getindex(ds.d, i, j+1, k)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,3}, j::Int, k::Int) where {T}
    @boundscheck checkbounds(ds.d, :, j+1, k)
    @inbounds r = getindex(ds.d, :, j+1, k)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,3}, k::Int) where {T}
    @boundscheck checkbounds(ds.d, :, :, k)
    @inbounds r = getindex(ds.d, :, :, k)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,4}, i::Int, j::Int, k::Int, m::Int) where {T}
    @boundscheck checkbounds(ds.d, i, j+1, k, m)
    @inbounds r = getindex(ds.d, i, j+1, k, m)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,4}, j::Int, k::Int, m::Int) where {T}
    @boundscheck checkbounds(ds.d, :, j+1, k, m)
    @inbounds r = getindex(ds.d, :, j+1, k, m)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,4}, j::Int, m::Int) where {T}
    @boundscheck checkbounds(ds.d, :, j+1, :, m)
    @inbounds r = getindex(ds.d, :, j+1, :, m)
    return r
end

@inline function Base.getindex(ds::StochasticDataSeries{T,4}, m::Int) where {T}
    @boundscheck checkbounds(ds.d, :, :, :, m)
    @inbounds r = getindex(ds.d, :, :, :, m)
    return r
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,1}, x, j::Int) where {T}
    @assert length(x) == 1
    @boundscheck checkbounds(ds.d, j+1)
    @inbounds setindex!(ds.d, x, j+1)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,2}, x, i::Int, j::Int) where {T}
    @assert length(x) == 1
    @boundscheck checkbounds(ds.d, i, j+1)
    @inbounds setindex!(ds.d, x, i, j+1)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,2}, x, j::Int) where {T}
    @assert ndims(x) == 1
    @assert length(x) == size(ds.d, 1)
    @boundscheck checkbounds(ds.d, :, j+1)
    @inbounds setindex!(ds.d, x, :, j+1)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,3}, x, i::Int, j::Int, k::Int) where {T}
    @assert length(x) == 1
    @boundscheck checkbounds(ds.d, i, j+1, k)
    @inbounds setindex!(ds.d, x, i, j+1, k)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,3}, x, j::Int, k::Int) where {T}
    @assert ndims(x) == 1
    @assert length(x) == size(ds.d, 1)
    @boundscheck checkbounds(ds.d, :, j+1, k)
    @inbounds setindex!(ds.d, x, :, j+1, k)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,3}, x, j::Int) where {T}
    @assert ndims(x) == 2
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @boundscheck checkbounds(ds.d, :, j+1, :)
    @inbounds setindex!(ds.d, x, :, j+1, :)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,4}, x, i::Int, j::Int, k::Int, m::Int) where {T}
    @assert length(x) == 1
    @boundscheck checkbounds(ds.d, i, j+1, k, m)
    @inbounds setindex!(ds.d, x, i, j+1, k, m)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,4}, x, j::Int, k::Int, m::Int) where {T}
    @assert ndims(x) == 1
    @assert length(x) == size(ds.d, 1)
    @boundscheck checkbounds(ds.d, :, j+1, k, m)
    @inbounds setindex!(ds.d, x, :, j+1, k, m)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,4}, x, j::Int, m::Int) where {T}
    @assert ndims(x) == 2
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @boundscheck checkbounds(ds.d, :, j+1, :, m)
    @inbounds setindex!(ds.d, x, :, j+1, :, m)
end

@inline function Base.setindex!(ds::StochasticDataSeries{T,4}, x, j::Int) where {T}
    @assert ndims(x) == 3
    @assert size(x,1) == size(ds.d, 1)
    @assert size(x,2) == size(ds.d, 3)
    @assert size(x,3) == size(ds.d, 4)
    @boundscheck checkbounds(ds.d, :, j+1, :, :)
    @inbounds setindex!(ds.d, x, :, j+1, :, :)
end

# TODO Implement convert() to/from array.


text:
import .MathOptInterface
const MOI = MathOptInterface

struct MOIOptimizationProblem{T,F<:OptimizationFunction,uType,P} <: MOI.AbstractNLPEvaluator
    f::F
    u0::uType
	p::P
    J::Matrix{T}
    H::Matrix{T}
    cons_H::Vector{Matrix{T}}
end

MOI.eval_objective(moiproblem::MOIOptimizationProblem, x) = moiproblem.f(x, moiproblem.p)

MOI.eval_constraint(moiproblem::MOIOptimizationProblem, g, x) = g .= moiproblem.f.cons(x)

MOI.eval_objective_gradient(moiproblem::MOIOptimizationProblem, G, x) = moiproblem.f.grad(G, x)

function MOI.eval_hessian_lagrangian(moiproblem::MOIOptimizationProblem{T}, h, x, σ, μ) where {T}
    n = length(moiproblem.u0)
    a = zeros(n, n)
    moiproblem.f.hess(a, x)
    if iszero(σ)
        fill!(h, zero(T))
    else
        moiproblem.f.hess(moiproblem.H, x)
        k = 0
        for col in 1:n
            for row in 1:col
                k += 1
                h[k] = σ * moiproblem.H[row, col]
            end
        end
    end
    if !isempty(μ) && !all(iszero, μ)
        moiproblem.f.cons_h(moiproblem.cons_H, x)
        for (μi, Hi) in zip(μ, moiproblem.cons_H)
            k = 0
            for col in 1:n
                for row in 1:col
                    k += 1
                    h[k] += μi * Hi[row, col]
                end
            end
        end
    end
    return
end

function MOI.eval_constraint_jacobian(moiproblem::MOIOptimizationProblem, j, x)
    isempty(j) && return
    moiproblem.f.cons_j === nothing && error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends")
    n = length(moiproblem.u0)
    moiproblem.f.cons_j(moiproblem.J, x)
    for i in eachindex(j)
        j[i] = moiproblem.J[i]
    end
end

function MOI.jacobian_structure(moiproblem::MOIOptimizationProblem)
    return Tuple{Int,Int}[(con, var) for con in 1:size(moiproblem.J,1) for var in 1:size(moiproblem.J,2)]
end

function MOI.hessian_lagrangian_structure(moiproblem::MOIOptimizationProblem)
    return Tuple{Int,Int}[(row, col) for col in 1:length(moiproblem.u0) for row in 1:col]
end

function MOI.initialize(moiproblem::MOIOptimizationProblem, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in MOI.features_available(moiproblem))
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

MOI.features_available(moiproblem::MOIOptimizationProblem) = [:Grad, :Hess, :Jac]

function make_moi_problem(prob::OptimizationProblem)
    num_cons = prob.ucons === nothing ? 0 : length(prob.ucons)
    f = instantiate_function(prob.f,prob.u0,prob.f.adtype,prob.p,num_cons)
    T = eltype(prob.u0)
    n = length(prob.u0)
    moiproblem = MOIOptimizationProblem(f,prob.u0,prob.p,zeros(T,num_cons,n),zeros(T,n,n),Matrix{T}[zeros(T,n,n) for i in 1:num_cons])
    return moiproblem
end

function __map_optimizer_args(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MOI.OptimizerWithAttributes};
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    kwargs...)

    mapped_args = Vector{Pair{String, Any}}[]
    mapped_args = [mapped_args..., [Pair(string(j.first),j.second) for j = kwargs]...]

    if isa(opt, MOI.AbstractOptimizer)
        if length(mapped_args) > 0
            opt = MOI.OptimizerWithAttributes(typeof(opt), mapped_args...)
        else
            opt = typeof(opt)
        end
    end

    optimizer = MOI.instantiate(opt)

    if !isnothing(maxtime)
        MOI.set(optimizer, MOI.TimeLimitSec(), maxtime)
    end

    if !isnothing(reltol)
        @warn "common reltol argument is currently not used by $(optimizer). Set tolerances via optimizer specific keyword aguments."
    end

    if !isnothing(abstol)
        @warn "common abstol argument is currently not used by $(optimizer). Set tolerances via optimizer specific keyword aguments."
    end

    if !isnothing(maxiters)
        @warn "common maxiters argument is currently not used by $(optimizer). Set number of interations via optimizer specific keyword aguments."
    end

    return optimizer
end

function __solve(prob::OptimizationProblem, opt::Union{MOI.AbstractOptimizer, MOI.OptimizerWithAttributes};
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    kwargs...)

    maxiters = _check_and_convert_maxiters(maxiters)
    maxtime = _check_and_convert_maxtime(maxtime)

    opt_setup = __map_optimizer_args(prob, opt; abstol=abstol, reltol=reltol, maxiters=maxiters, maxtime=maxtime, kwargs...)

    num_variables = length(prob.u0)
	θ = MOI.add_variables(opt_setup, num_variables)
	if prob.lb !== nothing
        @assert eachindex(prob.lb) == Base.OneTo(num_variables)
		for i in 1:num_variables
			MOI.add_constraint(opt_setup, MOI.SingleVariable(θ[i]), MOI.GreaterThan(prob.lb[i]))
        end
    end
	if prob.ub !== nothing
        @assert eachindex(prob.ub) == Base.OneTo(num_variables)
		for i in 1:num_variables
			MOI.add_constraint(opt_setup, MOI.SingleVariable(θ[i]), MOI.LessThan(prob.ub[i]))
        end
    end
    @assert eachindex(prob.u0) == Base.OneTo(num_variables)
	for i in 1:num_variables
		MOI.set(opt_setup, MOI.VariablePrimalStart(), θ[i], prob.u0[i])
	end
    MOI.set(opt_setup, MOI.ObjectiveSense(), prob.sense === MaxSense ? MOI.MAX_SENSE : MOI.MIN_SENSE)
    if prob.lcons === nothing
        @assert prob.ucons === nothing
        con_bounds = MOI.NLPBoundsPair[]
    else
        @assert prob.ucons !== nothing
        con_bounds = MOI.NLPBoundsPair.(prob.lcons, prob.ucons)
    end
	MOI.set(opt_setup, MOI.NLPBlock(), MOI.NLPBlockData(con_bounds, make_moi_problem(prob), true))

	MOI.optimize!(opt_setup)

    if MOI.get(opt_setup, MOI.ResultCount()) >= 1
        minimizer = MOI.get(opt_setup, MOI.VariablePrimal(), θ)
        minimum = MOI.get(opt_setup, MOI.ObjectiveValue())
        opt_ret = Symbol(string(MOI.get(opt_setup, MOI.TerminationStatus())))
    else
        minimizer = fill(NaN, num_variables)
        minimum = NaN
        opt_ret= :Default
    end
    SciMLBase.build_solution(prob, opt, minimizer, minimum; original=opt_setup, retcode=opt_ret)
end

ext: jl


text:
@testset "MPS" begin

D = 10
d = 4
sites = 5
T = ComplexF64

@testset "Random MPS with the same physical dimension" begin

    ψ = randn(MPS{T}, sites, D, d)

    @testset "has correct number of sites" begin
        @test length(ψ) == sites
        @test size(ψ) == (sites, )
    end

    @testset "has correct type" begin
        @test eltype(ψ) == T
    end

    @testset "has correct rank" begin
        @test rank(ψ) == Tuple(fill(d, 1:sites))
    end

    @testset "has correct bonds" begin
        @test bond_dimension(ψ) ≈ D
        @test verify_bonds(ψ) === nothing
    end

    @testset "is equal to itself" begin
        @test ψ == ψ
        @test ψ ≈ ψ
    end

    @testset "is equal to its copy" begin
        ϕ = copy(ψ)
        @test ϕ == ψ
        @test ϕ ≈ ψ
    end
end

@testset "Random MPS with varying physical dimension" begin

    dims = (3, 2, 5, 4)
    ψ = randn(MPS{T}, D, dims)

    @testset "has correct number of sites" begin
        n = length(dims)
        @test length(ψ) == n
        @test size(ψ) == (n, )
    end

    @testset "has correct type" begin
        @test eltype(ψ) == T
    end

    @testset "has correct rank" begin
        @test rank(ψ) == dims
    end

    @testset "has correct bonds" begin
        @test bond_dimension(ψ) ≈ D
        @test verify_bonds(ψ) === nothing
    end

    @testset "is equal to itself" begin
        @test ψ == ψ
        @test ψ ≈ ψ
    end

    @testset "is equal to its copy" begin
        ϕ = copy(ψ)
        @test ϕ == ψ
        @test ϕ ≈ ψ
    end
end

@testset "Random MPO with the same physical dimension" begin

    W = randn(MPO{T}, sites, D, d)

    @testset "has correct number of sites" begin
        @test length(W) == sites
        @test size(W) == (sites, )
    end

    @testset "has correct type" begin
        @test eltype(W) == T
    end

    @testset "is equal to itself" begin
        @test W == W
        @test W ≈ W
    end

    @testset "is equal to its copy" begin
        U = copy(W)
        @test U == W
        @test U ≈ W
    end
end

@testset "Random MPO with varying physical dimension" begin

     dims = (3, 2, 5, 4)
     W = randn(MPO{T}, D, dims)

     @testset "has correct number of sites" begin
         n = length(dims)
         @test length(W) == n
         @test size(W) == (n, )
     end

     @testset "has correct type" begin
         @test eltype(W) == T
     end

     @testset "is equal to itself" begin
         @test W == W
         @test W ≈ W
     end

     @testset "is equal to its copy" begin
         U = copy(W)
         @test U == W
         @test U ≈ W
     end
end

@testset "MPS from tensor" begin
    ϵ = 1E-14

    dims = (2,3,4,3,5)
    sites = length(dims)
    A = randn(T, dims)

    ψ = MPS(A, :right)

    @test norm(ψ) ≈ 1
    @test_nowarn verify_bonds(ψ)
    @test_nowarn verify_physical_dims(ψ, dims)
    @test is_right_normalized(ψ)

    B = randn(T, dims...)
    ϕ = MPS(B, :left)

    @test norm(ϕ) ≈ 1
    @test_nowarn verify_bonds(ϕ)
    @test_nowarn verify_physical_dims(ϕ, dims)
    @test is_left_normalized(ϕ)

    χ = MPS(A, :left)

    @test norm(χ) ≈ 1
    @test abs(1 - abs(dot(ψ, χ))) < ϵ
end

end

@testset "Objects with equal tensors have the same hash" begin
    D = 10
    d = 4
    sites = 5
    T = ComplexF64

    ψ = randn(MPS{T}, sites, D, d)
    ϕ = copy(ψ)

    W = randn(MPO{T}, sites, D, d)
    V = copy(W)

    @testset "Equal MPSs have the same hash" begin
        @test hash(ψ) == hash(ϕ)
    end

    @testset "Equal MPOs have the same hash" begin
        @test hash(W) == hash(V)
    end

    @testset "Equal tuples with MPS and MPO have the same hash" begin
        tuple_1 = (ψ, W, [1, 2, 3])
        tuple_2 = (ϕ, V, [1, 2, 3])
        @test tuple_1 == tuple_2
        @test hash(tuple_1) == hash(tuple_2)
    end
end


text:
import Pkg;
#Pkg.add("Documenter")
using Documenter, PostgresORM

makedocs(
        sitename = "PostgresORM documentation",
        modules = [PostgresORM],
        pages = ["Index" => "index.md",
                 "Getting started" => "getting-started.md",
                 "Modules" => [
                        "modules/PostgresORM.md",
                        # "modules/PostgresORM.Controller.md"
                        ]
                 ],
)
deploydocs(; repo = "github.com/JuliaPostgresORM/PostgresORM.jl.git",
             devbranch = "main")


text:
#
# PAINTER.jl --
#
# Polychromatic opticAl INTErferometric Reconstruction software in Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of PAINTER.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015, Antony Schutz, Andre Ferrari.
#
#------------------------------------------------------------------------------
module PAINTER
    using JLD
    using OIFITS
    using OptimPack
    using NFFT
    using Wavelets
    using SpecialFunctions

    export painter, painterdemo, mask, paintersave, painterload, painterfitsexport, painterplotfct, WT, optiminit

    include("paintertype.jl")

    export PAINTER_Input, PAINTER_Data

    include("paintertools.jl")
    include("painterio.jl")
    include("painteroifits.jl")
    include("painterconstmat.jl")
    include("paintercheckinit.jl")
    include("painterplot.jl")
    include("optpckpt.jl")
end


text:
using GPUAcceleratedTracking, DrWatson, DataFrames
@quickactivate "GPUAcceleratedTracking"

raw_data_df = collect_results(datadir("benchmarks/reduction"))

num_ants = 4; num_correlators = 3;

fig_bar = plot_reduction_benchmark(raw_data_df, num_ants, num_correlators)
save( plotsdir("reduction.pdf"), fig_bar)

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
"""
    struct ElectricFieldChargeDriftModel{T <: SSDFloat} <: AbstractChargeDriftModel{T}

Charge drift model in which the electrons and holes drift along the electric field with a mobility of ± 1m²/Vs.

This model is the default when no charge drift model is defined in the configuration file.
"""
struct ElectricFieldChargeDriftModel{T <: SSDFloat} <: AbstractChargeDriftModel{T} end

ElectricFieldChargeDriftModel(T::Type{<:SSDFloat}) = ElectricFieldChargeDriftModel{T}()
getVe(fv::SVector{3, T}, cdm::ElectricFieldChargeDriftModel) where {T <: SSDFloat} = -fv
getVh(fv::SVector{3, T}, cdm::ElectricFieldChargeDriftModel) where {T <: SSDFloat} = fv


text:
# # Transmission line circuit
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/TransmissionLine.ipynb)
#
#md # !!! note "Overview"
#md #     System type: linear continuous system\
#md #     State dimension: parametric, typically between 4 to 40\
#md #     Application domain: Power Systems Stability
#
# We consider the transmission line model used in [AKS11].
# The RLC circuit of the transmission line is shown below.
# In the circuit, ``U_{in}`` is the voltage at the sending end,
# and ``U_{out}`` is the voltage at the receiving end.

# ![](https://raw.githubusercontent.com/JuliaReach/ReachabilityAnalysis.jl/master/examples/TransmissionLine/tline.png)

# ## Model description

# For reference to understand this model we refer to any standard textbook on electrical circuits.
# The electrical elements law for resistors (R), inductors (L) and capacitors (C) are studied,
# for instance, in Chapter 3 of [^K15].

# Let's assume that the network consists of ``η > 2`` nodes.
# There are ``η`` dynamic equations for the capacitor voltages and ``\eta`` dynamic
# equations for the circuit currents. Therefore, the state vector can be represented as

# ```math
# x = [U_1, U_2, …, U_η, I_1, I_2, …, I_η]^T,
# ```
# and the state dimension is ``2η``.

# **Analysis of the first node.** When writing the equations for the voltages and currents,
# we should pay attention to the sign choices in given circuit, which do not exactly match
# the convention in text-books (in the sense that ``U_1``s positive terminal is at the bottom,
# as indicated by the arrows). Let ``R_d`` denote the driver resistance's current.
# By Kirchhoff's voltage law,

# ```math
# U_{in} = R_d I_1 + LI_1' - U_1,
# ```
# therefore
# ```math
# \boxed{I_1' = \dfrac{U_{in} + U_1}{L} - \dfrac{R_d}{L}I_1}.
# ```
#
# By Kirchhoff's current law, and if ``I_{1, C}`` denotes the current through the
# capacitor ``C``, connected to the first node,
#
# ```math
#    I_1 = I_2 + I_{C, 1},
# ```
# and
# ```math
# I_{C, 1} = -CU_1',
# ```
# then
# ```math
#    \boxed{U_1' = \frac{I_2 - I_1}{C}}.
# ```

# **Analysis of the intermediate nodes.** Other node's equations are obtained
# in a similar fashion. For instance, for the second loop one has

# ```math
# -U_1 = RI_2 + LI_2' - U_2 \Rightarrow I_2' = \dfrac{U_2-U_1}{L} - \frac{R}{L}I_2
# ```
# for the current's equation, and
#
# ```math
#    I_{2, C} = I_2 - I_3,\Rightarrow U_2' = \frac{I_3 - I_2}{C}.
# ```
# for the voltage's equation.
#
# Generalizing for arbitrary ``l = 2, \ldots, \eta - 1`` is trivial and gives:
#
# ```math
#     \boxed{I'_l = \dfrac{U_l - U_{l-1}}{L} - \frac{R}{L}I_l} \\
#     \boxed{U'_l = \dfrac{I_{l+1} - I_l}{C}}.
# ```

# **Analysis of the last node.** The last node corresponds to the case ``l = \eta``
#
# ```math
# \boxed{I'_\eta = \dfrac{U_\eta - U_{\eta-1}}{L} - \frac{R}{L}I_\eta}, \\
# \boxed{U'_{out} = U'_{out} = - \frac{I_{\eta}}{C}}.
# ```

# **System of linear ODEs.** Recall that that the stateset is ``\mathbb{R}^{2\eta}``,
# where the state variables are
# ``x = [U_1, U_2, \ldots, U_\eta, I_1, I_2, \ldots, I_\eta]^T``.
# The system can be written as a block-diagonal system of linear ODEs,
# ```math
# x'(t) = Ax(t) + BU_{in}(t),
# ```
# using the results in the previous section.

# The coefficients matrix ``A`` can be written as a block-diagonal matrix.
# There exist useful constructors in the base package `LinearAlgebra` that greatly simplify
# building matrices with special shape, as in our case, such as diagonal and band matrices,
# using the types `Diagonal` and `Bidiagonal` in our case.

using ReachabilityAnalysis, LinearAlgebra, SparseArrays

function tline(;η=3, R=1.00, Rd=10.0, L=1e-10, C=1e-13 * 4.00)
    A₁₁ = zeros(η, η)
    A₁₂ = Bidiagonal(fill(-1/C, η), fill(1/C, η-1), :U)
    A₂₁ = Bidiagonal(fill(1/L, η), fill(-1/L, η-1), :L)
    A₂₂ = Diagonal(vcat(-Rd/L, fill(-R/L, η-1)))
    A  = [A₁₁ A₁₂; A₂₁ A₂₂]
    B = sparse([η+1], [1], 1/L, 2η, 1)
    return A, B
end

# We can visualize the structure of the cofficients matrix ``A`` for the case
# $\eta = 20$ with `spy` plot:

using Plots

A, _ = tline(η=20)
spy(A, legend=nothing, markersize=2.0, title="Sparsity pattern of A", xlab="columns", ylab="rows")

# Notice that the matrix coefficients are rather big, hence it is convenient to rescale the system
# for numerical stability.
# Let ``α > 0`` be a scaling factor, and let ``\tilde{x}(t) = x(\alpha t)``.
# By the chain rule,
# ```math
# \tilde{x}'(t) = \alpha x'(\alpha t) = \alpha A x(\alpha t) + \alpha B U_{in}(\alpha t) = \tilde{A} \tilde{x}(t) + \tilde{B} \tilde{U}_{in}(t),
# ```
# where ``\tilde{A} := \alpha A`` and ``\tilde{B} := \alpha B``.

function scale!(s, α=1.0)
    s.A .*= α
    s.B .*= α
    return s
end

# Note that under this transformation, the time horizon has to be transformed as well,
# through ``\tilde{T} = \alpha T``.

# ## Reachability settings

# The transmission line parameters used in this model are displayed in the following table.

# |resistance in [Ω] | driver resistance in [Ω] | Inductance in [H] | Capacitance in [F] |
# |------------------|--------------------------|------------------|-------------------|
# |R = 1.00| Rdriver = 10.0| L = 1e−10 | C = 4e−13|

# The steady state is obtained by zeroing the left-hand side of the ODE, which gives
#
# ```math
#    0 = x' = Ax_\infty + Bu_0 \Rightarrow x_\infty = -A^{-1}B u_0
# ```
# The initial set under consideration corresponds to the steady state for input voltages
# ``U_{in, ss} := [-0.2, 0.2]``.
# Moreover, an uncertainty is added so that the initial currents are also uncertain.
# The set of initial states is then
# ```math
#    x(0) \in \mathcal{X}_0 := -A^{-1} B  U_{in, ss} \oplus □(0.001),
# ```
# where ``□(ϵ)`` is the infinity-norm ball of center zero and radius ``ϵ``.

# The time horizon is chosen as ``T = 0.7`` seconds. We consider a
# scaling factor ``\alpha = 1e-9``.

# We are interested in the step response to an input voltage ``U_{in}(t)``, constant
# for ``t \in [0, T]`` over the domain ``U_{in} = [0.99, 1.01]``.

# We consider the case of ``\eta = 20`` nodes as in [AKS11], such that the system
# has ``n = 40`` state variables.

η = 20 # order
n = 2η # state-space dimension
A, B = tline(η=η);

Uin_ss = Interval(-0.2, 0.2)
□(ϵ) = BallInf(zeros(n), ϵ)
X0 = -inv(Matrix(A)) * B * Uin_ss ⊕ □(0.001);

Uin = Interval(0.99, 1.01)
s = @system(x'= A*x + B*u, x ∈ Universe(n), u ∈ Uin)
α = 1e-9 # scaling factor
scale!(s, α);

T = 0.7 * 1e-11 # time horizon
P = InitialValueProblem(s, X0);

# ## Results

# We solve system using a step size of `δ=1e-3` and the box algorithm.

sol = solve(P, T=0.7, alg=BOX(δ=1e-3));

# To get the variable ``U_{out}`` we have to project onto the ``η``-th coordinate
# and invert the sign of the flowpipe.

Uout_vs_t = @. (-1.0) * project(sol, η);

#-

plot(Uout_vs_t, vars=(0, η), color=:blue, xlab="t", ylab="Uout", alpha=.5, lw=0.5)

#md # !!! tip "Technical note"
#md #     Since we are only interested in the behavior of ``U_{out}``, we can use
#md #     algorithm `BFFPSV18` with the options `alg=BFFPSV18(δ=1e-3, dim=statedim(P), vars=[η]))`,
#md #     which will use an interval (1D) decomposition of the state space and only compute
#md #     the flowpipe associated with variable ``\eta``.

# ## References

# [^AKS11]: Althoff, Matthias, Bruce H. Krogh, and Olaf Stursberg. [Analyzing reachability of linear dynamic systems with parametric uncertainties.](http://www.i6.in.tum.de/Main/Publications/Althoff2011b.pdf). Modeling, Design, and Simulation of Systems with Uncertainties. Springer, Berlin, Heidelberg, 2011. 69-94.

# [^K15]: Kluever, Craig A. [Dynamic systems: modeling, simulation, and control.](https://www.wiley.com/en-us/Dynamic+Systems%3A+Modeling%2C+Simulation%2C+and+Control-p-9781118289457). John Wiley & Sons, 2015.

# [^GLG08]: Girard, Antoine, and Colas Le Guernic. [Efficient reachability analysis for linear systems using support functions.]() IFAC Proceedings Volumes 41.2 (2008): 8966-8971.


text:
using BinDeps
using Compat.Libdl
using Base.Sys: WORD_SIZE
import Compat.Sys
@BinDeps.setup

function validate_lz4(name,handle)
    f = Libdl.dlsym_e(handle, "LZ4F_getVersion")
    return f != C_NULL
end

liblz4 = library_dependency("liblz4", validate = validate_lz4)

long_version = "1.8.1.2"
short_version = "1.8.1"
zipname = "lz4_v$(replace(short_version, "." => "_"))_win$(WORD_SIZE).zip"

suffix = "so.$short_version"
if Sys.isapple()
    suffix = "$short_version.$(Libdl.dlext)"
end

# Best practice to use a fixed long_version here, either a long_version number tag or a git sha
# Please don't download "latest master" because the version that works today might not work tomorrow

provides(Sources, URI("https://github.com/lz4/lz4/archive/v$long_version.tar.gz"),
    liblz4, unpacked_dir="lz4-$long_version", os = :Unix)

srcdir = joinpath(BinDeps.srcdir(liblz4), "lz4-$long_version")

provides(SimpleBuild,
    (@build_steps begin
        GetSources(liblz4)
        CreateDirectory(BinDeps.libdir(liblz4))
        @build_steps begin
            ChangeDirectory(srcdir)
            MAKE_CMD
            `mv lib/liblz4.$suffix "$(BinDeps.libdir(liblz4))/liblz4.$(Libdl.dlext)"`
        end
    end), liblz4, os = :Unix)

provides(BuildProcess,
    (@build_steps begin
        FileDownloader("https://github.com/lz4/lz4/releases/download/v$long_version/$zipname",
        joinpath(BinDeps.downloadsdir(liblz4), zipname))
        CreateDirectory(BinDeps.srcdir(liblz4), true)
        FileUnpacker(joinpath(BinDeps.downloadsdir(liblz4), zipname), BinDeps.srcdir(liblz4), "dll")
        CreateDirectory(BinDeps.libdir(liblz4), true)
        @build_steps begin
            ChangeDirectory(joinpath(BinDeps.srcdir(liblz4), "dll"))
            `powershell cp "liblz4.$suffix.$(Libdl.dlext)" $(joinpath(BinDeps.libdir(liblz4),"liblz4.$(Libdl.dlext)"))`
        end
    end), liblz4, os = :Windows)

if Sys.iswindows()
    push!(BinDeps.defaults, BuildProcess)
    @BinDeps.install Dict(:liblz4 => :liblz4)
    pop!(BinDeps.defaults)
else
    @BinDeps.install Dict(:liblz4 => :liblz4)
end


text:
"""
$(DocStringExtensions.README)
"""
module DiffEqBayes

using DocStringExtensions
using DiffEqBase, Distributions, Turing, MacroTools
using RecursiveArrayTools, ModelingToolkit
using Parameters, Distributions, Optim, Requires
using Distances, DocStringExtensions, Random, StanSample

STANDARD_PROB_GENERATOR(prob,p) = remake(prob;u0=eltype(p).(prob.u0),p=p)
STANDARD_PROB_GENERATOR(prob::EnsembleProblem,p) = EnsembleProblem(remake(prob.prob;u0=eltype(p).(prob.prob.u0),p=p))

include("turing_inference.jl")
# include("abc_inference.jl")
include("stan_string.jl")
include("stan_inference.jl")

function __init__()
    @require DynamicHMC="bbc10e6e-7c05-544b-b16e-64fede858acb" begin
        using .DynamicHMC, TransformVariables, LogDensityProblems
        include("dynamichmc_inference.jl")
        export dynamichmc_inference
    end
end

export turing_inference, stan_inference ,abc_inference
end # module


text:
"""
    y = delay(x, n=1)

Delay signal `x` by `n` samples
"""
function delay(x, n=1)
    nx = length(x)
    if nx <= n
        return zeros(size(x))
    else
        return [zeros(n); x[1:nx-n]]
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
#import Winston

@testset "PostProcessTests" begin
  @testset "refine_valid_triangles" begin
    local mesh = DelaunayMeshes.Mesh()
    DelaunayMeshes.setboundingbox(mesh, [-15.0, 15.0, -15.0, 15.0])
    # seed rng and store seed in file
    local seed = 2628265373
    #local seed = rand(UInt32)
    local fs = open("refine_valid_triangles.seed", "w")
    write(fs, string(seed))
    close(fs)

    # push some random vertices
    srand(seed)
    local nvertices = 200
    local points = rand(Float64, nvertices, 2)*30. - 15.
    push!(mesh, points)

    local nringvertices = 40

    # push inner ring constraint vertices
    local θ = linspace(0., 2.*π, nringvertices + 1)[1:end-1]
    local rinner = 5.
    local router = 10.
    local innerRing = rinner * [cos(θ'); sin(θ')]'
    local outerRing = router * [cos(θ'); sin(θ')]'

    push!(mesh, innerRing)
    push!(mesh, outerRing)

    # add constraints
    local innerConstraintVertexList = [x for x in (nvertices+3+nringvertices):-1:(nvertices+3+1)]
    local outerConstraintVertexList = [x for x in (nvertices+3+nringvertices+1):(nvertices + 3 + 2*nringvertices)]
    DelaunayMeshes.addconstraint!(mesh, innerConstraintVertexList)
    DelaunayMeshes.addconstraint!(mesh, outerConstraintVertexList)

    # refine
    DelaunayMeshes.refine_valid_triangles!(mesh)

    # check face locations
    local innerBoundaryEdge = findedgeconnectingvertices(mesh.tesselation, nvertices+3+1, nvertices+3+2)
    @test mesh.tesselation.edges[innerBoundaryEdge].constraintType == DelaunayMeshes.OnLeftSide
    local insideFaces = DelaunayMeshes.getfacesinsideregion(mesh, DelaunayMeshes.rot(innerBoundaryEdge))
    local expLoc = fill(false, length(mesh.tesselation.faces))
    expLoc[insideFaces] = true
    @test expLoc == mesh.faceLocation

    # check conformity of triangles
    check_conformity(mesh.tesselation)

    # check that all triangles are high quality
    for tri in mesh.tesselation
      if (mesh.faceLocation[tri.face])
        @test DelaunayMeshes.check_triangle(mesh, tri)
      end
    end

    # usual checks
    testvertexcache(mesh.tesselation)

    #=
    # plot it
    xc, yc = DelaunayMeshes.getdelaunaycoordinates(mesh.tesselation)
    p = Winston.FramedPlot(aspect_ratio=1
      #,xrange=[1.496 1.50], yrange=[1.256, 1.26]
      )
    Winston.add(p, Winston.Curve(xc, yc))
    Winston.savefig(p, "refine_valid_triangles.svg")
    =#

  end

  @testset "refine_triangle" begin
    local mesh = DelaunayMeshes.Mesh()
    DelaunayMeshes.setboundingbox(mesh, [-15.0, 15.0, -15.0, 15.0])

    push!(mesh, [-5.0 -5.0; 5.0 -5.0; 0.0 5.0])

    # find face foce for vertices
    local testTriangle = DelaunayMeshes.Triangle(mesh.tesselation, 4, 5, 6)

    # refine
    DelaunayMeshes.refine_triangle(mesh, testTriangle)

    # check that there is a new vertex and that it's at the offcenter
    @test length(mesh.tesselation.vertices) == 7
    local offcenter = DelaunayMeshes.compute_offcenter(mesh, 4, 5, 6)
    @test_approx_eq(mesh.tesselation.vertices[7]', offcenter')

    #=
    # plot it
    xc, yc = getplotedges(mesh.tesselation)
    local vertexPoints = reduce(hcat, mesh.tesselation.vertices)'
    p = Winston.FramedPlot(aspect_ratio=1
      ,xrange=[1.45, 1.55], yrange=[1.2, 1.3]
      )
    Winston.add(p, Winston.Points(vertexPoints[:,1], vertexPoints[:,2], kind="circle", color="red"))
    Winston.add(p, Winston.Curve(xc, yc))
    Winston.savefig(p, "refine_triangle.svg")
    =#
  end

  @testset "refine_triangle_constraints" begin
    local mesh = DelaunayMeshes.Mesh()
    DelaunayMeshes.setboundingbox(mesh, [-15.0, 15.0, -15.0, 15.0])

    push!(mesh, [-5.0 -5.0; 5.0 -5.0; 0.0 5.0])

    # find face foce for vertices
    local tri = DelaunayMeshes.Triangle(mesh.tesselation, 4, 5, 6)

    # set edges as constraints
    DelaunayMeshes.setconstraint!(mesh.tesselation, tri.edges[1], DelaunayMeshes.OnLeftSide)
    DelaunayMeshes.setconstraint!(mesh.tesselation, tri.edges[2], DelaunayMeshes.OnLeftSide)

    # refine
    DelaunayMeshes.refine_triangle(mesh, tri)

    # check that there is a new vertex and that it's at the edge centers
    @test length(mesh.tesselation.vertices) == 8
    midpoint1 = 0.5*(mesh.tesselation.vertices[4] + mesh.tesselation.vertices[5])
    midpoint2 = 0.5*(mesh.tesselation.vertices[5] + mesh.tesselation.vertices[6])
    @test_approx_eq(mesh.tesselation.vertices[7]', midpoint1')
    @test_approx_eq(mesh.tesselation.vertices[8]', midpoint2')

    # check that the edges going to and from these vertices are marked as
    # segments
    local e47 = findedgeconnectingvertices(mesh.tesselation, 4, 7)
    local e75 = findedgeconnectingvertices(mesh.tesselation, 7, 5)
    @test mesh.tesselation.edges[e47].constraintType == DelaunayMeshes.OnLeftSide
    @test mesh.tesselation.edges[e75].constraintType == DelaunayMeshes.OnLeftSide
    @test mesh.tesselation.edges[DelaunayMeshes.sym(e47)].constraintType == DelaunayMeshes.OnRightSide
    @test mesh.tesselation.edges[DelaunayMeshes.sym(e75)].constraintType == DelaunayMeshes.OnRightSide
    local e58 = findedgeconnectingvertices(mesh.tesselation, 5, 8)
    local e86 = findedgeconnectingvertices(mesh.tesselation, 8, 6)
    @test mesh.tesselation.edges[e58].constraintType == DelaunayMeshes.OnLeftSide
    @test mesh.tesselation.edges[e86].constraintType == DelaunayMeshes.OnLeftSide
    @test mesh.tesselation.edges[DelaunayMeshes.sym(e58)].constraintType == DelaunayMeshes.OnRightSide
    @test mesh.tesselation.edges[DelaunayMeshes.sym(e86)].constraintType == DelaunayMeshes.OnRightSide

    #=
    # plot it
    xc, yc = getplotedges(mesh.tesselation)
    local vertexPoints = reduce(hcat, mesh.tesselation.vertices)'
    p = Winston.FramedPlot(aspect_ratio=1
      ,xrange=[1.45, 1.55], yrange=[1.2, 1.3]
      )
    Winston.add(p, Winston.Points(vertexPoints[:,1], vertexPoints[:,2], kind="circle", color="red"))
    Winston.add(p, Winston.Curve(xc, yc))
    Winston.savefig(p, "refine_triangle.svg")
    =#
  end

  @testset "get_shortest_edge" for par in [[-1.0, 2.0, 1.3, 5.0, 7.0, 7.5], [1.0, -2.0, 3.141, 10.0, 12.0, 21.5], [0.0, 0.0, 0.0, 1.0, 2.0, 2.5]] begin
    cx, cy, alpha0, a, b, c = par

    mesh = DelaunayMeshes.Mesh()

    s = 0.5*(a+b+c)
    r = sqrt((s-a)*(s-b)*(s-c)/s)
    alpha = 2.0*(atan(r/(s-a)))
    beta = 2.0*(atan(r/(s-b)))
    gamma = 2.0*(atan(r/(s-c)))
    (v1x, v1y) = (cx + b * cos(alpha + alpha0), cy + b * sin(alpha + alpha0))
    (v2x, v2y) = (cx + c * cos(alpha0), cy + c * sin(alpha0))
    mesh.tesselation.vertices = [[cx, cy],[v2x, v2y],[v1x, v1y]]

    p, q, r = DelaunayMeshes.get_shortest_edge(mesh, 1, 2, 3)
    @test p == 1
    @test q == 2
    @test r == 3
    p, q, r = DelaunayMeshes.get_shortest_edge(mesh, 2, 3, 1)
    @test p == 1
    @test q == 2
    @test r == 3
    p, q, r = DelaunayMeshes.get_shortest_edge(mesh, 3, 1, 2)
    @test p == 1
    @test q == 2
    @test r == 3
  end
  end

  @testset "compute_circumcenter_radius" for r in [0.1, 1.1], c in [[1.6, 1.3], [1.25, 1.7]], t in [[0, pi/2., pi], [pi/8, pi-pi/16, pi+pi/16]] begin
    # generate triangle on given circumcircle
    local cx = c[1]
    local cy = c[2]
    local t1 = t[1]
    local t2 = t[2]
    local t3 = t[3]
    mesh = create_mesh_with_triangle(cx, cy, r, t1, t2, t3)

    sx, sy, rs = DelaunayMeshes.compute_circumcenter_radius(mesh, 1, 2, 3)

    @test(abs(sx-cx) <= 1e-4)
    @test(abs(sy-cy) <= 1e-4)
    @test(abs(rs-r*r) <= 1e-4)
  end
  end

  @testset "compute_offcenter" for c in [[1.6, 1.3], [1.25, 1.7]] , a0 in [pi/24., pi/16., pi/8.], offset in [0., pi/3., pi, 1.5*pi] begin
    local cx = c[1]
    local cy = c[2]
    local r = 2.
    local t1 = pi - a0 + offset
    local t2 = pi + a0 + offset
    local t3 = pi/4. + offset
    mesh = create_mesh_with_triangle(cx, cy, r, t1, t2, t3)

    local offcenter = DelaunayMeshes.compute_offcenter(mesh, 1, 2, 3)

    # check that new triangle p - q - oc has correct edge/circumcircle radius
    mesh.tesselation.vertices[3] = offcenter
    sx, sy, rs = DelaunayMeshes.compute_circumcenter_radius(mesh, 1, 2, 3)
    local lpq = norm(mesh.tesselation.vertices[1] - mesh.tesselation.vertices[2])
    @test sqrt(rs)/lpq <= DelaunayMeshes.beta
  end
  end

  @testset "vertex_encroaches_segment" for rd in [0.8, 1.2] begin
    # seed rng and store seed in file
    #local seed =2410381395
    local seed = rand(UInt32)
    local fs = open("vertex_encroaches_segment.seed", "w")
    write(fs, string(seed))
    close(fs)
    srand(seed)

    mesh = DelaunayMeshes.Mesh()

    # mark first edge as boundary
    DelaunayMeshes.setconstraint!(mesh.tesselation, 1, DelaunayMeshes.OnLeftSide)
    mesh.tesselation.edges[5].constraintType = DelaunayMeshes.NoConstraint
    mesh.tesselation.edges[7].constraintType = DelaunayMeshes.NoConstraint

    # create potentially encroaching vertex for that edge
    me = 0.5*(mesh.tesselation.vertices[1] + mesh.tesselation.vertices[2])
    l = norm(mesh.tesselation.vertices[1] - mesh.tesselation.vertices[2])
    θ = rand(Float64) * 2. * pi
    vp = convert(DelaunayMeshes.Point, me + rd * 0.5 * l * [cos(θ), sin(θ)])
    @test DelaunayMeshes.vertex_encroaches_segment(mesh, 1, vp) == (rd < 1.0)
    @test DelaunayMeshes.vertex_encroaches_segment(mesh, 3, vp) == (rd < 1.0)

    # create possibly encroaching vertex for non boundary
    me = 0.5*(mesh.tesselation.vertices[2] + mesh.tesselation.vertices[3])
    l = norm(mesh.tesselation.vertices[2] - mesh.tesselation.vertices[3])
    θ = rand(Float64) * 2. * pi
    vp = convert(DelaunayMeshes.Point, me + rd * 0.5 * l * [cos(θ), sin(θ)])
    @test DelaunayMeshes.vertex_encroaches_segment(mesh, 5, vp) == false
    @test DelaunayMeshes.vertex_encroaches_segment(mesh, 7, vp) == false
  end
  end
end


text:
printstyled("  SAC support extensions\n", color=:light_green)

printstyled("    writesac from SeisEvent\n", color=:light_green)
writesac(rse_wb(8))

sac_file    = path*"/SampleFiles/SAC/test_le.sac"
f_stub      = "1981.088.10.38.23.460"
f_out       = f_stub * "..CDV...R.SAC"
sacv7_out   = "v7_out.sac"

test_fs   = 50.0
test_lat  = 41.5764
test_lon  = -87.6194
test_ot   = -0.75
test_stla = 41.54
test_stlo = -87.64
test_mag  = 3.1f0

# Read test file
C = verified_read_data("sac", sac_file, full=true)[1]
ev_lat = C.misc["evla"]
ev_lon = C.misc["evlo"]
ev_dep = C.misc["evdp"]
ev_id  = "8108838" # taken from kevnm

Ev = SeisEvent()
Ev.data = convert(EventTraceData, SeisData(C))
Ev.hdr.mag.val = test_mag
Ev.hdr.loc.dep = ev_dep
Ev.hdr.loc.lat = ev_lat
Ev.hdr.loc.lon = ev_lon
Ev.hdr.id = ev_id

# Write to file
writesac(Ev)
@test safe_isfile(f_out)

# read unmodified file, check for preserved values
C = verified_read_data("sac", f_out, full=true)[1]
@test C.misc["evla"] == ev_lat
@test C.misc["evlo"] == ev_lon
@test C.misc["evdp"] == ev_dep
@test C.misc["mag"] == test_mag
@test string(C.misc["nevid"]) == ev_id

# Modify the original file
printstyled("    SAC v7\n", color=:light_green)
io = open(sac_file, "r")
sac_raw = read(io)
close(io)

# Set version to 7
sac_raw[305] = 0x07

# Set magnitude to test_mag
mag = reinterpret(UInt8, [test_mag])
sac_raw[157:160] .= mag

# Change some values
reset_sacbuf()
dv     = BUF.sac_dv
dv[1]  = 1.0/test_fs
dv[4]  = test_ot
dv[17] = test_lon
dv[18] = test_lat
dv[19] = test_stla
dv[20] = test_stlo
dv2 = deepcopy(dv)
sac_dbl_buf = reinterpret(UInt8, dv)
io = open(sacv7_out, "w")
write(io, sac_raw)
write(io, sac_dbl_buf)
close(io)

C = read_data("sac", sacv7_out, full=true)[1]
@test C.fs == test_fs

printstyled("    fill_sac_evh!\n", color=:light_green)
Ev = SeisEvent()
Ev.data = convert(EventTraceData, SeisData(C))

fill_sac_evh!(Ev, sacv7_out, k=1)
@test Ev.hdr.loc.lat == test_lat
@test Ev.hdr.loc.lon == test_lon
@test Ev.hdr.mag.val == test_mag
@test isapprox(d2u(Ev.hdr.ot) - Ev.data.t[1][1,2]*μs, test_ot, atol=0.001)

printstyled("      big-endian\n", color=:light_green)
sac_be_file = path*"/SampleFiles/SAC/test_be.sac"
io = open(sac_be_file, "r")
sac_raw = read(io)
close(io)
sac_raw[308] = 0x07

reset_sacbuf()
dv .= bswap.(dv2)
sac_dbl_buf = reinterpret(UInt8, dv)

io = open(sacv7_out, "w")
write(io, sac_raw)
write(io, sac_dbl_buf)
close(io)

C = read_data("sac", sacv7_out, full=true)[1]
Ev = SeisEvent()
Ev.data = convert(EventTraceData, SeisData(C))
fill_sac_evh!(Ev, sacv7_out, k=1)
@test Ev.hdr.loc.lat == test_lat
@test Ev.hdr.loc.lon == test_lon
@test isapprox(d2u(Ev.hdr.ot) - Ev.data.t[1][1,2]*μs, test_ot, atol=0.001)


text:
import ACE
using LinearAlgebra, StaticArrays, BenchmarkTools, Test, Printf
using ACE: evaluate, evaluate_d!, evaluate_ed, alloc_B, release!,
           evaluate!, acquire!, acquire_B!, acquire_dB!,
           release_B!, release_dB!, evaluate_d
using ACE.Transforms: PolyTransform
using ACE.Testing
using BenchmarkTools

##

trans = PolyTransform(2, 1.0)
P = ACE.OrthPolys.transformed_jacobi(20, trans, 2.0, 0.5; pcut = 2, pin = 2)

B = acquire_B!(P)
dB = acquire_dB!(P)
r = ACE.rand_radial(P)

@info("non-allocating (for reference)")
@btime evaluate!($B, $P, $r)
@btime evaluate_d!($dB, $P, $r)

@info("fully allocating")
@btime evaluate($P, $r)
@btime evaluate_d($P, $r)

@info("allocate and release")
@info("In this case there seems to be limited advantage?")
_ev(P, r) = (B = evaluate!(acquire_B!(P, r), P, r); release_B!(P, B))
_ev_d(P, r) = (dB = evaluate_d!(acquire_dB!(P, r), P, r); release_dB!(P, dB))
@btime _ev($P, $r)
@btime _ev_d($P, $r)


text:
module RulesBIFMHelperOutTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "rules:BIFMHelper:out" begin

     @testset "Variational Message Passing: (q_in::Any, )" begin

        @test_rules [ with_float_conversions = true ] BIFMHelper(:out, Marginalisation) [
            (input = (q_in = MvNormalMeanCovariance([1.0, 2.0], [3.0 0; 0 2.0]), ),         output = ProdFinal(MvNormalMeanCovariance([1.0, 2.0], [3.0 0; 0 2.0])) ),
            (input = (q_in = NormalMeanVariance(9.5, 3.2), ),                               output = ProdFinal(NormalMeanVariance(9.5, 3.2)) ),
            (input = (q_in = MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 0; 0 2.0]), ),  output = ProdFinal(MvNormalWeightedMeanPrecision([1.0, 2.0], [3.0 0; 0 2.0])) )
        ]

    end

end

end


text:
export SimulatedTemperatureSensor, SimulatedTemperatureSensorParams, numChannels, getTemperature

Base.@kwdef struct SimulatedTemperatureSensorParams <: DeviceParams

end
SimulatedTemperatureSensorParams(dict::Dict) = params_from_dict(SimulatedTemperatureSensorParams, dict)

Base.@kwdef mutable struct SimulatedTemperatureSensor <: TemperatureSensor
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::SimulatedTemperatureSensorParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}
end

init(sensor::SimulatedTemperatureSensor) = sensor.present = true
neededDependencies(::SimulatedTemperatureSensor) = []
optionalDependencies(::SimulatedTemperatureSensor) = []
Base.close(sensor::SimulatedTemperatureSensor) = nothing

numChannels(sensor::SimulatedTemperatureSensor) = 1
getTemperature(sensor::SimulatedTemperatureSensor)::Vector{typeof(1u"°C")} = [42u"°C"]
getTemperature(sensor::SimulatedTemperatureSensor, channel::Int)::typeof(1u"°C") = 42u"°C"


text:
#io = open("C:/Users/Andrew.Bannerman/Desktop/Julia/My_First_Neural_Network/Data/train-images-idx3-ubyte/train-images.idx3-ubyte", n=28)
#write(io)

#open("C:/Users/Andrew.Bannerman/Desktop/Julia/My_First_Neural_Network/Data/train-images-idx3-ubyte/train-images.idx3-ubyte") do file
#    for ln in eachline(file)
#        println("$(length(ln)), $(ln)")
#    end
#end

# Back propogation neural network
# 'Hello World' of neural networks
# Andrew Bannerman 4.21.2019

# Resources
# Cheat Sheet = https://ml-cheatsheet.readthedocs.io/en/latest/forwardpropagation.html
# Matrix multiplication Rules - https://www.mathsisfun.com/algebra/matrix-multiplying.html
# How to get MNIST data as a CSV https://pjreddie.com/projects/mnist-in-csv/
# Kaggle - https://www.kaggle.com/ngbolin/mnist-dataset-digit-recognizer
# drop out = https://wiseodd.github.io/techblog/2016/06/25/dropout/
# Neural network: https://github.com/wiseodd/hipsternet/blob/master/hipsternet/neuralnet.py

# To do
# root mena prop: https://engmrk.com/rmsprop/?utm_campaign=News&utm_medium=Community&utm_source=DataCamp.com
# softmax
# numpy from scratch: https://www.kaggle.com/scaomath/simple-mnist-numpy-from-scratch
# Activation function == https://github.com/Kulbear/deep-learning-nano-foundation/wiki/ReLU-and-Softmax-Activation-Functions
# MATH AI = https://hackernoon.com/learning-ai-if-you-suck-at-math-p6-math-notation-made-easy-1277d76a1fe5

# Machine learning cheat sheet: https://ml-cheatsheet.readthedocs.io/en/latest/activation_functions.html
# Drop out video: https://www.youtube.com/watch?v=UcKPdAM8cnI /
# Deeplearning.ai
# https://www.youtube.com/watch?v=D8PJAL-MZv8

# Bias = https://medium.com/@erikhallstrm/backpropagation-from-the-beginning-77356edf427d
# Full code for updating bias: https://medium.com/coinmonks/implement-back-propagation-in-neural-networks-ed09897593e7 #--- numpy

# To Do
# 4.22.19 - drop out
# 4.22.19 - bias

# Needed packages

using DataFrames
using CSV
using Distributions
using Plots
using Random
using Images, CoordinateTransformations, TestImages, OffsetArrays, CoordinateTransformations
using PaddedViews

lineplot(sin, 1:.5:20, labels=false)

# Train and test set data preparation
# Load MNIST Data set from .csv format
train_set = CSV.read("C:/Users/Andrew.Bannerman/Desktop/Julia/First_Neural_Network.jl/kaggle/MNIST/train.csv",header=true,types=fill(Float64,785)) # "C:/Users/Andrew.Bannerman/Desktop/Julia/My_First_Neural_Network/Data/mnist_train.csv"
test_set = CSV.read("C:/Users/Andrew.Bannerman/Desktop/Julia/First_Neural_Network.jl/kaggle/MNIST/test.csv",header=true,types=fill(Float64,784)) # "C:/Users/Andrew.Bannerman/Desktop/Julia/My_First_Neural_Network/Data/mnist_test.csv"

# Plot Example image
example = Float64.(train_set[4,2:size(train_set,2)])
example = reshape(example,28,28)
example = rotl90(reshape(example,28,28))
heatmap(example)

# Extract labels
train_label = train_set[:,1]
test_label = test_set[:,1]

# Extract Data and exclude the label in the first column
train_set = train_set[:,2:size(train_set,2)]
test_set = test_set[:,1:size(test_set,2)]

# Pre process - Scale the inputs to a range of 0.01 to 1.00
m_train = Array{Float64}(zeros(size(train_set,1),size(train_set,2))) # empty matrix
@inbounds for i in 1:size(train_set,1)
    @inbounds for j in 1:size(train_set,2)
    m_train[i,j] = (train_set[i,j] / 255.0 * 0.99) + 0.01
    end
end

m_test = Array{Float64}(zeros(size(test_set,1),size(test_set,2))) # empty matrix
@inbounds for i in 1:size(test_set,1)
    @inbounds for j in 1:size(test_set,2)
    m_test[i,j] = (test_set[i,j] / 255.0 * 0.99) + 0.01
    end
end

# Convert matrix to DataFrame and append the labels back to the train and test set
scaled_orig_train_inputs = hcat(train_label,m_train)
scaled_test_inputs = DataFrame(hcat(m_test))

# Plot the new scale example
example = Float64.(scaled_test_inputs[7,1:size(scaled_test_inputs,2)])
example = rotl90(reshape(example,28,28))
heatmap(example)

# Data Augmentation
# Rotate images
training_set_aug_minus = zeros(size(scaled_orig_train_inputs,1),size(scaled_orig_train_inputs,2))
training_set_aug_plus = zeros(size(scaled_orig_train_inputs,1),size(scaled_orig_train_inputs,2))
degree_rot_minus = -.3
degree_rot_plus = .3
i=1
for i in 1:size(scaled_orig_train_inputs,1)
        label = scaled_orig_train_inputs[i,1]
        example = Float64.(scaled_orig_train_inputs[i,2:size(scaled_orig_train_inputs,2)])
        example = reshape(example,28,28)
        #rot_minus = LinearMap(RotMatrix(degree_rot_minus))
        #rot_plus = LinearMap(RotMatrix(degree_rot_plus))
        #imgw_minus = warp(example_minus, rot_minus, axes(example_plus)) #  axes(example_minus) = crops
        #imgw_plus = warp(example_plus, rot_plus, axes(example_plus)) # axes() = crops
        imgw_minus = imrotate(example, degree_rot_minus, axes(example))
        imgw_plus = imrotate(example, degree_rot_plus, axes(example))
        imgw_minus = reshape(imgw_minus,1,784)
        imgw_plus = reshape(imgw_plus,1,784)
            for k in 1:size(imgw_minus,2)
                if isnan(imgw_minus[k]) # Convert NaN to 0.0
                    imgw_minus[k] = 0.01
                end
                if isnan(imgw_plus[k]) # Convert NaN to 0.0
                    imgw_plus[k] = 0.01
                end
            end
            #imgw_plus .= (imgw_plus ./ 255.0 .* 0.99) .+ 0.01 # Scale the outputs
            #imgw_minus .= (imgw_minus ./ 255.0 .* 0.99) .+ 0.01 # Scale the outputs
            training_set_aug_minus[i,1] = label
            training_set_aug_minus[i,2:size(scaled_orig_train_inputs,2)] = imgw_minus
            training_set_aug_plus[i,1] = label
            training_set_aug_plus[i,2:size(scaled_orig_train_inputs,2)] = imgw_plus
        print("This is iteration i",i,"\n")
end

    maximum(imgw_plus)

# Combine augmented data with the original
# note may want to add shuffling?
scaled_train_inputs = DataFrame(vcat(scaled_orig_train_inputs,training_set_aug_minus,training_set_aug_plus))
#scaled_train_inputs = DataFrame(training_set_aug_plus)

# Shuffle the augmented + original data
scaled_train_inputs = scaled_train_inputs[shuffle(1:end), :]

# Build the back propogation neural network

# Create sigmoid activation function
function sigmoid(x)
    return 1 / (1.0 + exp(-x))
end
#return maximum(0.0,x)
function ReLU(x)
    #return x*(x>0)
    return maximum([0.0,x])
end

x = -500
ReLU(x)

"""
Leaky ReLU
f(x)=1(x<0)(αx)+1(x>=0)(x) where α = some constant
"""
function ReLuL(x)
    if (x > 0.0)
        return x
    elseif (x <= 0.0)
        return x * 0.01
    end
end

ReLuL(20)

function softmax(x)
    return exp(x) / sum(exp(x))
end

class Softmax(Activation):
    @staticmethod
    def evaluate(z):
        z = np.exp(z)
        return z / np.sum(z, axis=0)

# Create the training function
# Step 1 - feed the data through the 3 layers whilst adjusting the inputs by link weights + squishing with sigmoid function at each output node.
# Layer 1 = inputs (inputs * weights)
# Layer 2 = hidden layer = squish with sigmoid the output of layer 1 - the output of hidden layer then is * by link weights (who) which then becomes the input to the output layer
# Layer 3 = output = squish output of hidden layer by sigmoid to obtain final output
# Step 2
# Calculate the error - the network output (what the network thinks the answer is) - the target
# Step 3 - Update the link weights in proporation to the error (ie if a link weight is higher than another link weight to same node - adjust the larger weight by the bigger fraciton)
# Output Errors = targets - final_outputs
# Hidden layer output errors = who' * output_errors
# Using chain rule - update the link weights from the hidden layer to output (who) and again for the input to the hidden layer (wih)
# This updating of weights is how the network learns!
drop_out = false
drop_out_p = .8
function train(wih::Array{Float64,2}, who::Array{Float64,2},inputs::Array{Float64,2}, targets::Array{Float64,2}; lr::Float64=0.1, drop_out::Bool=true, drop_out_p::Float64=.8) #wih, who, lr, inputs, targets)
    # Matrix multiplication Rules - https://www.mathsisfun.com/algebra/matrix-multiplying.html
    # [m * n] * [n * p] = [m * p] # the last n must match the first n
    # [row, column]
    # wih dim = [200,784]
    # input dim = [784, 1]
    hidden_inputs = (wih * inputs + bih) # hidd_input dim = [200, 1]

    # Apply sigmoid function to each element the hidden layer output
    # hidden_inputs dim = [200, 1]
    if drop_out == false
    hidden_outputs = sigmoid.(hidden_inputs)
    elseif drop_out == true
    hidden_outputs = sigmoid.(hidden_inputs)
    hidden_outputs .= (hidden_outputs .* rand(Binomial(1,drop_out_p), hnodes)) ./ drop_out_p # Apply drop out by hidden output * 0,1's / p (1/p - scaled so we dont need to apply at test time)(inverted dropout)
    end

    # Apply the weights to the hidden layer outputs which then are the final layers inputs
    # who dims = [10,200]
    # hidden_outputs dims = [200, 1]
    final_inputs = (who * hidden_outputs) + bho

    # Apply sigmoid function to each element the final layer output
    # final_inputs dim = [10,1]
    if drop_out == false
    final_outputs = sigmoid.(final_inputs)
    elseif drop_out == true
    final_outputs = sigmoid.(final_inputs)
    final_outputs .= (final_outputs .* rand(Binomial(1,drop_out_p), onodes)) ./ drop_out_p # Apply drop out by hidden output * 0,1's / p (1/p - scaled so we dont need to apply at test time)(inverted dropout)
end

    # output layer error is the (target - actual)
    # output_errors dim = [10,1]
    output_errors = targets - final_outputs
    # hidden_errors dim = [200,1]
    # hidden layer error is the output_errors, split by weights, recombined at hidden nodes
    hidden_errors = who' * output_errors

    # Update weights of the neural network - this is how the network learns by adjusting the link weights in propotation to the link weight and error.
    # update the weights for the links between the hidden and output layers
    who .= who .+ (lr .* (output_errors .* final_outputs .* (1.0 .- final_outputs)) * hidden_outputs')
    bho .= bho .+ (lr .* (output_errors .* final_outputs .* (1.0 .- final_outputs)) * hidden_outputs')

    # update the weights for the links between the input and hidden layers
    wih .= wih .+ (lr .* (hidden_errors .* hidden_outputs .* (1.0 .- hidden_outputs)) * inputs')
    bih .= bih .+ (lr .* (hidden_errors .* hidden_outputs .* (1.0 .- hidden_outputs)) * inputs')

    return(wih, who, final_outputs)
end

lr = .01
#=
# Example of link weight updating (prior to using calclus chain rule)
# Refine weight on output node (1x node example with two w11 and w21 connected links)
# Splits the error in proportion to the error ie if w11 is contributing more to error, adjust that by bigger weight.
# Fraction of e1 to refine w11
w11 = 6.0  # weight 1, 1 means node connection 1 to node 1
w21 = 3.0 # weight 2, 1 means node connection 2 to 1
refine_w11 = w11 / (w11+w21) # 2.3

# Fraction of e1 to refine w21
refine_w11 = w21 / (w21+w11) # 1/3

# Recombine split errors for the links using the error backpropogation
# ehidden,1 = sum of split errors on links w11 and w12
# Link weights beween hidden and output layer example
w11 = 2.0
w12 = 1.0
w21 = 3.0
w22 = 4.0
e1 = 0.8 # final output error
e2 = 0.5 # final output error
# Find the sum of the errors on e1 and e2
error_out_e1 = e1 * (w11 / (w11 + w21)) + e2 * (w12 / (w12 + w22)) # Sum of split errors
# Flip to find the other summed link error
error_out_e2 = e2 * (w22 / (w22 + w12)) + e1 * (w21 / (w21 + w11)) # Sum of split errors

# Back propogate error_out_e1 and error_out_e2 to the hidden layer input weights
w11 = 3.0
w12 = 1.0
w21 = 2.0
w22 = 7.0
e1 = error_out_e1 # Hidden layer output sum of proportional errors
e2 = error_out_e2 # Hidden layer output sum of proportional errors
error_out_e1 = e1 * (w11 / (w11 + w21)) + e2 * (w12 / (w12 + w22)) # Sum of split errors
error_out_e2 = e2 * (w22 / (w22 + w12)) + e1 * (w21 / (w21 + w11))

# simply in matrix
(w11 / (w11 + w21)) (w12 / (w12 + w22))     e1
(w21 / (w21 + w11)) (w22 / (w22 + w12))     e2
=#

# Function for applying the network to unseen or test data not seen (trained) by the network before
function query(wih::Array{Float64,2}, who::Array{Float64,2},inputs::Array{Float64,2})
   # calculate signals into hidden layer
   hidden_inputs = wih * inputs
   # calculate the signals emerging from hidden layer
   hidden_outputs = sigmoid.(hidden_inputs)

   # calculate signals into final output layer
   final_inputs = who * hidden_outputs
   # calculate the signals emerging from final output layer
   final_outputs = sigmoid.(final_inputs)

   return(final_outputs)
end

# Number of input, hidden and output nodes
input_nodes = size(training_set_aug_plus,2)-1 # equal to the size of the 28 * 28 image (784 elements)
hidden_nodes = 256 # 200
output_nodes = 10 # Equal to the number of labels in this case numbers 0 to 9 (10 elements total)

# learning rate
learning_rate = 0.01 # small learning rate for ReLU

# Create instance of neural network
inodes = input_nodes
hnodes = hidden_nodes
onodes = output_nodes
# As randn generates numbers from a standard normal distribution we have to multiply by the standard deviation
wih = randn(hnodes,inodes) * inodes^(-0.5)
who = randn(onodes,hnodes) * hnodes^(-0.5)
bih = zero(wih)
bho = zero(who)

# Train the neural network
# 4.22.2019 - Add functinality to save the numbers which the network did not get correct
# Epochs is the number of times the training data set is used for training
epochs = 15
# Initialize training progress
learning_curve = fill(0,size(scaled_train_inputs,1))
training_accuracy = fill(0.0,size(scaled_train_inputs,1))
# Check that the neural network is converging
@inbounds for i in 1:epochs
   # go through all records in the training data set
   @inbounds for j in 1:size(scaled_train_inputs,1)
       # Subset by row [j,1] = label [j,2:length(scaled_train_inputs)] is the scaled image data
       label = Int64.(scaled_train_inputs[j,1])
       inputs = reshape(Float64.(scaled_train_inputs[j,2:size(scaled_train_inputs,2)]),784,1) # rotr90() for changing dim to [784,1] from [1,784] to meet matrix multiplication rules
       # create the target output values (all 0.01, except the desired label which is 0.99)
       targets = zeros(10,1)
       targets .= 0.01 .+ targets
       # Set the target .99 to the correct index position within the target array
       targets[label+1] = 0.99 # We +1 because the number convention starts at 0 to 9. Thus array position 6 is actually target number 5.
       global wih, who, final_outputs = train(wih, who, inputs, targets; lr=learning_rate,drop_out=true,drop_out_p=.8)
       # Check what the network thought the number was
       network_train_output = argmax(final_outputs)[1]-1 # adjust index position to account from 0 start
       if (label == network_train_output)
           # network's answer matches correct answer, add 1 to learning curve
           learning_curve[j] = 1
       else
           # network's answer doesn't match correct answer, add 0 to learning curve
           learning_curve[j] = 0
       end
       if (j == 1)
           global learning_curve = fill(0,size(scaled_train_inputs,1))
       end
       print("This is epoch ",i,"\nThis is iteration ",j,"\nCheck hidden layer output weight links are updating ", who[1],"\n")
       print("Training Accuracy ",(sum(learning_curve) / j),"\n")
       training_accuracy[j] = (sum(learning_curve) / j)
   end
end

# Plot the learning curve
outie = training_accuracy[1:500]
plot(outie,title="Training Accuracy")

# Call the network on an indivdual entry
example_call = Float64.(scaled_train_inputs[500,2:size(scaled_train_inputs,2)])
example_call_image = rotl90(reshape(example_call,28,28))
# Plot Number
heatmap(example_call_image)
#------------------------------
using Plots; pgfplots()
test = randn(3,3)
heatmap(test, c = :greys)
plot(test)
#-------------------------------------
# Grab the known correct label
correct_label = scaled_train_inputs[500,1]
# Apply the trained network to the input example
outputs = query(wih, who, example_call)
# Check what the label thought the number was
network_label = argmax(outputs)[1]-1
print("Network thinks the number is ",network_label)
if network_label ==  correct_label
    print("The neural network correctly identified the hand written number - ", network_label)
else
    print("The neural network incorrectly identified the hand written number")
end

# Test the neural network
# Initialize the scorecard output the same size as the number of test set labels
scorecard = fill(0,size(scaled_test_inputs,1))
test_accuracy = fill(0.0,size(scaled_test_inputs,1))
kaggle_out = fill(0.0,size(scaled_test_inputs,1))
kaggle_out = DataFrame(ImageId = collect(1:size(scaled_test_inputs,1)), Label = kaggle_out)

# Input the test set data over the trained neural network to see how well it does
i=1
@inbounds for i in 1:size(scaled_test_inputs,1)
    let wih = wih, who = who
    #correct_label = Int64.(scaled_test_inputs[i,1])
    inputs = reshape(Float64.(scaled_test_inputs[i,1:size(scaled_test_inputs,2)]),784,1) # rotr90() for changing dim to [784,1] from [1,784] to meet matrix multiplication rules
    # query the network
    outputs = query(wih, who, inputs) # Test the network over the test set data (wih and who already set during training - no weight updating happens this time round)
    # Find element position of what the network thinks the output is ie if output is a .96 in element position 6 then the network thinks the label was 5.
    # This is because we have hand written numbers 0 to 9, total 10 outputs. Starting from index position 1 = label 0 hence index position 6 then = label 5
    label = argmax(outputs)[1]-1 # adjust index position to account from 0 start
    # Kaggle out
    kaggle_out[i,2] = label
    test_accuracy[i] = sum(scorecard) / i
    # append correct or incorrect to list
    #if (label == correct_label)
        # network's answer matches correct answer, add 1 to scorecard
    #    scorecard[i] = 1
    #else
        # network's answer doesn't match correct answer, add 0 to scorecard
    #    scorecard[i] = 0
    #    end
    #    print("Accuracy ",sum(scorecard) / i,"\n")
    end
end

 CSV.write("C:/Users/Andrew.Bannerman/Desktop/Julia/andrew_bannerman_submission_twentythree.csv", kaggle_out;delim=',')

# Percentage Correct
percentage = sum(scorecard) / size(scorecard,1)

# Change log
# 4.22.2019 - added a % print function to check train / test accuracy


text:
function interpolate(A::AbstractArray, ::NoInterp)
    interpolate(Int, eltype(A), A, NoInterp())
end

# How many non-NoInterp dimensions are there?
count_interp_dims(::Type{NoInterp}) = 0

interpdegree(::NoInterp) = NoInterp()

iscomplete(::NoInterp) = true

prefilter(::Type{TWeights}, ::Type{TC}, A::AbstractArray, ::NoInterp) where {TWeights, TC} = A

lbound(ax, ::NoInterp) = first(ax)
ubound(ax, ::NoInterp) = last(ax)

weightedindex_parts(fs, ::NoInterp, ax, x::Number) = Int(x)

# positions(::NoInterp, ax, x) = (Int(x),), 0

# value_weights(::NoInterp, δx) = (oneunit(δx),)
# gradient_weights(::NoInterp, δx) = (NoInterp(),)
# hessian_weights(::NoInterp, δx) = (NoInterp(),)

padded_axis(ax::AbstractUnitRange, ::NoInterp) = ax


text:
#
println("Bb) Tests of the parity non-conservation, Schiff moment and anapole moment amplitudes.")
#
grid = JAC.Radial.Grid(true)
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(26.),
                        configs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")] )

wxa  = perform(wa; output=true)
wma  = wxa["multiplet:"]

nModel = Nuclear.Model(26.);    flow = 6;    fup = 8;   ilow = 1;   iup = 3

println("\n\nWeak-charge and Schiff-moment amplitudes:\n")
for  finalLevel in wma.levels
    for  initialLevel in wma.levels
        if  flow <= finalLevel.index <= fup   &&    ilow <= initialLevel.index <= iup
            ##x println("Levels with indices f = $(finalLevel.index) and i = $(initialLevel.index)")
            ParityNonConservation.weakChargeAmplitude(finalLevel, initialLevel, nModel, grid; display=true)
            ParityNonConservation.schiffMomentAmplitude(finalLevel, initialLevel, nModel, grid; display=true)
        end
    end
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL

"""
	lar4mesh(points,faces)

Transform a LAR model (verts, cells) in a GLMesh

# Example

```
julia> using LinearAlgebraicRepresentation

julia> Lar = LinearAlgebraicRepresentation
LinearAlgebraicRepresentation

julia> V,FV = Lar.simplexGrid([2,2])
([0.0 1.0 … 1.0 2.0; 0.0 0.0 … 2.0 2.0],
Array{Int64,1}[[1, 2, 4], [2, 4, 5], [2, 3, 5], [3, 5, 6], [4, 5, 7], [5, 7, 8], [5, 6, 8], [6, 8, 9]])

julia> vertices,normals = GL.lar4mesh(GL.two2three(V'),FV);
(Float32[0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0  …  0.0, 4.0, 4.0, 0.0, 3.0, 4.0, 0.0, 4.0, 3.0, 0.0],
Float32[0.0, 0.0, -0.5, 0.0, 0.0, -0.5, 0.0, 0.0, -0.5, 0.0  …  -0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.0, 0.5])
```
"""
function lar4mesh(verts,cells) # cells are triangles
	vertices=Vector{Float32}()
	normals =Vector{Float32}()
	for cell in cells
		p3,p2,p1 = verts[cell[1],:],verts[cell[2],:],verts[cell[3],:]
		p1 = convert(Point3d, p1)
		p2 = convert(Point3d, p2)
		p3 = convert(Point3d, p3)
		n=0.5*computeNormal(p1::Point3d,p2::Point3d,p3::Point3d)

		append!(vertices,p1); append!(normals,n)
		append!(vertices,p2); append!(normals,n)
		append!(vertices,p3); append!(normals,n)
	end
	return vertices,normals
end

"""
	two2three(points)

Transform a ``n×2`` `Array{Float64,2}` of 2D points to an ``nx3`` array of 3D points, adding a column of zeros.
# Example
```
julia> V
2×9 Array{Float64,2}:
 0.0  1.0  2.0  0.0  1.0  2.0  0.0  1.0  2.0
 0.0  0.0  0.0  1.0  1.0  1.0  2.0  2.0  2.0

 julia> GL.two2three(V')
 9×3 Array{Float64,2}:
  0.0  0.0  0.0
  1.0  0.0  0.0
  2.0  0.0  0.0
  0.0  1.0  0.0
  1.0  1.0  0.0
  2.0  1.0  0.0
  0.0  2.0  0.0
  1.0  2.0  0.0
  2.0  2.0  0.0
```
"""
function two2three(points)
      return [points zeros(Float64,size(points,1),1)]
end

"""
	glGenBuffer()

Return 1 buffer object name.
"""
function glGenBuffer()
	id = GLuint[0]
	glGenBuffers(1, id)
	glCheckError("generating a buffer, array, or texture")
	id[]
end

"""
	glGenVertexArray()

Specifies the number (one) of vertex array object names to generate.
"""
function glGenVertexArray()
	id = GLuint[0]
	glGenVertexArrays(1, id)
	glCheckError("generating a buffer, array, or texture")
	id[]
end

"""
	glLineWidth()

Define the width of the line primitive
"""
function glLineWidth()
	width = convert(UInt32, 2)
	ModernGL.glLineWidth(width)
end

"""
	glPointSize()

Define the size of the point primitive
"""
function glPointSize()
	size = convert(UInt32, 8) # multiples of 2 give sprites!!
	ModernGL.glPointSize(size)
end

"""
	glCheckError(actionName="")

`glGetError` returns the value of the error flag. Each detectable error is assigned a numeric code and symbolic name. When an error occurs, the error flag is set to the appropriate error code value.
"""
function glCheckError(actionName="")
	message = glErrorMessage()
	if length(message) > 0
		if length(actionName) > 0
		error("Error ", actionName, ": ", message)
		else
		error("Error: ", message)
		end
	end
end

"""
	glErrorMessage()

Error message management.
"""
function glErrorMessage()
	err = glGetError()
	if err == GL_NO_ERROR return "" end
	if err == GL_INVALID_ENUM return "GL_INVALID_ENUM" end
	if err == GL_INVALID_VALUE return "GL_INVALID_VALUE"  end
	if err == GL_INVALID_OPERATION return "GL_INVALID_OPERATION"  end
	if err == GL_INVALID_FRAMEBUFFER_OPERATION return "GL_INVALID_FRAMEBUFFER_OPERATION"  end
	if err == GL_OUT_OF_MEMORY return "GL_OUT_OF_MEMORY" end
	return "Unknown OpenGL error with error code"
end

const __release_gpu_resources__=[]

"""
	glDeleteLater(fun::Function)
"""
function glDeleteLater(fun::Function)
	global __release_gpu_resources__
	append!(__release_gpu_resources__,[fun])
end

"""
	glDeleteNow()
"""
function glDeleteNow()
	global __release_gpu_resources__
	for fun in __release_gpu_resources__
		fun()
	end
end

"""
	normalize2(V::Lar.Points; flag=true::Bool)::Lar.Points

2D normalization transformation (isomorphic by defaults) of model
vertices to normalized coordinates ``[0,1]^2``. Used with SVG importing.
"""
function normalize2(V::Lar.Points; flag=true, fold=-1)
	m,n = size(V)
	if m > n # V by rows
		V = convert(Lar.Points, V')
	end

	xmin = minimum(V[1,:]); ymin = minimum(V[2,:]);
	xmax = maximum(V[1,:]); ymax = maximum(V[2,:]);
	box = [[xmin; ymin] [xmax; ymax]]	# containment box
	aspectratio = (xmax-xmin)/(ymax-ymin)
	if flag
		if aspectratio > 1
			umin = 0; umax = 1
			vmin = 0; vmax = 1/aspectratio; ty = vmax
		elseif aspectratio < 1
			umin = 0; umax = aspectratio
			vmin = 0; vmax = 1; ty = vmax
		end
		T = fold==-1 ? Lar.t(0,ty) : Lar.t(0,0) *
		 	Lar.s(1,fold) * Lar.s((umax-umin), (vmax-vmin)) *
			Lar.s(1/(xmax-xmin),1/(ymax-ymin)) * Lar.t(-xmin,-ymin)
	else
		T = Lar.t(0, ymax-ymin) * Lar.s(1,-1)
	end
	dim = size(V,1)
	W = T[1:dim,:] * [V;ones(1,size(V,2))]
	#V = map( x->round(x,digits=8), W )
	V = map(Lar.approxVal(8), W)

	if m > n # V by rows
		V = convert(Lar.Points, V')
	end
	return V
end

function normalize3(V::Lar.Points, flag=true)
	Vxy = normalize2(V[1:2,:], flag=flag, fold=1)
	if size(V,1) == 3
		zmin = minimum(V[3,:]); zmax = maximum(V[3,:]);
		Vz = (V[3,:] .- zmin) / (zmax - zmin)
	elseif size(V,1) == 2
		Vz = zeros(1,size(V,2))
	end
	return [Vxy; Vz]
end

function explodecells(V,FVs,sx=1.2,sy=1.2,sz=1.2)
	V,FVs
	outVerts, outCells = [],[]
	grids = []
	for FV in FVs
		vertidx = sort(collect(Set(cat(FV))))
		vcell = V[:,vertidx]
		vdict = Dict(zip(vertidx,1:length(vertidx)))

		center = sum(vcell,dims=2)/size(vcell,2)
		scaled_center = size(center,1)==2 ? center .* [sx;sy] : center .* [sx;sy;sz]
		translation_vector = scaled_center - center
		cellverts = vcell .+ translation_vector
		newcells = [[vdict[v] for v in cell] for cell in FV]

		mesh = cellverts, newcells
		push!(grids, mesh)
	end
	return grids
end


text:
# initialize hidden and cell arrays
function initstate(atype, hidden, batchsize)
    state = Array(Any, 2)
    state[1] = zeros(batchsize, hidden)
    state[2] = zeros(batchsize, hidden)
    return map(s->convert(atype,s), state)
end

# initialize all weights of decoder network
function initweights(o::Dict)
    w = Dict()
    L,D = o[:visual]
    w["wdec"] = o[:winit]*randn(o[:embed]+o[:hidden]+D, 4*o[:hidden])
    w["bdec"] = zeros(1, 4*o[:hidden])
    w["bdec"][1:o[:hidden]] = 1 # forget gate bias
    w["wsoft"] = o[:winit]*randn(o[:hidden], o[:vocabsize])
    w["bsoft"] = zeros(1, o[:vocabsize])
    w["wemb"] = o[:winit]*randn(o[:vocabsize], o[:embed])
    w["wce"] = o[:winit]*randn(D, D)
    w["whe"] = o[:winit]*randn(o[:hidden], D)
    w["watt"] = o[:winit]*randn(D,1)
    return convert_weight(o[:atype], w)
end

max_stars_repo_name: ilkerkesen/show-attend-and-tell

max_issues_repo_name: ilkerkesen/show-attend-and-tell

max_forks_repo_name: ilkerkesen/show-attend-and-tell


text:
function gmsh_mesh(filename)

    if !isfile(filename)
        println("Error in gmsh.read():")
        println("  Could not find " * filename * "!")
        exit(1)
    end

    gmsh = open(filename)

    line = ""

    find_keyword(name) = begin
        while line != "\$"*name
            line = readline(gmsh)
        end
    end

    find_keyword("MeshFormat")
    mesh_format = readline(gmsh)

    physical_lines=Dict{Int,String}()
    find_keyword("PhysicalNames")
    nphys = parse(Int, readline(gmsh) )
    for i=1:nphys
        words = split(readline(gmsh))
        d   = parse(Int, words[1])
        num = parse(Int, words[2])
        name = words[3][2:end-1]
        if d==1
            physical_lines[num] = name
        end
    end

    find_keyword("Nodes")
    nnodes = parse(Int, readline(gmsh))
    nodes = Array{Vector,1}()
    for i=1:nnodes
        words = split(readline(gmsh))
        x = parse(Float64, words[2])
        y = parse(Float64, words[3])
        push!(nodes, Vector(x,y))
    end

    Edge = Tuple{Int64,Int64}

    find_keyword("Elements")
    cell2point = LabelListList()
    boundary_lines = Dict{Int, Array{Edge,1}}()
    for (k,v) in physical_lines
        boundary_lines[k] = []
    end
    n = parse(Int, readline(gmsh))
    for i = 1:n
        numbers = [parse(Int,w) for w in split(readline(gmsh))]
        eid   = numbers[1]
        etype = numbers[2]
        ntags = numbers[3]

        if etype==1  # line
            push!( boundary_lines[numbers[4]], (numbers[4+ntags],numbers[5+ntags]) )
        elseif etype==2 || etype==3  # triangle or quad
            push!(cell2point, numbers[4+ntags:end])
        else
            error("Unknown element type ", etype)
        end
    end

    close(gmsh)

    #### Build data structures for Mesh

    edge_owner = Dict{Edge,Label}()

    owner = LabelList()
    neighbor = LabelList()
    face2point = LabelListList()

    push_face(p1, p2, o, n=nothing) = begin
        push!(owner, o)
        if n != nothing; push!(neighbor, n); end
        push!(face2point, [p1, p2])
    end

    # Process all faces and push internal faces into owner/neighbor/face2point
    for o = 1:length(cell2point)
        pts = [cell2point[o]; cell2point[o][1]]
        for i = 1:length(pts)-1
            p1, p2 = pts[i], pts[i+1]

            if haskey(edge_owner, (p2,p1))
                n = edge_owner[(p2,p1)]
                delete!(edge_owner, (p2,p1))
                push_face(p1, p2, o, n)
            else
                edge_owner[(p1,p2)] = o
            end
        end
    end

    patch = PatchInfoList()
    for (tag, edges) in boundary_lines
        name = physical_lines[tag]
        start = length(owner)+1
        stop  = start + length(edges) - 1
        push!(patch, PatchInfo(name, start:stop))

        for e in edges
            p1, p2 = e[1], e[2]
            if haskey(edge_owner, (p1, p2))
                o = edge_owner[(p1, p2)]
                push_face(p1, p2, o)
            elseif haskey(edge_owner, (p2, p1))
                o = edge_owner[(p2, p1)]
                push_face(p2, p1, o)
            else
                error("Error in processing mesh faces!")
            end
        end
    end

    return Mesh(nodes, owner, neighbor, patch, cell2point, face2point)
end

#mesh = gmsh_mesh("ctverec.msh");
#println(mesh)

max_stars_repo_name: furstj/NSP

max_issues_repo_name: furstj/NSP

max_forks_repo_name: furstj/NSP


text:
module BayesianSparseFactorGmatrix

	using Statistics
	using Random
	using Distributions
	using LinearAlgebra
	using BlockDiagonals
	using Kronecker
	using ProgressMeter

	include("utils.jl")
	include("types.jl")
	include("sampling_functions.jl")
	include("posterior_utils.jl")
	include("BSFG.jl")

	## Initialize run parameters & prior hyper-parameters
	# run parameters
	# _burn = 1000;
	# _sp = 1000;
	# _thin = 10

	# _b0 = 1.;
	# _b1 = 0.0005

	# _epsilon = 1e-2
	# _h2_divisions = 100

	# # prior hyperparamters

	# _k_init = 20; # initial number of factors to initialize
	# _as = 2.; _bs = 1. / 10; # inverse gamma hyperparameters for model residuals; as well as non-factor random effects
	# _df = 3; # degrees of freedom for t-distribution of ARD prior on factor loadings
	# _ad1 = 2.1; _bd1 = 1. /20; # inverse gamma hyperparamters for first factor shrinkage multiplier [/delta_1]
	# _ad2 = 3.; _bd2 = 1.; # inverse gamma hyperparamters for remaining factor shrinkage multiplier [/delta_i, i \in 2...k]

	# _k_min = 1e-1; # minimum factor loading size to report in running status
	# _prop = 1.00;  # proportion of redundant elements within columns necessary to drop column
	function runBSFGModel(Y, X, A, Z_1; burn::Int64=2000, sp::Int64=1000, thin::Int64=20,
	                      b0::Float64=1., b1::Float64=0.0005, epsilon::Float64=1e-2,
	                      h2_divisions::Int64=100, k_init::Int64=20,
	                      as::Float64=2., bs::Float64=0.1, df::Float64=3.,
	                      ad1::Float64=2.1, bd1::Float64=1. / 20,
	                      ad2::Float64=3.0, bd2::Float64=1.,
	                      k_min::Float64=0.1, prop::Float64=1.)
	    Pr = Priors(burn, sp, thin, b0, b1, epsilon, h2_divisions, k_init,
	                as, bs, df, ad1, bd1, ad2, bd2, k_min, prop)
	    D = InputData(Y, X, A, Z_1);
	    resid = Residuals(Pr.as, Pr.bs, D.Y, D.p, rand(Gamma(Pr.as, 1. / Pr.bs), D.p))
	    Factors = LatentFactors(D.r1, D.n, D.p, Pr.k_init, Pr.df, Pr.ad1, Pr.bd1, Pr.ad2, Pr.bd2,
	                            Pr.h2_divisions, Pr.sp, Pr.nrun)
	    genetic_effects = GeneticEffects(D.n, Pr.as, Pr.bs, D.r1, D.p, Factors.k, Factors.h2)
	    interaction_effects = InteractionEffects(Pr.as, Pr.bs, D.p, D.r2)
	    fixed_effects = FixedEffects(D.b, D.p)
	    Factors.scores = genetic_effects.U * D.Z_1 + broadcast(*, randn(Factors.k, D.n), sqrt.(1. .- Factors.h2))
	    Posterior = PosteriorSample(D.n, D.p, Pr.k_init, Pr.sp, D.r1, D.r2, D.b)

	    sampleBSF_G!(Posterior::PosteriorSample, genetic_effects::GeneticEffects, Factors::LatentFactors,
	             resid::Residuals, interaction_effects::InteractionEffects, fixed_effects::FixedEffects,
	             D::InputData, Pr::Priors)

	    posterior_mean = PosteriorMeans(Posterior, D, Pr)

	    [posterior_mean, Posterior, D, Pr, fixed_effects]
	end

	export runBSFGModel, cov2cor;

end


text:
using FEMAssembler

# scalar parameters
E  = 1.0 # Young modulus in Pa
ν  = .3  # Poisson's ratio
p  = 2   # Tension load in Pa
Lx = 2   # Dimension in x of the boxD in m
Ly = 1   # Dimension in y of the boxD in m
Lz = 1   # Dimension in z of the boxD in m

# Lamé parameters
λ = E*ν / ( (1+ν) * (1-2*ν) )
μ = E / ( 2 * (1+ν) )

# set material structs
SVKmat = materialsData( E, ν );

# set element structs
triangl = elementsData( "face", 0 )
tetra   = elementsData( "tetrahedron" , 0 )

# Set the loads
loadsTimeFact= t-> p*t

# boundary conditions
BCs = [ boundaryCondsData( [ ], [ ], [1], [ loadsTimeFact( 0 ) ] ),  # tension
        boundaryCondsData( [1], [0], [ ], [ ] ),  # non-friction constraint x=0
        boundaryCondsData( [3], [0], [ ], [ ] ),  # non-friction constraint y=0
        boundaryCondsData( [5], [0], [ ], [ ] ) ] # non-friction constraint z=0

# the connectivity matrix is given by a vector of vectors.
elemNodalConnec =[ [ 5 8 6   ] ,  # loaded face
                   [ 6 8 7   ] ,  # loaded face
                   [ 4 1 2   ] ,  # x=0 supp face
                   [ 4 2 3   ] ,  # x=0 supp face
                   [ 6 2 1   ] ,  # y=0 supp face
                   [ 6 1 5   ] ,  # y=0 supp face
                   [ 1 4 5   ] ,  # z=0 supp face
                   [ 4 8 5   ] ,  # z=0 supp face
                   [ 1 4 2 6 ] ,  # tetrahedron
                   [ 6 2 3 4 ] ,  # tetrahedron
                   [ 4 3 6 7 ] ,  # tetrahedron
                   [ 4 1 5 6 ] ,  # tetrahedron
                   [ 4 6 5 8 ] ,  # tetrahedron
                   [ 4 7 6 8 ]  ] # tetrahedron

# matrix with MEBI indexes of each property (on each row)
# MEBI = Materials / Elements / BoundaryConditions / InitialConditions
MEBIValsMat = [ 0 1 1 0 ;  # BC 1
                0 1 2 0 ;  # BC 2
                0 1 3 0 ;  # BC 3
                0 1 4 0 ;  # BC 4
                1 2 0 0 ]

MEBIVec = Int.( vcat( ones(2), 2*ones(2), 3*ones(2), 4*ones(2), 5*ones(6) ) )

# The coordinate matrix is given by
nodalCoords = [ 0    0    0 ;
                0    0   Lz ;
                0   Ly   Lz ;
                0   Ly    0 ;
                Lx   0    0 ;
                Lx   0   Lz ;
                Lx  Ly   Lz ;
                Lx  Ly    0 ]

                finalTime = 1;
                deltaTime = .5 ;

print("Starting time loop.\n")
currTime = 0.0
# while currTime < finalTime
#     global currTime = currTime + deltaTime
#     print("Solving time: ", currTime,"\n")

#     BCs[1].NeumannNodalVals = [loadsTimeFact( currTime )]
#     print( "  current pressure ", BCs[1].NeumannNodalVals, "\n" )

#     # Run FEMAssembler
#     KG, FG, Kmatrices = assembler( nodalCoords, elemNodalConnec, MEBIValsMat, MEBIVec, [SVKmat], [ triangl, tetra ], BCs, [] )

# end


text:
# Converts a MAT file to an HDF5 file
function mat_to_h5(fn, out)
    matdata = matread(fn)
    h5open(out, "w") do file
        for (k, v) in matdata
            write(file, k, size(v, 2) == 1 ? vec(v) : v)
        end
    end
end


text:
using Knet.KnetArrays: KnetPtr, KnetArray, Cptr
using AutoGrad: Param
using CUDA: CUDA, functional, CuPtr, CuArray

const JLDMODE=Val(0)
const GPUMODE=Val(1)
const CPUMODE=Val(2)

# Do not use type asserts because type may change
jld2serialize(x) = _ser(x,IdDict(),JLDMODE) # Regular JLD2 functions now work, this is no longer needed
gpucopy(x)   = _ser(x,IdDict(),GPUMODE)
cpucopy(x)   = _ser(x,IdDict(),CPUMODE)

function _ser(x::KnetPtr,s::IdDict,::typeof(JLDMODE))
    if !haskey(s,x)
        if isa(x.ptr, Cptr) && (x.dev >= 0)
            a = Array{UInt8}(undef,x.len)
            unsafe_copyto!(pointer(a), CuPtr{UInt8}(UInt(x.ptr)), x.len)
            s[x] = KnetPtr(a,x.len,-1,nothing)
        elseif isa(x.ptr, Array{UInt8,1})
            if CUDA.functional()
                s[x] = KnetPtr(x.len)
                unsafe_copyto!(CuPtr{UInt8}(UInt(s[x].ptr)), pointer(x.ptr), x.len)
            else
                s[x] = x  # Leave conversion to array to KnetArray
            end
        else
            error("Unrecognized KnetPtr")
        end
    end
    return s[x]
end

function _ser(x::KnetArray{T,N},s::IdDict,m::typeof(JLDMODE)) where {T,N}
    if !haskey(s,x)
        if isa(x.ptr.ptr, Array) && !CUDA.functional()
            s[x] = copy(reshape(reinterpret(eltype(x),view(x.ptr.ptr,1:sizeof(T)*length(x))),size(x)))
        else
            s[x] = KnetArray{T,N}(_ser(x.ptr,s,m),x.dims)
        end
    end
    return s[x]
end

_ser(x::KnetArray,s::IdDict,::typeof(GPUMODE))=x
_ser(x::KnetArray,s::IdDict,::typeof(CPUMODE))=(haskey(s,x) ? s[x] : s[x]=Array(x))

_ser(x::CuArray,s::IdDict,::typeof(GPUMODE))=x
_ser(x::CuArray,s::IdDict,::typeof(CPUMODE))=(haskey(s,x) ? s[x] : s[x]=Array(x))

# Partially fixes the issue: when KA converts to A because no gpu, surrounding parametric types remain Param{KA}.
# However other container types that include KnetArray may still have an inconsistent parametric type problem.
_ser(x::Param, s::IdDict, m::Val)=(haskey(s,x) ? s[x] : s[x]=Param(_ser(x.value,s,m),_ser(x.opt,s,m)))

_ser(x::Array, s::IdDict, m::Val) = (haskey(s, x) ? s[x] : s[x] = _ser_array_t(x, eltype(x), s, m))

function _ser_array_t(@nospecialize(x), T, s::IdDict, m::Val)
    if !isbitstype(T)
        # map(xi->_ser(xi,s,m), x) # this fails with unassigned values
        dest = similar(x,Any)   # convert eltype to Any because it may change
        # stackdict[x] = dest # we do this in the caller
        for i = 1:(length(x)::Int)
            if ccall(:jl_array_isassigned, Cint, (Any, Csize_t), x, i-1) != 0
                xi = ccall(:jl_arrayref, Any, (Any, Csize_t), x, i-1)
                if !isbits(xi)
                    xi = _ser(xi, s, m) # deepcopy_internal(xi, stackdict)::typeof(xi)
                end
                ccall(:jl_arrayset, Cvoid, (Any, Any, Csize_t), dest, xi, i-1)
            end
        end
        return dest
    elseif m === GPUMODE
        KnetArray(x)
    else
        x  # Note: deepcopy actually copies here, we don't
    end
end

# Generic serialization rules from deepcopy.jl in julia/base:
_ser(x::Union{Symbol,Core.MethodInstance,Method,GlobalRef,DataType,Union,UnionAll,Task,Regex},::IdDict,::Val) = x
_ser(x::Tuple, s::IdDict, m::Val) = ntuple(i->_ser(x[i], s, m), length(x))
_ser(x::Module, ::IdDict, ::Val) = error("serialize of Modules not supported")
_ser(x::Core.SimpleVector, s::IdDict,m::Val) = (haskey(s, x) ? s[x] : s[x] = Core.svec(Any[_ser(x[i], s, m) for i = 1:length(x)]...))
_ser(x::String, s::IdDict,::Val) = (haskey(s, x) ? s[x] : s[x] = (GC.@preserve x unsafe_string(pointer(x), sizeof(x))))

function _ser(@nospecialize(x), stackdict::IdDict, m::Val)
    T = typeof(x)::DataType
    nf = nfields(x)
    if ismutable(x)
        if haskey(stackdict, x)
            return stackdict[x]
        end
        y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
        stackdict[x] = y
        for i in 1:nf
            if isdefined(x, i)
                xi = getfield(x, i)
                xi = _ser(xi, stackdict, m)  # removed ::typeof(xi)
                ccall(:jl_set_nth_field, Cvoid, (Any, Csize_t, Any), y, i-1, xi)
            end
        end
    elseif nf == 0 || isbitstype(T)
        y = x
    else
        flds = Vector{Any}(undef, nf)
        for i in 1:nf
            if isdefined(x, i)
                xi = getfield(x, i)
                xi = _ser(xi, stackdict, m)  # removed ::typeof(xi)
                flds[i] = xi
            else
                nf = i - 1 # rest of tail must be undefined values
                break
            end
        end
        y = ccall(:jl_new_structv, Any, (Any, Ptr{Any}, UInt32), T, flds, nf)
    end
    return y
end

function _ser(x::Union{Dict{K,V},IdDict{K,V}}, s::IdDict, m::Val) where {K,V}
    if haskey(s, x)
        return s[x] # removed ::typeof(x)
    end
    if isbitstype(K) && isbitstype(V)
        return (s[x] = x)
    end
    # Type of key or value may change unless isbitstype
    if isbitstype(K)
        dest = empty(x, K, Any)
    elseif isbitstype(V)
        dest = empty(x, Any, V)
    else
        dest = empty(x, Any, Any)
    end
    s[x] = dest
    for (k, v) in x
        dest[_ser(k, s, m)] = _ser(v, s, m)
    end
    dest
end