text:
using SparseArrays
using TimerOutputs

"""
    lift_and_project

Compute lift-and-project cuts.

Returns a list of violated split cuts and K* cuts.
"""
function lift_and_project(
    x_, sf::StandardProblem, cgcp_optimizer;
    nb_rounds_max::Int=typemax(Int),
    normalization::Symbol=:Conic,
    kcut_pre_check::Bool=true,
    kcut_post_check::Bool=true,
    strengthen_flag::Bool=true,
    ϵ_integer::Float64=1e-4,
    ϵ_cut_violation::Float64=1e-4,
    timer=TimerOutput(),
    xref::Vector{Float64}=Float64[]
)
    Kcuts = Tuple{Vector{Int}, Vector{Float64}, Float64}[]  # K*    cuts
    Scuts = Tuple{Vector{Int}, Vector{Float64}, Float64}[]  # Split cuts

    # Problem size
    m, n = size(sf.A)

    # Current lower bound
    z_ = dot(x_, sf.c)

    # Clean ̄x
    @timeit timer "Clean x_" for j in 1:n
        if abs(x_[j]) <= 1e-7
            x_[j] = 0.0
        end
    end

    # Compute fractional parts
    f = min.(ceil.(x_) .- x_, x_ .- floor.(x_))
    f .*= sf.vartypes
    p = sortperm(f, rev=true)
    @debug extrema(f)

    # Check if ̄x is conic feasible
    @timeit timer "Conic infeas." H = extract_conic_infeasibilities(x_, sf.cones) ./ 2
    η = minimum(H)

    # Cut generation loop
    nkcuts_pre = 0
    nkcuts_post = 0
    nscuts = 0
    @timeit timer "Cut-generation" for j in p

        # TODO: check stopping criterion

        # Skip non-fractional variables
        f[j] >= ϵ_integer || continue

        kcut_detected = false
        # Early check for K* cut
        @timeit timer "K*-cut pre-check" if kcut_pre_check && normalization == :Standard
            η_ = min(-f[j], f[j] - 1) / 2

            kcut_detected = (η <= -f[j] / 2) && (η <= (f[j]-1)/2)
            kcut_detected && @debug("K*-cut detected", η)

            if kcut_detected
                for (ηi, (kidx, k)) in zip(H, sf.cones)
                    if ηi <= η
                        # Polyhedral cones are
                        @debug "K* cut found" x_[kidx] k kidx ηi
                        isa(k, POLYHEDRAL_CONE) && continue

                        @assert isa(k, MOI.SecondOrderCone) "Only SOC are supported (is $k)"

                        # Most violated K* cut according to normalization
                        λ = copy(x_[kidx])
                        λ[1] = norm(λ[2:end])
                        λ[2:end] .*= -1
                        λ ./= λ[1]

                        # Record K* cut
                        push!(Kcuts, (kidx, λ, 0.0))
                    end
                end
                @timeit timer "K*-pre" nkcuts_pre += 1
            end
        end

        kcut_detected && continue

        # Look for a split cut
        pi = zeros(n)
        pi[j] = 1.0
        pi0 = floor(x_[j])

        # Build CGCP
        @timeit timer "CGCP-build" cgcp = build_cgcp(
            cgcp_optimizer,
            m, n, sf.A, sf.b, sf.cones,
            x_, pi, pi0,
            nrm=normalization,
            bridge_type=Float64,
            xref=xref
        )
        cgcp_moi = backend(cgcp)

        # Solve CGCP
        @timeit timer "CGCP-solve" MOI.optimize!(cgcp_moi)

        # TODO: track number of solves && IPM iterations

        # Check termination status
        st = MOI.get(cgcp_moi, MOI.TerminationStatus())
        pst = MOI.get(cgcp_moi, MOI.PrimalStatus())
        if !(
            (pst == MOI.FEASIBLE_POINT)
            || (st == MOI.DUAL_INFEASIBLE && pst == MOI.INFEASIBILITY_CERTIFICATE)
        )
            @warn "CGCP exited with status" st pst
            continue
        end

        # Check if found violated cut
        δ = objective_value(cgcp)
        @debug "Cut violation: $δ"
        if δ >= -ϵ_cut_violation && st != MOI.DUAL_INFEASIBLE
            # Cut is not violated
            @debug "No violated cut" δ st pst
            continue
        end

        # Get Farkas multipliers
        @timeit timer "Get Farkas" begin
            v = value.(cgcp[:v])
            norm(v, Inf) >= 1e2 && @warn "|v| is large: $(extrema(v))"
            u0 = value(cgcp[:u0])
            v0 = value(cgcp[:v0])
            λ = value.(cgcp[:λ])
            μ = value.(cgcp[:μ])
        end

        # Check that u0, v0 are not "too negative", otherwise reject the cut
        (u0 >= -1e-5 && v0 >= -1e-5) || continue

        # Clean u0, v0
        u0 = max(0.0, u0)
        v0 = max(0.0, v0)
        (abs(u0) <= 1e-6) && (u0 = 0.0)
        (abs(v0) <= 1e-6) && (v0 = 0.0)

        # Clean conic multipliers
        @timeit timer "Clean λ" for (kidx, k) in sf.cones
            kd = MOI.dual_set(k)
            if isa(kd, MOI.Nonnegatives)
                λ[kidx[1]] = max(0.0, λ[kidx[1]])
                μ[kidx[1]] = max(0.0, μ[kidx[1]])

            elseif isa(kd, MOI.Nonpositives)
                λ[kidx[1]] = min(0.0, λ[kidx[1]])
                μ[kidx[1]] = min(0.0, μ[kidx[1]])

            elseif isa(kd, MOI.Reals)
                # Nothing to do
                # α[kidx[1]] = 0.0

            elseif isa(kd, MOI.Zeros)
                # TODO
                # abs(a1[kidx[1]] - a2[kidx[1]]) <= 1e-5 || @warn "|a1 - a2| = $(abs(a1[kidx[1]] - a2[kidx[1]])) but should be 0"
                λ[kidx[1]] = 0.0
                μ[kidx[1]] = 0.0

            elseif isa(kd, MOI.SecondOrderCone)
                # Shift dual multipliers if needed
                λ[kidx[1]] = max(λ[kidx[1]], norm(λ[kidx[2:end]]))
                μ[kidx[1]] = max(μ[kidx[1]], norm(μ[kidx[2:end]]))

                # TODO: clean multipliers

            elseif isa(kd, MOI.RotatedSecondOrderCone)
                # start by shifting initial coordinates
                λ[kidx[1]] = max(0.0, λ[kidx[1]])
                λ[kidx[2]] = max(0.0, λ[kidx[2]])
                # Check norm
                η = 2 * λ[kidx[1]] * λ[kidx[2]] - sum(λ[kidx[3:end]] .^ 2)
                if η < 0.0
                    # TODO
                    # Increase first two coordinates
                end
                @assert 2 * λ[kidx[1]] * λ[kidx[2]] - sum(λ[kidx[3:end]] .^ 2) >= -1e-7
                @assert 2 * μ[kidx[1]] * μ[kidx[2]] - sum(μ[kidx[3:end]] .^ 2) >= -1e-7

            else
                @warn "Cone $(typeof(k)) not recognized in feasibility check"
            end
        end

        α = λ .- (u0 .* pi)
        β = min(-u0 * pi0, dot(sf.b, v) + v0 * (pi0 + 1))
        δ = dot(α, x_) - β

        # Check for K* cut
        kcut_flag = (iszero(u0) && β <= 1e-6) || (iszero(v0) && (β - dot(sf.b, v)) <= 1e-6)
        @timeit timer "K*-cut post-check" if kcut_post_check && kcut_flag
            @debug "K*-cut" u0 v0 η kcut_detected δ f[j] j
            # Dis-aggregate the cut
            for (kidx, k) in sf.cones
                if iszero(u0)
                    λ_ = λ[kidx]
                else
                    λ_ = μ[kidx]
                end

                # Skip small components
                norm(λ_) >= 1e-5 || continue

                # Record K* cut
                # TODO: add at most one K* cut per cone
                push!(Kcuts, (kidx, λ_, 0.0))
            end
            @timeit timer "count" nkcuts_post += 1
            continue
        end

        # TODO: strengthening
        @timeit timer "Strengthen" if strengthen_flag && (u0+v0) >= 1e-4

            Atv = sf.A'v
            for (kidx, k) in sf.cones

                isa(k, NONLINEAR_CONE) && all(sf.vartypes_implied[kidx] .*  iszero.(x_[kidx])) >= 1 && @info "Could have strengthened a cone"

                isa(k, Union{MOI.Nonpositives, MOI.Nonnegatives}) || continue
                j_ = kidx[1]

                # Only strengthen integer variables that are at zero
                (sf.vartypes_implied[j_] && iszero(x_[j_])) || continue

                α[j_] = strengthen(Atv[j_], u0, v0, k)
            end
        end

        # TODO: re-scale the cut

        # TODO: book-keeping

        # TODO: validity check (make it optional)

        # Record the cut
        for j_ in 1:n
            if abs(α[j_]) <= 1e-7
                α[j_] = 0.0
            end
        end
        abs(β) <= 1e-8 && (β = 0.0)

        # Record split cut
        α_ = sparsevec(α)
        push!(Scuts, (α_.nzind, α_.nzval, β))

        nscuts += 1
    end

    return Kcuts, Scuts
end

function extract_conic_infeasibilities(x, cones)
    # TODO: also return corresponding K* cuts
    ncones = length(cones)
    H = zeros(ncones)

    for (i, (kidx, k)) in enumerate(cones)
        isa(k, MOI.Zeros) && continue  # Skip zero cones
        H[i] = epigraph_violation(x[kidx], k)
    end

    return H
end

function epigraph_violation(x, k::MOI.Nonnegatives)
    length(x) == MOI.dimension(k) || throw(DimensionMismatch(
        "x has length $(length(x)) but cone has dimension $(MOI.dimension(k))"
    ))

    return minimum(x)
end

function epigraph_violation(x, k::MOI.Nonpositives)
    length(x) == MOI.dimension(k) || throw(DimensionMismatch(
        "x has length $(length(x)) but cone has dimension $(MOI.dimension(k))"
    ))

    return -maximum(x)
end

function epigraph_violation(x, k::MOI.Reals)
    length(x) == MOI.dimension(k) || throw(DimensionMismatch(
        "x has length $(length(x)) but cone has dimension $(MOI.dimension(k))"
    ))

    return zero(eltype(x))
end

function epigraph_violation(x, k::MOI.Zeros)
    length(x) == MOI.dimension(k) || throw(DimensionMismatch(
        "x has length $(length(x)) but cone has dimension $(MOI.dimension(k))"
    ))
    error("Cannot compute epigraph violation: zero cone has empty interior")
end

function epigraph_violation(x, k::MOI.SecondOrderCone)
    length(x) == MOI.dimension(k) || throw(DimensionMismatch(
        "x has length $(length(x)) but cone has dimension $(MOI.dimension(k))"
    ))

    return x[1] - norm(x[2:end], 2)
end

max_stars_repo_name: mtanneau/CLaP

max_issues_repo_name: mtanneau/CLaP

max_forks_repo_name: mtanneau/CLaP


text:
using Test

# H_ref = H[1:27:end]
tests_1d = Dict("diffusion_1D_expl.jl"=> [1.687700671616363e-11, 2.026554839900115e-7, 0.00026304660891236103,
                                          0.03690782595065184, 0.5597777471603071,
                                          0.9177496036501366, 0.16264630103928718, 0.003115843137683058,
                                          6.452357026295877e-6, 1.4443505264260712e-9],
                "diffusion_1D_impl.jl"=> [1.687700671616363e-11, 2.026554839900115e-7, 0.00026304660891236103,
                                          0.03690782595065184, 0.5597777471603071, 0.9177496036501366,
                                          0.16264630103928718, 0.003115843137683058, 6.452357026295877e-6,
                                          1.4443505264260712e-9],
                "diffusion_1D_damp.jl"=>[1.687700671616363e-11, 2.026554839900115e-7, 0.00026304660891236103,
                                         0.03690782595065184, 0.5597777471603071, 0.9177496036501366,
                                         0.16264630103928718, 0.003115843137683058, 6.452357026295877e-6,
                                         1.4443505264260712e-9]
                )

# H_ref = repr(H[1:33:end, 30:33:end].data)
tests_2d = Dict("iceflow.jl"=> [0.0 0.0 0.0 0.0 0.0;
                                2107.470713773525 2300.5901089255344 1839.9641691700776 1968.1291004159364 668.7717653362265;
                                1602.810338109451 2392.6458316876924 2480.2918550322306 0.0 0.0],
                 "iceflow_xpu.jl"=> [0.0 0.0 0.0 0.0 0.0;
                                     2107.470713773525 2300.5901089255344 1839.9641691700776 1968.1291004159364 668.7717653362265;
                                     1602.810338109451 2392.6458316876924 2480.2918550322306 0.0 0.0],
                 "iceflow_xpu_evo.jl"=> [0.0 0.0 0.0 0.0 0.0;
                                         0.0 0.0 0.0 0.0 0.0;
                                         0.0 251.15273565047443 0.0 0.0 0.0]
                 )

@testset "helpers.jl" begin
    include(joinpath(@__DIR__, "../scripts/helpers.jl"))
    # TODO: maybe add actual tests...
end

## Test the scripts
# note: this cannot be run in a test-set as the include acts in global scope
# note: all XPU tests relying on ParallelStencil.jl must have same dims and backend to succeed (else requires a reset)

do_visu, do_save = false, false # disable plotting & saving

for (fl, H_ref) in tests_1d
    lx = 10.0        # domain size
    nx = 256         # numerical grid resolution
    dx = lx/nx       # grid size

    println("Runninng $fl")
    xc, H = include(joinpath(@__DIR__, "../scripts", fl))
    @test xc == LinRange(dx/2, lx-dx/2, nx)
    @test H[1:27:end] ≈ H_ref
end

ref_data = load_data()
for (fl, H_ref) in tests_2d
    println("Runninng $fl")
    include(joinpath(@__DIR__, "../scripts", fl))
    @test size(ref_data[1]) == size(H)
    @test dims(ref_data[1]) == dims(H)
    @test H_ref ≈ H[1:33:end, 30:33:end].data
end

ext: jl


text:
using Yao
using Yao.Intrinsics

"""
    binary_basis(geometry...) -> BitArray

Return bases in reshaped binary form (e.g. binary images)
"""
binary_basis(geometry...) = reshape(bitarray(geometry|>prod|>basis|>collect, num_bit=prod(geometry)), (geometry..., :))

# a sample is a bar or a stripe.
is_bar(mat::AbstractMatrix) = sum(abs.(diff(mat, 1)), (1, 2))[] == 0
is_stripe(mat::AbstractMatrix) = sum(abs.(diff(mat, 2)), (1, 2))[] == 0
is_bs(mat::AbstractMatrix) = is_bar(mat) || is_stripe(mat)

"""
    barstripe_pdf(M::Int, N::Int) -> Vector{Float64}

get bar and stripes PDF
"""
function barstripe_pdf(M::Int, N::Int)
    bss = binary_basis(M, N)
    [bss[:,:,i] |> is_bs for i in indices(bss, 3)] |> normalize
end

"""
    bs_configs(M::Int, N::Int) -> BitArray{3}

the total configuration space for bar and stripes.
"""
function bs_configs(M::Int, N::Int)
    bss = binary_basis(M, N)
    slicedim(bss, 3, [bss[:,:,i] |> is_bs for i in indices(bss, 3)])
end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
"""
VWN5 LDA correlation according to Vosko, Wilk, and Nusair, (DOI 10.1139/p80-159).
"""
function energy_per_particle(::Val{:lda_c_vwn}, ρ)
    # Ugly automatically generated code. TODO code from the equations
    x0  = sqrt(5913592351)
    x1  = 2^(2/3)*3^(1/6)
    x4  = (π*ρ)^(1/6)
    x5  = 468052052839577158149*x0/928058053502134909615000000
    x7  = 12500*x4*6^(1/3)
    x11 = 1 / (46593 * x1 * cbrt(π*ρ) + x7 + 323380 * sqrt(π*ρ))

    (  x5*atan(x0*x4/(12500*x1 + 46593*x4)) + 310907*log(x11*x7)/10000000
     + 76037485627899*log(x11*x4*(25000*x1 + 5249*x4)^2/100000)/78468213432500000)
end


text:
############################################################
# joProjection #############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export joProjection, joProjectionException

############################################################

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections
struct joProjection{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
	name::String
	m::Integer
	n::Integer
	info::Info
	geometry::Geometry
end

mutable struct joProjectionException <: Exception
	msg :: String
end

############################################################
## Constructor
"""
    joProjection(info, geometry)

Projection operator for sources/receivers to restrict or inject data at specified locations.\\
`info` is an `Info` structure and `geometry` is a `Geometry` structure with either source or\\
receiver locations.

Examples
========

`F` is a modeling operator of type `joModeling` and `q` is a seismic source of type `joData`:

    Pr = joProjection(info, rec_geometry)
    Ps = joProjection(info, q.geometry)

    dobs = Pr*F*Ps'*q
    qad = Ps*F'*Pr'*dobs

"""
function joProjection(info::Info, geometry::GeometryIC; DDT::DataType=Float32, RDT::DataType=DDT)
	(DDT == Float32 && RDT == Float32) || throw(joProjectionException("Domain and range types not supported"))
	m = 0
	for j=1:length(geometry.xloc)
		m += length(geometry.xloc[j])*geometry.nt[j]
	end
	n = info.n * sum(info.nt)

	return joProjection{Float32,Float32}("restriction operator",m,n,info,geometry)
end

function joProjection(info::Info, geometry::GeometryOOC; DDT::DataType=Float32, RDT::DataType=DDT)
	(DDT == Float32 && RDT == Float32) || throw(joProjectionException("Domain and range types not supported"))
	m = sum(geometry.nsamples)
	n = info.n * sum(info.nt)

	return joProjection{Float32,Float32}("restriction operator",m,n,info,geometry)
end

############################################################
## overloaded Base functions

# conj(joProjection)
conj{DDT,RDT}(A::joProjection{DDT,RDT}) =
	joProjection{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.geometry)

# transpose(joProjection)
transpose{DDT,RDT}(A::joProjection{DDT,RDT}) =
	joProjection{DDT,RDT}("injection operator",A.n,A.m,A.info,A.geometry)

# ctranspose(joProjection)
ctranspose{DDT,RDT}(A::joProjection{DDT,RDT}) =
	joProjection{DDT,RDT}("injection operator",A.n,A.m,A.info,A.geometry)

############################################################
## overloaded Base *(...joProjection...)

# *(joProjection,joData)
function *{ADDT,ARDT,vDT}(A::joProjection{ADDT,ARDT},v::joData{vDT})
	A.n == size(v,1) || throw(joProjectionException("shape mismatch"))
	compareGeometry(A.geometry,v.geometry) == true || throw(joProjectionException("geometry mismatch"))
	jo_check_type_match(ADDT,vDT,join(["DDT for *(joProjection,joData):",A.name,typeof(A),vDT]," / "))
	V = joRHS(A.info,v.geometry,v.data)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joProjection,joData):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(joProjection,joModeling)
function *{ARDT,BDDT,CDT}(A::joProjection{CDT,ARDT},B::joModeling{BDDT,CDT})
	A.n == size(B,1) || throw(joProjectionException("shape mismatch"))
	compareInfo(A.info, B.info) == true || throw(joProjectionException("info mismatch"))
	if typeof(A.geometry) == GeometryOOC
		m = sum(A.geometry.nsamples)
	else
		m = 0; for j=1:B.info.nsrc m+= length(A.geometry.xloc[j])*A.geometry.nt[j] end
	end
	n = B.info.n * sum(B.info.nt)
	return joPDE("joProjection*joModeling",B.info,B.model,A.geometry;options=B.options,DDT=CDT,RDT=ARDT)
end

function *{ARDT,BDDT,CDT}(A::joProjection{CDT,ARDT},B::joModelingAdjoint{BDDT,CDT})
	A.n == size(B,1) || throw(joProjectionException("shape mismatch"))
	compareInfo(A.info, B.info) == true || throw(joProjectionException("info mismatch"))
	if typeof(A.geometry) == GeometryOOC
		m = sum(A.geometry.nsamples)
	else
		m = 0; for j=1:B.info.nsrc m+= length(A.geometry.xloc[j])*A.geometry.nt[j] end
	end
	n = B.info.n * sum(B.info.nt)
	return joPDEadjoint("joProjection*joModelingAdjoint",B.info,B.model,A.geometry;options=B.options,DDT=CDT,RDT=ARDT)
end

############################################################
## Additional overloaded functions

# Subsample Modeling operator
function subsample{ADDT,ARDT}(P::joProjection{ADDT,ARDT}, srcnum)
	geometry = subsample(P.geometry,srcnum)		# Geometry of subsampled data container
	info = Info(P.info.n, length(srcnum), P.info.nt[srcnum])
	return joProjection(info, geometry;DDT=ADDT,RDT=ARDT)
end

getindex(P::joProjection,a) = subsample(P,a)

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
using JLD
using DataFrames
using Dates
using LinearAlgebra
using DifferentialEquations
using Distributions
using Random
using CSV
using Printf
using RCall
using CPUTime
using StatsPlots
include("utils.jl")
include("sir_ode.jl")
include("abc.jl")
include("abcMCMC.jl")
include("abcSMC.jl")

df_post = DataFrame(CSV.File("../data/covid19_bc_post_vaccine.csv"))

 """ Post-vaccine analysis """
 tot_days = size(df_post)[1]
 time_interval = 7

 y_post = hcat([df_post[:S],df_post[:I],df_post[:R]]...)[1:7:tot_days,:]
 y_post = y_post./1e4
 y_dim = size(y_post)[1]

 function data_generator_week(p)
     initial_state = y_post[1,:]; time_window=(1,Float64(tot_days)/time_interval)
     solve_ode(initial_state,time_window,p)
 end

 algo_parameters_abc = (prior = (Truncated(Normal(0.1,0.5),0,Inf),Truncated(Normal(0.2,0.2),0,Inf)),epsilon = 0.1,
 eta= identity_mapping, d= distanceFunction)
 output_abc=ABC(y_post', data_generator_week, algo_parameters_abc, 500)

 algo_parameter_mcmc = (prior = (Truncated(Normal(0.1,0.5),0,Inf),Truncated(Normal(0.2,0.2),0,Inf)),epsilon = 1, eta = identity_mapping,
                         d= compute_norm, proposal = proposal_Normal,
                         proposalRatio = proposalRatio_Normal,sd = (0.25,0.25),
                         thinning = 10000, burn_in=100, verbose=true)

 output_mcmc = ABC_MCMC(y_post', data_generator_week, algo_parameter_mcmc, 500)

 algo_parameters_smc = (prior = (Truncated(Normal(0.1,0.5),0,Inf),Truncated(Normal(0.2,0.2),0,Inf)), time_final=5, eps_list = [10, 5, 1, 0.5, 0.1],
 eta= identity_mapping, d= distanceFunction, kernel=proposal_Normal,
 kernel_density=proposal_Normal_density,
 sd=(0.5,0.5), resample_method=systematic_resample, verbose=true)
 output_smc = ABC_SMC(y_post', data_generator_week, algo_parameters_smc, 500)

 abc_solutions=zeros(3,y_dim,500)
 abc_mcmc_solutions=zeros(3,y_dim,500)
 abc_smc_solutions=zeros(3,y_dim,500)
 for i=1:500
     abc_solutions[:,:,i] = data_generator_week(output_abc[1][i,:])[1:3,:]
     abc_mcmc_solutions[:,:,i] = data_generator_week(output_mcmc[1][i,:])[1:3,:]
     abc_smc_solutions[:,:,i] = data_generator_week(output_smc[1][i,:])[1:3,:]
 end

 save("../results/abc_solution_500_post2.jld","abc_solutions",abc_solutions)
 save("../results/abc_mcmc_solution_500_post2.jld","abc_mcmc_solutions",abc_mcmc_solutions)
 save("../results/abc_smc_solution_500_post2.jld","abc_smc_solutions",abc_smc_solutions)

 alpha_level=1
 beta_plot = density(output_abc[1][:,1],alpha=alpha_level, label="ABC",  title="Predicted β Posterior", xlabel="β", ylabel="Density")
 density!(output_mcmc[1][:,1],alpha=alpha_level, label="ABC-MCMC")
 density!(output_smc[1][:,1],alpha=alpha_level, label="ABC-SMC")
 savefig(beta_plot,"../figs/beta_plot_post2.png")

 gamma_plot = density(output_abc[1][:,2],alpha=alpha_level, label="ABC",  title="Predicted γ Posterior", xlabel="γ", ylabel="Density")
 density!(output_mcmc[1][:,2],alpha=alpha_level, label="ABC-MCMC")
 density!(output_smc[1][:,2],alpha=alpha_level, label="ABC-SMC")
 savefig(gamma_plot,"../figs/gamma_plot_post2.png")

 # ess
ess_output = zeros(3)
ess_output_time = copy(ess_output)
ess_output[1]  = try rcopy(R"mean(mcmcse::ess($(output_abc[1])))") catch; 0 end
ess_output_time[1]  = ess_output[1] / output_abc[2]
ess_output[2]  = try rcopy(R"mean(mcmcse::ess($(output_mcmc[1])))") catch; 0 end
ess_output_time[2]  = ess_output[2] / output_mcmc[2]
ess_output[3]  = try rcopy(R"mean(mcmcse::ess($(output_smc[1])))") catch; 0 end
ess_output_time[3]  = ess_output[3] / output_smc[2]

save("../results/covid/ess_post_2.jld","ess_output",ess_output)
save("../results/covid/ess_cputime_post_2.jld","ess_output_time",ess_output_time)


text:
const ZOTERO = "https://api.zotero.org/"

struct ZoteroClient
    APIKey::String
    user_groupID::String
    user_group_prefix::String
end

function ZoteroClient(APIKey::String, userID::String, user::Bool=true)
    new(APIKey, userID, user ? "users" : "groups")
end

function ZoteroClient(path_to_ids::String=pwd())
    # TODO gives an information message for errors
    path_to_ids = isdir(path_to_ids) ? joinpath(path_to_ids, ".token") : path_to_ids
    return ZoteroClient(readlines(path_to_ids)...)
end

function base_url(client::ZoteroClient)
    return joinpath(ZOTERO, client.user_group_prefix, client.user_groupID)
end

function HTTP.request(client::ZoteroClient, verb::String, url::String, headers::Dict=Dict(), body::String=""; kwargs...)
    merge!(headers, Dict(
                    "Zotero-API-Version" => 3,
                    "Zotero-API-Key" => client.APIKey,
                    )
                )
    response = HTTP.request(verb, joinpath(base_url(client), url), headers, body; kwargs...)
    return response.body
end

function request_json(client::ZoteroClient, verb::String, url::String, headers::Dict=Dict(), body::String=""; kwargs...)
    headers = merge(headers, Dict("Content-Type" => "application/json"))
    body = HTTP.request(client, verb, url, headers, body; kwargs...)
    return JSON.parse(String(body))
end


text:
include("../src/HssMatrices.jl")
using .HssMatrices
using LinearAlgebra
using AbstractTrees
using Plots

# generate Cauchy matrix
K(x,y) = (x-y) > 0 ? 0.001/(x-y) : 2.
#K(x,y) = (x-y) != 0 ? 1/(x-y) : 1.
A = [ K(x,y) for x=-1:0.001:1, y=-1:0.001:1];
C = inv(A);

# test the simple implementation of cluster trees
m, n = size(A)
HssMatrices.setopts(leafsize=64)
rcl = bisection_cluster(1:m)
ccl = bisection_cluster(1:n)

# test compression
hssA = compress(A, rcl, ccl);
println("rel. approximation error with direct compression: ", norm(A - full(hssA))/norm(A))
println("abs. approximation error with direct compression: ", norm(A - full(hssA)))
println("hss-rank with direct compression: ", hssrank(hssA))

# test randomized compression
hssB = randcompress_adaptive(A, rcl, ccl);
println("rel. approximation error with randomized compression: ", norm(A - full(hssB))/norm(A))
println("abs. approximation error with randomized compression: ", norm(A - full(hssB)))
println("hss-rank with randomized compression: ", hssrank(hssB))

# test recompression
hssB = recompress!(copy(hssA))
println("rel. approximation error after recompression: ", norm(A - full(hssB))/norm(A))
println("abs. approximation error after recompression: ", norm(A - full(hssB)))
println("hss-rank after recompression: ", hssrank(hssB))

# test mat-vec
x = randn(size(A,2), 3);
println("rel. error in the matrix-vector products: ", norm(A*x - hssA*x)/norm(A*x))
println("abs. error in the matrix-vector products: ", norm(A*x - hssA*x))

# test the ULV based solver
b = randn(size(A,2), 5);
x = ulvfactsolve(hssA, b);
xcor = A\b;
println("rel. error in the solution of Ax = b: ", norm(x-xcor)/norm(xcor))
println("abs. error in the solution of Ax = b: ", norm(x-xcor))

# test left HSS division
Id(i,j) = Matrix{Float64}(i.*ones(length(j))' .== ones(length(i)).*j')
IdOp = LinearMap{Float64}(n, n, (y,_,x) -> x, (y,_,x) -> x, (i,j) -> Id(i,j))
hssI = randcompress(IdOp, ccl, ccl, 0)
hssC = ldiv!(hssA, hssI)
println("rel. error in the solution of AX = I: ", norm(full(hssC) - C) / norm(C) )
println("abs. error in the solution of AX = I: ", norm(full(hssC) - C) )

# test left HSS division
Id(i,j) = Matrix{Float64}(i.*ones(length(j))' .== ones(length(i)).*j')
IdOp = LinearMap{Float64}(n, n, (y,_,x) -> x, (y,_,x) -> x, (i,j) -> Id(i,j))
hssI = randcompress(IdOp, ccl, ccl, 0)
hssC = rdiv!(hssI, hssA)
println("rel. error in the solution of XA = I: ", norm(full(hssC) - C) / norm(C) )
println("abs. error in the solution of XA = I: ", norm(full(hssC) - C) )

# desequilibrate cluster trees and try again
hssA.A11 = prune_leaves!(hssA.A11)
hssI = randcompress(IdOp, ccl, ccl, 0)
hssI.A11 = prune_leaves!(hssI.A11)
hssC = ldiv!(hssA, hssI)
println("rel. error in the solution of XA = I: ", norm(full(hssC) - C) / norm(C) )
println("abs. error in the solution of XA = I: ", norm(full(hssC) - C) )

# # test computation of generators
# U1, V2 = generators(hssA, (1,2))
# A12 = A[1:hssA.m1,hssA.n1+1:end];
# println(norm(A12 - U1*hssA.B12*V2'))
# U2, V1 = generators(hssA, (2,1))
# A21 = A[hssA.m1+1:end,1:hssA.n1];|
# println(norm(A21 - U2*hssA.B21*V1'))

# #println(typeof(hssA))

# test the ULV based solver
# x = ulvfactsolve(hssA, b);

# @time x = ulvfactsolve(hssA, b);

# # test plotting
# plot = plotranks(hssA)


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")
##

include("sp_model.jl")
include("plot_utils.jl")

##

using DataFrames
using CSV

timesteps = 1:365*24

pv = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps, 3]
wind = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps, 4]
demand = CSV.read("timeseries/basic_example.csv", DataFrame)[timesteps, 2]

##

pars = copy(default_es_pars)

##

using Statistics
average_hourly_demand = mean(demand)

pars[:flexible_demand] = average_hourly_demand * 100. # Over the year, we have 100 demand hours that are shiftable freely...
pars[:recovery_time] = 24
pars[:c_storage] = 100.
pars[:c_pv] = 300.
pars[:c_wind] = 800.
pars[:c_sto_op] = 0.00001

es = define_energy_system(pv, wind, demand; p = pars)

##
n = 100
F_max = average_hourly_demand * 0.1 # Have ancillary services equal to 10% of our typical demand
t_max = length(pv) - es.parameters[2].defaults[:recovery_time]
scens = simple_flex_sampler(n, F_max, t_max)

##

using Cbc

sp = instantiate(es, scens, optimizer = Cbc.Optimizer)

# using GLPK

# sp = instantiate(es, scens, optimizer = GLPK.Optimizer)

##

optimize!(sp)

##

println("Objective value: $(objective_value(sp))")

println("PV: $(round(value.(sp[1, :u_pv]); digits = 2)) -- Wind: $(round(value.(sp[1, :u_wind]); digits = 2)) -- Battery: $(round(value.(sp[1, :u_storage]); digits = 2))")

# for s in 1:minimum((n, 10))
#     println("Second stage objective value in scenario $s: $(objective_value(sp, s))")
# end

##

plot_results(sp, pv, wind, demand, s = 8)

##


text:
module Mex

using MATLAB, Libdl

export jl_mex, input, call_matlab, is_interrupt_pending, check_for_interrupt

# --- Initialization --- #

# MATLAB mex interface
const libmex = Ref{Ptr{Cvoid}}()
const mex_call_matlab_with_trap = Ref{Ptr{Cvoid}}()

# MATLAB interrupt handling (from undocumented MATLAB)
const libut = Ref{Ptr{Cvoid}}()
const ut_is_interrupt_pending = Ref{Ptr{Cvoid}}()

function __init__()
    libmex[] = Libdl.dlopen(joinpath(MATLAB.matlab_libpath(), "libmex"), Libdl.RTLD_GLOBAL)
    mex_call_matlab_with_trap[] = Libdl.dlsym(libmex[], :mexCallMATLABWithTrap)

    libut[]  = Libdl.dlopen(joinpath(MATLAB.matlab_libpath(), "libut"), Libdl.RTLD_GLOBAL)
    ut_is_interrupt_pending[] = Libdl.dlsym(libut[], :utIsInterruptPending)
end

include("exceptions.jl")
include("juliatomatlab.jl")
include("matlabtojulia.jl")

# Patch MATLAB.jl's handling of tuples
function MATLAB.mxarray(t::Tuple)
    pm = mxcellarray(length(t))
    for i = 1:length(t)
        set_cell(pm, i, mxarray(t[i]))
    end
    return pm
end

# Patch MATLAB.jl's handling of reinterpreted arrays
function MATLAB.mxarray(a::Base.ReinterpretArray{T}) where T<:MATLAB.MxRealNum
    mx = mxarray(T, size(a))
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt), data_ptr(mx), parent(a), length(a)*sizeof(T))
    return mx
end

function MATLAB.mxarray(a::Base.ReinterpretArray{T}) where T<:MATLAB.MxComplexNum
    mx = mxarray(T, size(a))
    na = length(a)
    rdat = unsafe_wrap(Array, MATLAB.real_ptr(mx), na)
    idat = unsafe_wrap(Array, MATLAB.imag_ptr(mx), na)
    @inbounds for i = 1:na
        rdat[i] = real(a[i])
        idat[i] = imag(a[i])
    end
    return mx
end

# Patch MATLAB.jl's handling of reshaped arrays
function MATLAB.mxarray(a::Base.ReshapedArray{T}) where T<:MATLAB.MxRealNum
    mx = mxarray(T, size(a))
    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt), data_ptr(mx), parent(a), length(a)*sizeof(T))
    return mx
end

function MATLAB.mxarray(a::Base.ReshapedArray{T}) where T<:MATLAB.MxComplexNum
    mx = mxarray(T, size(a))
    na = length(a)
    rdat = unsafe_wrap(Array, MATLAB.real_ptr(mx), na)
    idat = unsafe_wrap(Array, MATLAB.imag_ptr(mx), na)
    @inbounds for i = 1:na
        rdat[i] = real(a[i])
        idat[i] = imag(a[i])
    end
    return mx
end

end # module


text:
using Pkg
using LinearAlgebra
using DataFrames
import Statistics.mean
import Statistics.std
import Statistics.cor
using Printf
using CSV

include(pwd()*"/src/SDEMEM OU neuron data/ou_sdemem.jl")
include(pwd()*"/src/SDEMEM OU neuron data/mcmc.jl")

seed = 100

M_subjects = 1

y, prior_parameters_η, prior_parameters_σ_ϵ = set_up(M_subjects, seed)

startval_ϕ = ones(M_subjects,3)

for j = 1:3

    μ_0_j, M_0_j, α_j, β_j = prior_parameters_η[j,:]

    startval_ϕ[:,j,1] = μ_0_j*ones(M_subjects)

end

startval_log_σ_ϵ = log(0.5)

# set time step
dt = diff(y[1].mSec)[1]

# posterior mean ests from kalman filter
log_σ_ϵ = -3.60203
ϕ  = [-3.26986 -0.888105 -0.780767]
nbr_particles_cor = 10

# estimate correlation
nbr_pf_eval_corr = 100

ll1 = zeros(nbr_pf_eval_corr)
ll2 = zeros(nbr_pf_eval_corr)

ρ = 0.99

for i = 1:nbr_pf_eval_corr

    nbr_particles = nbr_particles_cor*ones(Int64, M_subjects)

    # set number of obs for each data set
    T_vec = zeros(Int64, M_subjects)
    for i = 1:M_subjects; T_vec[i] = length(y[i].mV); end

    u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
    u_resample_old = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
        append!(u_resample_old, [randn(T_vec[i],2)])
    end

    u_prop_new = [randn(nbr_particles[1], T_vec[1]+1)]
    u_resample_new = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_new, [randn(nbr_particles[i],T_vec[i]+1)])
        append!(u_resample_new, [randn(T_vec[i],2)])
    end

    # correlat with old standard normal numbers
    u_prop_tilde = ρ*u_prop_old + sqrt(1-ρ^2)*u_prop_new
    u_resample_tilde = ρ*u_resample_old + sqrt(1-ρ^2)*u_resample_new

    ll1[i] = sum(bridge_pf(y, exp(log_σ_ϵ), ϕ, dt, u_prop_old, u_resample_old, nbr_particles,true))
    ll2[i] = sum(bridge_pf(y, exp(log_σ_ϵ), ϕ, dt, u_prop_tilde, u_resample_tilde, nbr_particles,true))

end

corr_ll = cor(ll1,ll2)

var_loglik_target = 2.16^2/(1-corr_ll^2)

#nbr_pf_eval = 100
#nbr_particles = 50

function run_bridge_pf(particels)

    nbr_particles = particels*ones(Int64, M_subjects)
    ρ = 0.0

    # set number of obs for each data set
    T_vec = zeros(Int64, M_subjects)
    for i = 1:M_subjects; T_vec[i] = length(y[i].mV); end

    u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
    u_resample_old = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
        append!(u_resample_old, [randn(T_vec[i],2)])
    end

    run_time = @elapsed ll1 =  bridge_pf(y, exp(log_σ_ϵ), ϕ, dt, u_prop_old, u_resample_old, nbr_particles, true)

    return ll1, run_time

end

function loglik_est(nbr_particles_pf, nbr_pf_eval)

    ll = zeros(nbr_pf_eval)

    for i in 1:nbr_pf_eval
        ll_part, ~ = run_bridge_pf(nbr_particles_pf)
        ll[i] = sum(ll_part)
    end

    return ll

end

nbr_pf_eval = 100
nbr_particles_pf = [1, 2, 5, 10, 20, 50, 100]

loglik_matrix = zeros(nbr_pf_eval, length(nbr_particles_pf))
run_times = zeros(length(nbr_particles_pf))

i = 1
run_times[i]  = @elapsed loglik_matrix[:,i] = loglik_est(nbr_particles_pf[i], nbr_pf_eval)

Random.seed!(1)

for i in 1:length(nbr_particles_pf)
    @printf "Current nbr of particels: %.2f\n" nbr_particles_pf[i]
    run_times[i]  = @elapsed loglik_matrix[:,i]  = loglik_est(nbr_particles_pf[i], nbr_pf_eval)
end

using PyPlot

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, mean(loglik_matrix, dims = 1)[:], "--*")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Avg. loglik est. val.")
#PyPlot.savefig("figures/avg_loglik_val_vs_nbr_particels_cpf_0999.png")

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, var(loglik_matrix, dims = 1)[:], "--*")
PyPlot.plot(nbr_particles_pf, var_loglik_target*ones(length(nbr_particles_pf),1), "k")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Var of loglik est.")
#PyPlot.savefig("figures/loglik_var_vs_nbr_particles_cpf_0999.png")

PyPlot.figure(figsize=(20,10))
PyPlot.plot(nbr_particles_pf, run_times/100, "--*")
PyPlot.xlabel("Nbr particles.")
PyPlot.ylabel("Run time (sec)")
#PyPlot.savefig("figures/run_time_vs_nbr_particls_cpf_0999.png")

max_stars_repo_path: src/SDEMEM OU neuron data/tune_nbr_particles_cpf_bridge_one_subject.jl


text:
function clip{T<:FloatingPoint}(val::T; low::T=1e-10, high::T=1.0-1e-10)
	if val < low
		return low
	elseif val > high
		return high
	else
		return val
	end
end


text:
module PkgOptionsDemo

using PkgOptionsLoader

const pkgoptions = PkgOptionsLoader.@load

end # module


text:
module SV
using Econometrics, Random, Statistics, LinearAlgebra
# Utilities
include("HAR.jl")
include("aux_stat.jl")
include("SVmodel.jl")
include("SVmoments.jl")
include("logL.jl")
include("prior.jl")
include("proposal1.jl")
include("proposal2.jl")
export HAR, aux_stat, SVmodel, SVmoments, logL, prior, proposal1, proposal2
end

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Examples that demonstrate the diffs reported by EnhancedTestSet

Acknowledgements
----------------
* This example was taken more or less directly from taken directly from the
  TestSetExtensions package developed by Spencer Russell and his collaborators
  (https://github.com/ssfrr/TestSetExtensions.jl).
"""

# --- Imports

# Standard library
using Test

# Local modules
using TestTools.jltest: EnhancedTestSet

# --- Tests

@testset EnhancedTestSet "Diff Examples" begin
    @test [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]

    @test """Lorem ipsum dolor sit amet,
             consectetur adipiscing elit, sed do
             eiusmod tempor incididunt ut
             labore et dolore magna aliqua.
             Ut enim ad minim veniam, quis nostrud
             exercitation ullamco aboris.""" == """Lorem ipsum dolor sit amet,
                                                   consectetur adipiscing elit, sed do
                                                   eiusmod temper incididunt ut
                                                   labore et dolore magna aliqua.
                                                   Ut enim ad minim veniam, quis nostrud
                                                   exercitation ullamco aboris."""

    @test Dict(:foo => "bar", :baz => [1, 4, 5], :biz => nothing) ==
        Dict(:baz => [1, 7, 5], :biz => 42)
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))
push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using ClimaComms
using ClimaCommsMPI

include(joinpath(@__DIR__, "..", "..", "..", "test", "stencil.jl"))

stencil_test(ClimaCommsMPI.MPICommsContext)


text:
struct SymEGRQSMatrix{T,UT<:AbstractArray,VT<:AbstractArray,dT<:AbstractArray} <: AbstractMatrix{T}
    Ut::UT
    Vt::VT
    d::dT
	n::Int
	p::Int
	function SymEGRQSMatrix{T,UT,VT,dT}(Ut,Vt,d,n,p) where
			 {T,UT<:AbstractArray,VT<:AbstractArray,dT<:AbstractArray}
        Up, Un = size(Ut)
		Vp, Vn = size(Vt)
		(Un == Vn && Up == Vp && Vn == length(d)) || throw(DimensionMismatch())
    	new(Ut,Vt,d,n,p)
	end
end

# Constuctor
SymEGRQSMatrix(Ut::AbstractArray{T,N}, Vt::AbstractArray{T,N}, d::AbstractArray{T,M}) where {T,N,M}=
	SymEGRQSMatrix{T,typeof(Ut),typeof(Vt),typeof(d)}(Ut,Vt,d,size(Vt,2),size(Vt,1));

########################################################################
#### Helpful properties. Not nessecarily computionally efficient    ####
########################################################################
Matrix(K::SymEGRQSMatrix) = tril(K.Ut'*K.Vt) + triu(K.Vt'*K.Ut,1) + Diagonal(K.d)

size(K::SymEGRQSMatrix) = (K.n, K.n)

function getindex(K::SymEGRQSMatrix{T}, i::Int, j::Int) where T
	i > j && return dot(K.Ut[:,i],K.Vt[:,j])
	j == i && return dot(K.Vt[:,i],K.Ut[:,j]) + K.d[i]
	return dot(K.Vt[:,i],K.Ut[:,j])
end

Base.propertynames(F::SymEGRQSMatrix, private::Bool=false) =
    (private ? fieldnames(typeof(F)) : ())

########################################################################
#### Linear Algebra routines	                                    ####
########################################################################

#### Matrix-matrix product ####
function dss_mul_mat!(Y::Array, K::SymEGRQSMatrix, X::Array)
	Ut = K.Ut
	Vt = K.Vt
	d = getfield(K,:d)
    p, n = size(Ut);
    mx = size(X,2);
    Vbar = zeros(p,mx);
    Ubar = Ut*X;
    @inbounds for i = 1:n
        tmpV = Vt[:,i];
        tmpU = Ut[:,i];
        tmpX = X[i:i,:];
        Ubar -= tmpU .* tmpX;
        Vbar += tmpV .* tmpX;
        Y[i,:] = tmpU'*Vbar + tmpV'*Ubar + d[i]*tmpX;
    end
	return Y
end

#### Log-determinant ####
logdet(K::SymEGRQSMatrix) = 2.0*logdet(cholesky(K))

#### Determinant ####
det(K::SymEGRQSMatrix) = det(cholesky(K))^2

mul!(y::AbstractVecOrMat, K::SymEGRQSMatrix, x::AbstractVecOrMat) =
	dss_mul_mat!(y,K,x)
mul!(y::AbstractVecOrMat, K::Adjoint{<:Any,<:SymEGRQSMatrix}, x::AbstractVecOrMat) =
	dss_mul_mat!(y,K.parent,x)
function (\)(K::SymEGRQSMatrix, x::AbstractVecOrMat)
	L = cholesky(K);
	return L'\(L\x)
end

max_stars_repo_name: mipals/SymEGRSSMatrices


text:
type HermiteRadialField{N, T, F <: Function} <: AbstractScalarField{N}
	alphas::Vector{T}
	betas::Vector{SVector{N, T}}
	points::Vector{SVector{N, T}}
	phi::F
end

function HermiteRadialField{Dimension, T, F <: Function}(points::Vector{SVector{Dimension, T}}, normals::Vector{SVector{Dimension, T}}, phi::F=XCubed())
	@assert length(points) == length(normals)

	num_points = length(points)

	A = Array{T}(num_points * (Dimension + 1), num_points * (Dimension + 1))
	b = Array{T}(num_points * (Dimension + 1))
	# u = Array{T}(Dimension)
	# v = Array{T}(Dimension)
    dphi = derivative(phi)
    ddphi = derivative(dphi)

	for point_index = 1:num_points
		row = (point_index - 1) * (1 + Dimension) + 1
		for k = 1:num_points
			col = (k - 1) * (1 + Dimension) + 1
			u = points[point_index] - points[k]
			# for i = 1:Dimension
			# 	u[i] = points[i, point_index] - points[i, k]
			# end
			# u = points[:,point_index] - points[:,k]
			n = norm(u)
			if n <= 1e-9
				A[row + (0:Dimension), col + (0:Dimension)] = 0
			else
				f = phi(n)
                df = dphi(n)
                ddf = ddphi(n)
				df_over_n = df / n
				v = df_over_n * u
				# for i = 1:Dimension
				# 	v[i] = df_over_n * u[i]
				# end

				A[row, col] = f
				for i = 1:Dimension
					A[row, col + i] = v[i]
					A[row + i, col] = v[i]
				end
				scaling = (ddf - df_over_n) / (n^2)
				for i = 1:Dimension
					for j = 1:Dimension
						A[row + i, col + j] = scaling * u[i] * u[j]
					end
				end
				# A[row + (1:Dimension), col + (1:Dimension)] = (ddf - df_over_n) / (n^2) * (u * u')
				for i = 1:Dimension
					A[row + i, col + i] += df_over_n
				end
			end
		end

		b[row] = 0
		for i = 1:Dimension
			b[row + i] = normals[point_index][i]
		end
	end

	y = A \ b
	y = reshape(y, Dimension + 1, num_points)
	alphas = vec(y[1,:])
	betas = [SVector{Dimension, T}(y[2:end, i]) for i in 1:num_points]
	HermiteRadialField{Dimension, T, F}(alphas, betas, points, phi)
end

function (field::HermiteRadialField{Dimension, T}){Dimension, T}(x::AbstractVector{T})
	value::T = zero(T)
	num_points = length(field.points)
	@assert length(x) == Dimension
	@assert length(field.alphas) == num_points
	@assert length(field.betas) == length(field.points)
    dphi = derivative(field.phi)
	@inbounds for i = eachindex(field.points)
		u = x - field.points[i]
		n = norm(u)
		if n > 0
			z = sum(field.betas[i] .* u)
			p = field.phi(n)
			dp = dphi(n)
			value += field.alphas[i] * p
			value += dp / n * z
		end
	end
	value
end

function gradient{Dimension, T, PhiType}(field::HermiteRadialField{Dimension, T, PhiType}, x::AbstractVector{T})
	num_points = length(field.points)
	@assert length(x) == Dimension
	@assert length(field.alphas) == num_points
	@assert length(field.betas) == num_points
	g = zeros(SVector{Dimension, T})
    dphi = derivative(field.phi)
    ddphi = derivative(dphi)
	for i = eachindex(field.points)
		u = x - field.points[i]
		n = norm(u)

		if n > 1e-5
			uhat = u ./ n
            df = dphi(n)
            ddf = ddphi(n)
			alpha_df = field.alphas[i] * df
			beta_uhat = sum(field.betas[i] .* uhat)

			g += alpha_df .* uhat + beta_uhat * (ddf * uhat - u * df / n^2) + field.betas[i] * df / n
		end
	end
	g
end

gradient(field::HermiteRadialField) = x -> gradient(field, x)

# function bounds(field::HermiteRadialField)
# 	lb = Vec(minimum(field.points))
# 	ub = Vec(maximum(field.points))
# 	HyperRectangle(lb, ub - lb)
# end


text:
using Documenter, TensorFunctions

makedocs(;
    modules=[TensorFunctions],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/ho-oto/TensorFunctions.jl/blob/{commit}{path}#L{line}",
    sitename="TensorFunctions.jl",
    authors="ho-oto",
    assets=String[],
)

deploydocs(;
    repo="github.com/ho-oto/TensorFunctions.jl",
)


text:
using PUDLData
using Test

@testset "PUDLData.jl" begin
    @test all(startswith("pudl"), pudl_datasets())
end


text:
"""
```
solve(m::AbstractDSGEModel; apply_altpolicy = false)
```

Driver to compute the model solution and augment transition matrices.

### Inputs

- `m`: the model object

## Keyword Arguments

- `apply_altpolicy::Bool`: whether or not to solve the model under the
  alternative policy. This should be `true` when we solve the model to
  forecast, but `false` when computing smoothed historical states (since
  the past was estimated under the baseline rule).

### Outputs
 - TTT, RRR, and CCC matrices of the state transition equation:
    ```
    S_t = TTT*S_{t-1} + RRR*ϵ_t + CCC
    ```
"""
function solve(m::AbstractDSGEModel; apply_altpolicy = false, verbose::Symbol = :high)

    altpolicy_solve = alternative_policy(m).solve

    if altpolicy_solve == solve || !apply_altpolicy
        # Get equilibrium condition matrices
        Γ0, Γ1, C, Ψ, Π  = eqcond(m)

        # Solve model
        TTT_gensys, CCC_gensys, RRR_gensys, eu = gensys(Γ0, Γ1, C, Ψ, Π, 1+1e-6,
                                                        verbose = verbose)

        # Check for LAPACK exception, existence and uniqueness
        if eu[1] != 1 || eu[2] != 1
            throw(GensysError())
        end

        TTT_gensys = real(TTT_gensys)
        RRR_gensys = real(RRR_gensys)
        CCC_gensys = real(CCC_gensys)

        # Augment states
        TTT, RRR, CCC = augment_states(m, TTT_gensys, RRR_gensys, CCC_gensys)

    else
        # Change the policy rule
        TTT, RRR, CCC = altpolicy_solve(m)
    end
    return TTT, RRR, CCC
end

"""
```
GensysError <: Exception
```
A `GensysError` is thrown when:

1. Gensys does not give existence and uniqueness, or
2. A LAPACK error was thrown while computing the Schur decomposition of Γ0 and Γ1

If a `GensysError`is thrown during Metropolis-Hastings, it is caught by
`posterior`.  `posterior` then returns a value of `-Inf`, which
Metropolis-Hastings always rejects.

### Fields

* `msg::String`: Info message. Default = \"Error in Gensys\"
"""
mutable struct GensysError <: Exception
    msg::String
end
GensysError() = GensysError("Error in Gensys")
Base.showerror(io::IO, ex::GensysError) = print(io, ex.msg)


text:
"""
    All objects that can be assigned degrees of freedom should be subtypes of this.
"""
abstract DofObject
    # id::Int64
    # dofs::Array{T}
    # processor_id::Int64

"Return the list of dofs on this DofObject"
function dofs(dof_object::DofObject)
    return dof_object.dofs
end

"Return the processor_id this DofObject is assigned to"
function processor_id(dof_object::DofObject)
    return dof_object.processor_id
end


text:
"Holds the tableau of a Runge-Kutta method."
abstract type AbstractTableauRK{T <: Real} <: AbstractTableau{T} end

"Holds the tableau of an explicit Runge-Kutta method."
abstract type AbstractTableauERK{T} <: AbstractTableauRK{T} end

"Holds the tableau of an implicit Runge-Kutta method."
abstract type AbstractTableauIRK{T} <: AbstractTableauRK{T} end

"Holds the tableau of a partitioned Runge-Kutta method."
abstract type AbstractTableauPRK{T} <: AbstractTableauRK{T} end

Base.:(==)(tab1::T1, tab2::T2) where {T1 <: Union{AbstractTableauERK,AbstractTableauIRK},
                                      T2 <: Union{AbstractTableauERK,AbstractTableauIRK}} = (tab1.q == tab2.q && order(tab1) == order(tab2) && nstages(tab1) == nstages(tab2))

Base.:(==)(tab1::AbstractTableauPRK, tab2::AbstractTableauPRK) = (tab1.q == tab2.q && tab1.p == tab2.p && order(tab1) == order(tab2) && nstages(tab1) == nstages(tab2))

Base.isequal(tab1::AbstractTableauRK{DT1}, tab2::AbstractTableauRK{DT2}) where {DT1, DT2} = (tab1 == tab2 && DT1 == DT2 && typeof(tab1) == typeof(tab2))

"Reads and parses Tableau metadata from file."
function readTableauRKHeaderFromFile(file)
    f = open(file, "r")
    header = readline(f)
    close(f)

    if header[1] == '#'
        header = split(header[2:end])
    else
        header = ()
    end

    if length(header) ≥ 1
        O = Base.parse(Int, header[1])
    else
        O = 0
    end

    if length(header) ≥ 2
        S = Base.parse(Int, header[2])
    else
        S = 0
    end

    if length(header) ≥ 3
        T = Core.eval(Main, Meta.parse(header[3]))
    else
        T = Float64
    end

    return O, S, T
end

"Write Runge-Kutta tableau to file."
function writeTableauToFile(dir::AbstractString, tab::AbstractTableauRK{T}) where {T}
    # tab_array = zeros(T, S+1, S+1)
    # tab_array[1:S, 2:S+1] = tab.q.a
    # tab_array[S+1, 2:S+1] = tab.q.b
    # tab_array[1:S, 1] = tab.q.c
    # tab_array[S+1, 1] = tab.q.order

    tab_array = zeros(T, tab.q.s, tab.q.s+2)
    tab_array[1:tab.q.s, 1] = tab.q.c
    tab_array[1:tab.q.s, 2] = tab.q.b
    tab_array[1:tab.q.s, 3:tab.q.s+2] = tab.q.a

    header = string("# ", tab.q.o, " ", tab.q.s, " ", T, "\n")
    file   = string(dir, "/", tab.q.name, ".tsv")

    get_config(:verbosity) > 1 ? @info("Writing Runge-Kutta tableau $(tab.q.name) with $(tab.q.s) stages and order $(tab.q.o) to file\n$(file)") : nothing

    f = open(file, "w")
    write(f, header)
    writedlm(f, float(tab_array))
    close(f)
    # TODO Write data in original format (e.g., Rational).
end

# TODO function writeTableauToFile{Name, T}(dir::AbstractString, tab::AbstractTableauPRK{Name, T})


text:
export sinkhorn, dist, empirical_sinkhorn, dtw
# roadmap: implement all related functions from https://github.com/rflamary/POT/blob/master/ot/bregman.py

"""
    sinkhorn(a::Union{PyObject, Array{Float64}}, b::Union{PyObject, Array{Float64}}, M::Union{PyObject, Array{Float64}};
    reg::Float64 = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn")

Computes the optimal transport with Sinkhorn algorithm.
The implementation are adapted from https://github.com/rflamary/POT.
"""
function sinkhorn(a::Union{PyObject, Array{Float64}}, b::Union{PyObject, Array{Float64}}, M::Union{PyObject, Array{Float64}};
    reg::Union{PyObject,Float64} = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn", return_optimal::Bool=false)
    if isa(a, Array)
        @assert sum(a)≈1.0
        @assert all(a .>= 0)
    end
    if isa(b, Array)
        @assert sum(b)≈1.0
        @assert all(b .>= 0)
    end
    a = convert_to_tensor(a)
    b = convert_to_tensor(b)
    M = convert_to_tensor(M)
    reg = convert_to_tensor(reg)
    iter = convert_to_tensor(iter)
    tol = convert_to_tensor(tol)
    sk = load_system_op(COLIB["sinkhorn_knopp"]...; multiple=true)
    if method=="sinkhorn"
        if return_optimal
            return sk(a,b,M,reg,iter,tol, constant(0))
        end
        return sk(a,b,M,reg,iter,tol, constant(0))[2]
    elseif method=="lp"
        if !haskey(COLIB, "ot_network")
            install("OTNetwork", force=true)
        end
        lp = load_system_op("ot_network"; multiple=true)
        if return_optimal
            return lp(a, b, M, iter)
        end
        return lp(a, b, M, iter)[2]
    else
        error("$method not implemented")
    end
end

"""
    empirical_sinkhorn(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}}, dist::Function;
    reg::Union{PyObject,Float64} = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn")

Computes the empirical Wasserstein distance with sinkhorn algorithm.
The implementation are adapted from https://github.com/rflamary/POT.
"""
function empirical_sinkhorn(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}};
    reg::Union{PyObject,Float64} = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn", dist::Function=dist, return_optimal::Bool=false)
    M = dist(x, y)
    a = tf.ones(tf.shape(x)[1], dtype=tf.float64)/cast(Float64, tf.shape(x)[1])
    b = tf.ones(tf.shape(y)[1], dtype=tf.float64)/cast(Float64, tf.shape(y)[1])
    sinkhorn(a, b, M; reg=reg, iter=iter, tol=tol, method=method, return_optimal=return_optimal)
end

"""
    dist(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}}, order::Union{Int64, PyObject}=2)

Computes the distance function with norm `order`. `dist` returns a ``n\\times m`` matrix, where ``x\\in \\mathbb{R}^{n\\times d}`` and
``y\\in \\mathbb{R}^{m\\times d}``, and the return ``M\\in \\mathbb{R}^{n\\times m}``
```math
M_{ij} = ||x_i - y_j||_{o}
```
"""
function dist(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}}, order::Union{Int64, PyObject}=2)
    x = convert_to_tensor(x)
    y = convert_to_tensor(y)
    order = convert_to_tensor(order, dtype=Float64)
    x = tf.expand_dims(x, axis=1)
    y = tf.expand_dims(y, axis=0)
    sum(abs(x-y)^order, dims=3)^(1.0/order)
end

"""
    dtw(s::Union{PyObject, Array{Float64}}, t::Union{PyObject, Array{Float64}},
        use_fast::Bool = false)

Computes the dynamic time wrapping (DTW) distance between two time series `s` and `t`.
Returns the distance and path. `use_fast` specifies whether fast algorithm is used. Note
fast algorithm may not be accurate.
"""
function dtw(s::Union{PyObject, Array{Float64}}, t::Union{PyObject, Array{Float64}},
            use_fast::Bool = false)
    use_fast = Int32(use_fast)
    if !haskey(COLIB, "fast_dtw")
        install("FastDTW", force=true)
    end
    dtw_ = load_system_op(COLIB["fast_dtw"]...; multiple=true)
    s,t, use_fast = convert_to_tensor([s,t,use_fast], [Float64,Float64, Int32])
    cost, path = dtw_(s,t,use_fast)
    return cost, path + 1
end


text:
using Test, TestSetExtensions, SafeTestsets

include(joinpath(@__DIR__,"test_models.jl"))
module TestSetupTeardown

  using SearchLight
  using SearchLightOracle

  export prepareDbConnection, tearDown, connection_file

  connection_file = "oracle_connection.yml"

  function prepareDbConnection()

    conn_info_oracle = SearchLight.Configuration.load(connection_file)
    conn = SearchLight.connect(conn_info_oracle)
    return conn
  end

  function tearDown(conn)
    if conn !== nothing
        ######## Dropping used tables
        SearchLight.Migration.drop_migrations_table()

        # # insert tables you use in tests here
        # tables = ["Book","BookWithIntern","Callback","Author","BookWithAuthor"]

        # # obtain tables exists or not, if they does drop it
        # wheres = join(map(x -> uppercase(string("'", uppercase(SearchLight.Inflector.to_plural(x)), "'")), tables), " , ", " , ")
        # queryString = string("SELECT table_name FROM user_tables where table_name in ($wheres)")
        # result = SearchLight.query(queryString)
        # for item in eachrow(result)
        #     try
        #         SearchLight.Migration.drop_table((item[1]))
        #     catch ex
        #         @show "Table $(item[1]) doesn't exist"
        #     end
        # end

        SearchLight.disconnect(conn)
        rm(SearchLight.config.db_migrations_folder, force=true, recursive=true)
    end
  end

end

@safetestset "Gesamt test" begin

    @safetestset "Connection test" begin
      using SearchLight
      using SearchLightOracle
      using Main.TestSetupTeardown

      connection_file = "oracle_connection.yml"

      conn = prepareDbConnection()
      @test conn !== nothing

      tearDown(conn)

    end

    @safetestset "Orcalce connection and infos" begin
        using SearchLight
        using SearchLightOracle
        using Main.TestSetupTeardown

        conn = prepareDbConnection()

        infoDB = SearchLightOracle.connectionInfo()

        conn_info_oracle = SearchLight.Configuration.load(connection_file)
        itemsToTest=["host","username","database"]

        try
          for info in itemsToTest
              @test uppercase(conn_info_oracle[info]) == uppercase(infoDB[info])
          end
        catch ex
            println("Fehler")
        end

        tearDown(conn)
    end

    @safetestset "Oracle query" begin
        using SearchLight
        using SearchLightOracle
        using SearchLight.Configuration
        using SearchLight.Migrations
        using Main.TestSetupTeardown

        conn = prepareDbConnection()

        queryString = string("select table_name from user_tables where table_name = upper('$(SearchLight.SEARCHLIGHT_MIGRATIONS_TABLE_NAME)')")

        @test isempty(SearchLight.query(queryString, conn)) == true

      # create migrations_table
        SearchLight.Migration.create_migrations_table()

        @test Array(SearchLight.query(queryString, conn))[1] == uppercase(SearchLight.SEARCHLIGHT_MIGRATIONS_TABLE_NAME)

        tearDown(conn)

    end;

    @safetestset "Models and tableMigration" begin
        using SearchLight
        using SearchLightOracle
        using Main.TestSetupTeardown
        using Main.TestModels

      ## establish the database-connection
        conn = prepareDbConnection()

      ## create migrations_table
        SearchLight.Migration.create_migrations_table()

      ## make Table "Book"
        SearchLight.Generator.new_table_migration(Book)
        SearchLight.Migration.up()

        testBook = Book(title="Faust", author="Goethe")

        @test testBook.author == "Goethe"
        @test testBook.title == "Faust"
        @test typeof(testBook) == Book
        @test isa(testBook, AbstractModel)

        try
          testBook |> SearchLight.save
        catch ex
          errorCode = ex
        end

        # @test (testBook |> SearchLight.save) === true

      ############ tearDown ##################
        SearchLight.Migration.down()
        tearDown(conn)

    end

    @safetestset "Model Store and Query models without intern variables" begin
        using SearchLight
        using SearchLightOracle
        using Main.TestModels

        using Main.TestSetupTeardown

      ## establish the database-connection
        conn = prepareDbConnection()

      ## create migrations_table
        SearchLight.Migration.create_migrations_table()

      ## make Table "Book"
        SearchLight.Generator.new_table_migration(Book)
        SearchLight.Migration.up()

        testBooks = Book[]

      ## prepare the TestBooks
        for book in TestModels.seed()
            push!(testBooks, Book(title=book[1], author=book[2]))
        end

        @test testBooks |> SearchLight.save == true

        booksReturn = SearchLight.find(Book)

        @test size(booksReturn) == (5,)

      ############ tearDown ##################
        SearchLight.Migration.all_down!!(confirm=false)
        tearDown(conn)

    end

    @safetestset "Query and Models with intern variables" begin
      using Test
      using SearchLight
      using SearchLightOracle
      using Main.TestSetupTeardown
      using Main.TestModels

      ## establish the database-connection
      conn = prepareDbConnection()

      ## make Table "BooksWithInterns"
      SearchLight.Migration.create_migrations_table()
      SearchLight.Generator.new_table_migration(BookWithInterns)
      SearchLight.Migration.up()

      booksWithInterns = BookWithInterns[]

      ## prepare the TestBooks
      for book in TestModels.seed()
          push!(booksWithInterns, BookWithInterns(title=book[1], author=book[2]))
      end

      testItem = BookWithInterns(author="Alexej Tolstoi", title="Krieg oder Frieden")

      savedTestItem = SearchLight.save(testItem)
      @test savedTestItem === true

      savedTestItems = booksWithInterns |> save
      @test savedTestItems === true

      idTestItem = SearchLight.save!(testItem)
      @test idTestItem.id !== nothing
      @test idTestItem.id.value  > 0

      resultBooksWithInterns = booksWithInterns |> save!

      fullTestBooks = find(BookWithInterns)
      @test isa(fullTestBooks, Array{BookWithInterns,1})
      @test length(fullTestBooks) > 0

      ############ tearDown ##################
      SearchLight.Migration.all_down!!(confirm=false)
      tearDown(conn)

    end

    @safetestset "Saving and Reading fields and datatbase columns are different" begin
        using SearchLight
        using SearchLightOracle
        using Main.TestSetupTeardown
        using Main.TestModels

        conn = prepareDbConnection()

        SearchLight.Migration.create_migrations_table()
        SearchLight.Generator.new_table_migration(BookWithAuthor)
        SearchLight.Migration.up()
        SearchLight.Generator.new_table_migration(Author)
        SearchLight.Migration.up()

        testAuthor = Author(firstname="Johann Wolfgang", lastname="Goethe")
        testId = testAuthor |> save!

        @test length(find(Author)) > 0

        ####### tearDown #########
        SearchLight.Migration.all_down!!(confirm=false)
        tearDown(conn)
    end

    @safetestset "Saving and Reading Models with fields containing submodels" begin
        using SearchLight
        using SearchLightOracle
        using Main.TestSetupTeardown
        using Main.TestModels

        conn = prepareDbConnection()
        SearchLight.Migration.create_migrations_table()
        SearchLight.Generator.new_table_migration(BookWithAuthor)
        SearchLight.Migration.up()
        SearchLight.Generator.new_table_migration(Author)
        SearchLight.Migration.up()

        #create an author
        testAuthor = Author(firstname="John", lastname="Grisham")
        #create books from the author above and bring it to them
        testAuthor.books = map(book -> BookWithAuthor(title=book), seedBook())

        testId = testAuthor |> save!

        idAuthor = testAuthor.id.value
        for book in testId.books
          @test book.id_author.value == idAuthor
        end

        result = find(Author)
        @test length(result) > 0

        @test !isempty(result[1].books)
        @test length(result[1].books) == length(seedBook())

        ####### tearDown #########
        SearchLight.Migration.all_down!!(confirm=false)
        tearDown(conn)
    end;

    @safetestset "functions findone_or_create, updateby_or_create etc" begin
      using SearchLight
      using SearchLightOracle
      using Main.TestSetupTeardown
      using Main.TestModels

      conn = prepareDbConnection()
      SearchLight.Migration.create_migrations_table()
      SearchLight.Generator.new_table_migration(BookWithAuthor)
      SearchLight.Migration.up()
      SearchLight.Generator.new_table_migration(Author)
      SearchLight.Migration.up()

      #create an author
      testAuthor = Author(firstname="John", lastname="Grisham")
      #create books from the author above and bring it to them
      testAuthor.books = map(book -> BookWithAuthor(title=book), seedBook())

      result = findone_or_create(typeof(testAuthor))

      @test result !== nothing
      @test result.first_name == ""
      @test result.last_name == ""
      @test isempty(result.books)

      ####### tearDown #########
      SearchLight.Migration.all_down!!(confirm=false)
      tearDown(conn)
    end;

    @safetestset "several individual functions"
      using SearchLight
      using SearchLightOracle
      using Main.TestSetupTeardown
      using Main.TestModels

      conn = prepareDbConnection()
      SearchLight.Migration.create_migrations_table()
      SearchLight.Generator.new_table_migration(BookWithAuthor)
      SearchLight.Migration.up()
      SearchLight.Generator.new_table_migration(Author)
      SearchLight.Migration.up()

      #prepare some data for testing
      authors = [Author(firstname= author[1], lastname = author[2]) for author in (split.(seedAuthor()))]
      tmpAuthors = deepcopy(authors)
      authors |> SearchLight.save
      for i in 1:5
        realTmp = deepcopy(tmpAuthors)
        realTmp |> SearchLight.save
      end

      @time find(Author)

      ### tearDown
      SearchLight.Migration.all_down!!(confirm=false)
      tearDown(conn)

    end

end


text:
# This file was generated, do not modify it. # hide
X = coerce(X, autotype(X, rules=(:discrete_to_continuous,)))

dtr_model = DecisionTreeRegressor()
dtr = machine(dtr_model, X, y)

fit!(dtr, rows=train)

ypred = predict(dtr, rows=test)
round(rms(ypred, y[test]), sigdigits=3)


text:
# Mittelman instances
include(joinpath(@__DIR__, "src/LP/clp.jl"))
include(joinpath(@__DIR__, "src/LP/cplex.jl"))
include(joinpath(@__DIR__, "src/LP/ecos.jl"))
include(joinpath(@__DIR__, "src/LP/glpk.jl"))
include(joinpath(@__DIR__, "src/LP/gurobi.jl"))
include(joinpath(@__DIR__, "src/LP/mosek.jl"))
include(joinpath(@__DIR__, "src/LP/tulip.jl"))

for finst in ["afiro.mps", "25fv47.mps", "pds-02.mps", "qap8.mps", "pilot.mps"]
    fname = joinpath(@__DIR__, "dat/netlib", finst)

    run_clp(fname)
    run_cplex(fname)
    run_ecos(fname)
    run_glpk(fname)
    run_gurobi(fname)
    run_mosek(fname)
    run_tulip(fname)
end

# Structured LP instances
include("src/RMP/cplex.jl")
include("src/RMP/gurobi.jl")
include("src/RMP/mosek.jl")
include("src/RMP/tulip.jl")
include("src/RMP/tulip_uba.jl")

for finst in ["DER_24_1024_10.mps"]
    frmp = joinpath(@__DIR__, "dat/rmp", finst)

    run_rmp_cplex(frmp)
    run_rmp_gurobi(frmp)
    run_rmp_mosek(frmp)
    run_rmp_tulip(frmp)
    run_rmp_tulip_uba(frmp)
end

# Extended precision
include("src/D64/tulipD64.jl")

tulip_d64(joinpath(@__DIR__, "dat/plato/brazil3.mps"), 1e-8)

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
"""
Self-contained wrapper of C library based on Source Extractor
(https://github.com/kbarbary/sep) for performing image background estimation
and segmentation (detection). This is used in initialization in Celeste.

In the future, this module could be replaced by a pure-Julia implementation.
The background could be made an explicit part of the model.
"""

module SEP

import Base: collect, broadcast!, broadcast, -, show

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Celeste.SEP not properly installed. Please run Pkg.build(\"Celeste\")")
end

# Julia Type -> numeric type code from sep.h
const SEPMaskType = Union{Bool, UInt8, Cint, Cfloat, Cdouble}
const PixType = Union{UInt8, Cint, Cfloat, Cdouble}
sep_typecode(::Type{Bool}) = Cint(11)
sep_typecode(::Type{UInt8}) = Cint(11)
sep_typecode(::Type{Cint}) = Cint(31)
sep_typecode(::Type{Cfloat}) = Cint(42)
sep_typecode(::Type{Cdouble}) = Cint(82)

# definitions from sep.h
const SEP_NOISE_NONE = Cshort(0)
const SEP_NOISE_STDDEV = Cshort(1)
const SEP_NOISE_VAR = Cshort(2)
const SEP_THRESH_REL = Cint(0)
const SEP_THRESH_ABS = Cint(1)
const SEP_FILTER_CONV = Cint(0)
const SEP_FILTER_MATCHED = Cint(1)

function sep_assert_ok(status::Cint)
    if status != 0
        msg = Vector{UInt8}(61)
        ccall((:sep_get_errmsg, libsep), Void, (Int32,Ptr{UInt8}), status, msg)
        msg[end] = 0  # ensure NULL-termination, just in case.
        error(unsafe_string(pointer(msg)))
    end
end

# internal use only: mirrors `sep_image` struct in sep.h
struct sep_image
data::Ptr{Void}
noise::Ptr{Void}
mask::Ptr{Void}
dtype::Cint
ndtype::Cint
mdtype::Cint
w::Cint
h::Cint
noiseval::Cdouble
noise_type::Cshort
gain::Cdouble
maskthresh::Cdouble
end

function sep_image(data::Array{T, 2};
                   noise=nothing, mask=nothing, noise_type=:stddev,
                   gain=0.0, mask_thresh=0.0) where {T<:PixType}
    sz = size(data)

    # data is required
    data_ptr = Ptr{Void}(pointer(data))
    dtype = sep_typecode(T)

    # mask options
    mask_ptr = C_NULL
    mdtype = Cint(0)
    if mask !== nothing
        isa(mask, Matrix) || error("mask must be a 2-d array")
        size(mask) == sz || error("data and mask must be same size")
        mask_ptr = Ptr{Void}(pointer(mask))
        mdtype = sep_typecode(eltype(mask))
    end

    # noise options
    ndtype = Cint(0)
    noise_ptr = C_NULL
    noiseval = 0.0
    noise_typecode = SEP_NOISE_NONE
    if noise !== nothing
        if isa(noise, Matrix)
            size(noise) == sz || error("noise array must be same size as data")
            noise_ptr = Ptr{Void}(pointer(noise))
            ndtype = sep_typecode(eltype(noise))
        elseif isa(noise, Real)
            noiseval = Cdouble(noise)
        else
            error("noise must be array or number")
        end
        noise_typecode = ((noise_type == :stddev) ? SEP_NOISE_STDDEV :
                          (noise_type == :var) ? SEP_NOISE_VAR :
                          error("noise_type must be :stddev or :var"))
    end

    return sep_image(data_ptr, noise_ptr, mask_ptr,
                     dtype, ndtype, mdtype,
                     sz[1], sz[2],
                     noiseval, noise_typecode,
                     Cdouble(gain),
                     Cdouble(mask_thresh))
end

# ---------------------------------------------------------------------------
# Background functions

"""
    Background(data::Matrix; <keyword arguments>)

Spline representation of the variable background and noise of an image.

# Arguments

- `data::Matrix`:
- `mask=nothing`:
- `boxsize=(64, 64)`:
- `filtersize=(3, 3)`:
- `filterthresh=0.0`:
- `mask_thresh=0`:
"""
mutable struct Background
    ptr::Ptr{Void}
    data_size::Tuple{Int, Int}

    function Background(data::Array{T, 2} where T<:PixType;
                        mask=nothing, boxsize=(64, 64), filtersize=(3, 3),
                        filterthresh=0.0, mask_thresh=0)
        im = sep_image(data; mask=mask, mask_thresh=mask_thresh)
        result = Ref{Ptr{Void}}(C_NULL)
        status = ccall((:sep_background, libsep), Cint,
                       (Ptr{sep_image}, Cint, Cint, Cint, Cint, Cdouble,
                        Ref{Ptr{Void}}),
                       &im, boxsize[1], boxsize[2], filtersize[1],
                       filtersize[2], filterthresh, result)
        sep_assert_ok(status)
        bkg = new(result[], size(data))
        finalizer(bkg, free!)
        return bkg
    end
end

free!(bkg::Background) = ccall((:sep_bkg_free, libsep), Void, (Ptr{Void},),
                               bkg.ptr)

global_mean(bkg::Background) = ccall((:sep_bkg_global, libsep), Cfloat,
                                     (Ptr{Void},), bkg.ptr)

global_rms(bkg::Background) = ccall((:sep_bkg_globalrms, libsep), Cfloat,
                                    (Ptr{Void},), bkg.ptr)

function show(io::IO, bkg::Background)
    print(io, "Background $(bkg.data_size[1])×$(bkg.data_size[2])\n")
    print(io, " - global mean: $(global_mean(bkg))\n")
    print(io, " - global rms : $(global_rms(bkg))\n")
end

function collect(::Type{T}, bkg::Background) where {T<:PixType}
    result = Array{T}(bkg.data_size)
    status = ccall((:sep_bkg_array, libsep), Cint, (Ptr{Void}, Ptr{Void}, Cint),
                   bkg.ptr, result, sep_typecode(T))
    sep_assert_ok(status)
    return result
end
# default collection type is Float32, because that's what's natively stored in
# background.
# TODO: make default the input array type in `background` instead?
collect(bkg::Background) = collect(Float32, bkg)

"""
    rms(bkg)
    rms(T, bkg)

Return an array of the standard deviation of the background. The result is
the size of the original image and is of type `T`, if given.
"""
function rms(::Type{T}, bkg::Background) where {T<:PixType}
    result = Array{T}(bkg.data_size)
    status = ccall((:sep_bkg_rmsarray, libsep), Cint,
                   (Ptr{Void}, Ptr{Void}, Cint),
                   bkg.ptr, result, sep_typecode(T))
    sep_assert_ok(status)
    return result
end
rms(bkg::Background) = rms(Float32, bkg)

# In-place background subtraction: A .-= bkg
function broadcast!(-, A::Array{T, 2},  ::Array{T, 2}, bkg::Background) where {T<:PixType}
    if size(A) != bkg.data_size
         throw(DimensionMismatch("dimensions must match"))
    end
    status = ccall((:sep_bkg_subarray, libsep), Cint,
                   (Ptr{Void}, Ptr{Void}, Cint),
                   bkg.ptr, A, sep_typecode(T))
    sep_assert_ok(status)
end
function broadcast(-, A::Array{T, 2} where T<:PixType, bkg::Background)
    B = copy(A)
    B .-= bkg
    return B
end

function (-)(A::Array{T, 2}, bkg::Background) where {T<:PixType}
    B = copy(A)
    B .-= bkg
    return B
end

# -----------------------------------------------------------------------------
# Source Extraction

# Internal use only: Mirror of C struct
struct sep_catalog
    nobj::Cint                 # number of objects (length of all arrays)
    thresh::Ptr{Cfloat}              # threshold (ADU)
    npix::Ptr{Cint}              # # pixels extracted (size of pix array)
    tnpix::Ptr{Cint}                # # pixels above thresh (unconvolved)
    xmin::Ptr{Cint}
    xmax::Ptr{Cint}
    ymin::Ptr{Cint}
    ymax::Ptr{Cint}
    x::Ptr{Cdouble}              # barycenter (first moments)
    y::Ptr{Cdouble}
    x2::Ptr{Cdouble}     # second moments
    y2::Ptr{Cdouble}
    xy::Ptr{Cdouble}
    errx2::Ptr{Cdouble}  # second moment errors
    erry2::Ptr{Cdouble}
    errxy::Ptr{Cdouble}
    a::Ptr{Cfloat}       # ellipse parameters
    b::Ptr{Cfloat}
    theta::Ptr{Cfloat}
    cxx::Ptr{Cfloat}     # alternative ellipse parameters
    cyy::Ptr{Cfloat}
    cxy::Ptr{Cfloat}
    cflux::Ptr{Cfloat}   # total flux of pixels (convolved)
    flux::Ptr{Cfloat}    # total flux of pixels (unconvolved)
    cpeak::Ptr{Cfloat}   # peak pixel flux (convolved)
    peak::Ptr{Cfloat}    # peak pixel flux (unconvolved)
    xcpeak::Ptr{Cint}    # x, y coords of peak (convolved) pixel
    ycpeak::Ptr{Cint}
    xpeak::Ptr{Cint}     # x, y coords of peak (unconvolved) pixel
    ypeak::Ptr{Cint}
    flag::Ptr{Cshort}    # extraction flags
    pix::Ptr{Ptr{Cint}}  # array giving indicies of object's pixels in
                         # image (linearly indexed). Length is `npix`.
                         # (pointer to within the `objectspix` buffer)
    objectspix::Ptr{Cint}  # buffer holding pixel indicies for all objects
end

struct Catalog
npix::Vector{Cint}
xmin::Vector{Cint}
xmax::Vector{Cint}
ymin::Vector{Cint}
ymax::Vector{Cint}
x::Vector{Cdouble}
y::Vector{Cdouble}
a::Vector{Cfloat}
b::Vector{Cfloat}
theta::Vector{Cfloat}
cxx::Vector{Cfloat}
cyy::Vector{Cfloat}
cxy::Vector{Cfloat}
flux::Vector{Cfloat}
pix::Vector{Vector{Cint}}
end

function show(io::IO, cat::Catalog)
    print(io, "SEP.Catalog with $(length(cat.x)) entries")
end

"""
    unsafe_copy(src::Ptr{T}, N)

Create a `Vector{T}` and copy `N` elements from `Ptr{T}` into it.
"""
function unsafe_copy(src::Ptr{T}, N) where {T}
    dest = Vector{T}(N)
    unsafe_copy!(pointer(dest), src, N)
    return dest
end

"""
    extract(data, thresh; <keyword arguments>)

Perform image segmentation on the data and return a catalog of sources.

# Arguments

- `data::Matrix`: Image data.
- `thresh::Real`: Threshold for detection in absolute units, or in standard
  deviations if noise is given.
- `noise=nothing`:
- `noise_type=:stddev`:
- `mask=nothing`:
- `mask_thresh=0`:
- `minarea=5`:
- `filter_kernel`:
- `filter_type=:conv`:
- `deblend_nthresh=32`:
- `deblend_cont=0.005`: Minimum contrast in deblending. Set to 1.0 to disable
  deblending.
- `clean=true`: Perform cleaning?
- `clean_param=1.0`
- `gain=0.0`
"""
function extract(data::Array{T, 2} where T, thresh::Real;
                 noise=nothing, mask=nothing, noise_type=:stddev,
                 minarea=5,
                 filter_kernel=Float32[1 2 1; 2 4 2; 1 2 1],
                 filter_type=:convolution,
                 deblend_nthresh=32, deblend_cont=0.005,
                 clean=true, clean_param=1.0,
                 gain=0.0, mask_thresh=0)

    im = sep_image(data; noise=noise, mask=mask, noise_type=noise_type,
                   gain=gain, mask_thresh=mask_thresh)

    thresh_typecode = noise === nothing ? SEP_THRESH_ABS : SEP_THRESH_REL

    filter_typecode = ((filter_type == :matched)? SEP_FILTER_MATCHED :
                       (filter_type == :convolution)? SEP_FILTER_CONV :
                       error("filter_type must be :matched or :convolution"))

    # convert filter kernel to Cfloat array
    filter_kernel_cfloat = convert(Array{Cfloat, 2}, filter_kernel)
    filter_size = size(filter_kernel_cfloat)

    ccatalog_ptr_ref = Ref{Ptr{sep_catalog}}(C_NULL)
    status = ccall((:sep_extract, libsep), Cint,
                   (Ptr{sep_image},
                    Cfloat, Cint, # thresh, thresh_type
                    Cint,  # minarea
                    Ptr{Cfloat}, Cint, Cint,  # conv, convw, convh
                    Cint, # filter_type
                    Cint, Cdouble,  # deblend_nthresh, deblend_cont
                    Cint, Cdouble,  # clean_flag, clean_param
                    Ref{Ptr{sep_catalog}}),
                   &im,
                   thresh, thresh_typecode,
                   minarea,
                   filter_kernel_cfloat, filter_size[1], filter_size[2],
                   filter_typecode,
                   deblend_nthresh, deblend_cont,
                   clean, clean_param,
                   ccatalog_ptr_ref)
    sep_assert_ok(status)

    # get pointer to C-allocated result
    ccatalog_ptr = ccatalog_ptr_ref[]

    # copy result arrays into jula-managed memory
    ccatalog = unsafe_load(ccatalog_ptr)
    nobj = ccatalog.nobj
    # translate pixel index vectors
    npix = unsafe_copy(ccatalog.npix, nobj)
    pix = Vector{Vector{Cint}}(nobj)
    for i in eachindex(pix)
        ptr = unsafe_load(ccatalog.pix, i)  # ptr to pixel indicies for i-th obj
        pix[i] = unsafe_copy(ptr, npix[i])
        pix[i] .+= Cint(1) # change to 1-indexing
    end

    result = Catalog(npix,
                     unsafe_copy(ccatalog.xmin, nobj),
                     unsafe_copy(ccatalog.xmax, nobj),
                     unsafe_copy(ccatalog.ymin, nobj),
                     unsafe_copy(ccatalog.ymax, nobj),
                     unsafe_copy(ccatalog.x, nobj),
                     unsafe_copy(ccatalog.y, nobj),
                     unsafe_copy(ccatalog.a, nobj),
                     unsafe_copy(ccatalog.b, nobj),
                     unsafe_copy(ccatalog.theta, nobj),
                     unsafe_copy(ccatalog.cxx, nobj),
                     unsafe_copy(ccatalog.cyy, nobj),
                     unsafe_copy(ccatalog.cxy, nobj),
                     unsafe_copy(ccatalog.flux, nobj),
                     pix)

    # switch to 1 indexing
    result.x .+= 1.0
    result.y .+= 1.0

    # free result allocated in sep_extract
    ccall((:sep_catalog_free, libsep), Void, (Ptr{sep_catalog},), ccatalog_ptr)

    return result
end

# mask ellipse
"""
    mask_ellipse(data, x, y, cxx, cyy, cxy; r=1.0)

Set values within ellipse to `true`.
"""
function mask_ellipse(data::Array{T, 2} where T <: Union{UInt8, Bool},
                      x::Real, y::Real, cxx::Real, cyy::Real, cxy::Real;
                      r=1.0)
    w, h = size(data)
    ccall((:sep_set_ellipse, libsep), Void,
          (Ptr{Cuchar}, Cint, Cint, Cdouble, Cdouble, Cdouble, Cdouble,
           Cdouble, Cdouble, Cuchar), data, w, h, x-1.0, y-1.0, cxx, cyy,
          cxy, r, true)
end

end # module