text:
module VectorHeatEquationTests

using Gridap
using ForwardDiff
using LinearAlgebra
using Test
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Gridap.FESpaces: get_algebraic_operator

# using GridapODEs.ODETools: ThetaMethodLinear

import Gridap: ∇
import GridapODEs.TransientFETools: ∂t

θ = 1.0

u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t
u(t::Real) = x -> u(x,t)
v(x) = t -> u(x,t)
f(t) = x -> ∂t(u)(t)(x)-Δ(u(t))(x)

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 2

reffe = ReferenceFE(lagrangian,Float64,order)
V0 = FESpace(
  model,
  reffe,
  conformity=:H1,
  dirichlet_tags="boundary"
)
U = TransientTrialFESpace(V0,u)

Ω = Triangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)

#
a(u,v) = ∫(∇(v)⊙∇(u))dΩ
b(v,t) = ∫(v⋅f(t))dΩ
m(ut,v) = ∫(ut⋅v)dΩ

X = TransientMultiFieldFESpace([U,U])
Y = MultiFieldFESpace([V0,V0])

_res(t,u,ut,v) = a(u,v) + m(ut,v) - b(v,t)

res(t,(u1,u2),(u1t,u2t),(v1,v2)) = _res(t,u1,u1t,v1) + _res(t,u2,u2t,v2)
jac(t,x,xt,(du1,du2),(v1,v2)) = a(du1,v1) + a(du2,v2)
jac_t(t,x,xt,(du1t,du2t),(v1,v2)) = m(du1t,v1) + m(du2t,v2)

op = TransientFEOperator(res,jac,jac_t,X,Y)

t0 = 0.0
tF = 1.0
dt = 0.1

U0 = U(0.0)
X0 = X(0.0)
uh0 = interpolate_everywhere(u(0.0),U0)
xh0 = interpolate_everywhere([uh0,uh0],X0)

ls = LUSolver()
# using Gridap.Algebra: NewtonRaphsonSolver
# nls = NLSolver(ls;show_trace=true,method=:newton) #linesearch=BackTracking())
odes = ThetaMethod(ls,dt,θ)
solver = TransientFESolver(odes)

sol_t = solve(solver,op,xh0,t0,tF)

l2(w) = w⋅w

tol = 1.0e-6
_t_n = t0

result = Base.iterate(sol_t)

for (xh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  uh_tn = xh_tn[1]
  e = u(tn) - uh_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
  @test el2 < tol
end

end #module


text:
module MSM

    #---------------------------------------------------------------------------
    # Dependencies
    #---------------------------------------------------------------------------
    using BlackBoxOptim
    using Optim
    using BSON
    using CSV
    using FiniteDifferences
    using DataFrames
    using DataStructures
    using OrderedCollections
    using Dates
    using Distributed
    using Random
    using Logging
    using Statistics
    using Distributions
    using LinearAlgebra
    using SharedArrays
    using LatinHypercubeSampling
    using Sobol
    using StatsBase

    # Exports from BlackBoxOptim
    #---------------------------
    export best_candidate

    #---------------------------------------------------------------------------
    # Includes
    #---------------------------------------------------------------------------
    # Types
    #------
    include("types.jl");

    # API
    #----
    include("api.jl")

    # General functions, useful at several places
    #---------------------------------------------
    include("generic.jl")

    # Functions to load and save
    #--------------------------
    include("save_load.jl")

    # Functions to minimize the objective function
    #---------------------------------------------
    include("optimize.jl")

    # Functions to do inference
    #---------------------------------------------
    include("econometrics.jl")

    # Functions to do plots
    #---------------------------------------------
    include("analysis.jl")

    # Exports
    #--------
    # Functions and types in types.jl
    #----------------------------------
    export MSMOptions, MSMProblem
    export default_function, rosenbrock2d
    export is_global_optimizer, is_local_optimizer
    export convert_to_optim_algo, convert_to_fminbox
    export is_bb_optimizer, is_optim_optimizer

    # Functions and types in api.jl
    #-------------------------------

    # Functions and types in generic.jl
    #----------------------------------
    export set_simulate_empirical_moments!, construct_objective_function!
    export set_priors!, set_weight_matrix!, set_empirical_moments!, set_Sigma0!
    export set_bbSetup!, generate_bbSearchRange
    export create_lower_bound, create_upper_bound
    export set_global_optimizer!
    export latin_hypercube_sampling, sobol_sampling
    export get_now, info, linspace

    # Functions and types in save_load.jl
    #------------------------------------
    export read_priors, read_empirical_moments

    #export saveMSMOptim, loadMSMOptim

    # Functions and types in optimize.jl
    #-----------------------------------
    export msm_optimize!
    export msm_minimizer, msm_minimum
    export msm_refine_globalmin!
    export msm_local_minimizer, msm_local_minimum
    export msm_localmin, msm_multistart!
    export msm_multistart_minimizer, msm_multistart_minimum

    # Functions in econometrics.jl
    #-----------------------------
    export calculate_D, calculate_Avar!, calculate_se, calculate_t, calculate_pvalue, calculate_CI
    export summary_table, cov_NW, J_test

    # Functions in analysis.jl
    #-------------------------
    export msm_slices

end

ext: jl


text:
"""
function get_EU_LV_feeder_data()::Dict{String,Any}
    Gets the MATHEMATICAL data dictionary for the European Low Voltage Test Feeder.
"""
get_EU_LV_feeder_data(;power_base::Float64=1.0) = get_feeder_data(1, 1; power_base = power_base)

function get_feeder_data(ntw::Int64, fdr::Int64; power_base::Float64=1.0)::Dict{String,Any}
    #settings that can be changed

    rm_transfo = true
    rd_lines = true
    season = "winter"
    time_step = 144
    elm = ["load"]
    pfs = [0.95]

    #"fixed" part below
    data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr),data_model=_PMD.ENGINEERING)
    if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
    if rd_lines _PMDSE.reduce_enwl_lines_eng!(data) end
    data["settings"]["sbase_default"] = power_base

    _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

    return _PMD.transform_data_model(data)
end

max_stars_repo_name: MartaVanin/PaperData

max_issues_repo_name: MartaVanin/PaperData

max_forks_repo_name: MartaVanin/PaperData


text:
module TestForward
using CompressedSensing
using CompressedSensing: sparse_data, fr, foba, perturb, FR

using Test
using LinearAlgebra
n, m, k = 32, 48, 3
A, x, b = sparse_data(n = n, m = m, k = k) # needs to be determined
δ = 1e-2
y = perturb(b, δ)

@testset "forward regression" begin
    # noiseless
    xfr = fr(A, b, sparsity = k)
    @test xfr.nzind == x.nzind
    @test xfr.nzval ≈ x.nzval

    # with small perturbation
    xfr = fr(A, y, sparsity = k)
    @test xfr.nzind == x.nzind # support recovery
    @test isapprox(xfr.nzval, x.nzval, atol = 2δ)

    # test initialization with support set
    nzind = [1, 2, 3]
    P = FR(A, y, nzind)
    xi = P.AiQR \ y
    @test length(xi) == 3
    @test xi ≈ P.A[:, nzind] \ y
end

@testset "FoBa" begin
    xfr = foba(A, y, δ)
    @test xfr.nzind == x.nzind # support recovery
    @test isapprox(xfr.nzval, x.nzval, atol = 2δ)
end

end # TestForward


text: #Phasing is tested in test_data_processing.jl


text:
__precompile__()
"""
This is the main module and provides interfaces for specifying and solving reachability problems.
"""
module Reachability

using Reexport, RecipesBase, Memento, MathematicalSystems,
      MathematicalPredicates, HybridSystems, LinearAlgebra, Suppressor,
      SparseArrays, Printf

@reexport using LazySets
import LazySets: use_precise_ρ
using LazySets: LinearMap, AffineMap, ResetMap, Interval

import LazySets.Approximations: project

include("logging.jl")
include("Options/dictionary.jl")
include("Options/validation.jl")
include("Options/default_options.jl")
include("Utils/Utils.jl")
include("ReachSets/ReachSets.jl")

@reexport using Reachability.Utils,
                Reachability.ReachSets

include("solve.jl")
include("plot_recipes.jl")

# specify behavior of line search algorithm by dispatching on post operator
global discrete_post_operator = ConcreteDiscretePost()

@suppress_err begin
    function use_precise_ρ(cap::Intersection{N})::Bool where N<:Real
        return use_precise_ρ(discrete_post_operator, cap)
    end
end

end # module


text:
# Gun: variant S

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test

include(joinpath("..", "rk_helper", "gun_test_utils.jl"))

@bench @testset "NLEIGS: Gun variant S" begin
    verbose = displaylevel

    nep, Σ, Ξ, v, nodes, funres = gun_init()

    # solve nlep
    lambda, X, res, solution_info = nleigs(nep, Σ, Ξ=Ξ, logger=verbose > 0 ? 1 : 0, minit=70, maxit=100, v=v, nodes=nodes, static=true, errmeasure=funres, return_details=verbose > 1)

    verify_lambdas(21, nep, lambda, X)

    if verbose > 1
        include("nleigs_residual_plot.jl")
        nleigs_residual_plot("Gun: variant S", solution_info, Σ; ylims=[1e-17, 1e-1])
    end
end


text:
"Exception thrown when an error occurred calling MXNet API."
immutable MXError <: Exception
  msg :: AbstractString
end

################################################################################
# Common types used in MXNet API
################################################################################
typealias MX_uint Cuint
typealias MX_float Cfloat
typealias MX_handle Ptr{Void}

typealias char_p Ptr{UInt8}
typealias char_pp Ptr{char_p}

################################################################################
# Initialization and library API entrance
################################################################################
const MXNET_LIB = Libdl.find_library(["libmxnet.so","libmxnet.dll"],
                                     [joinpath("$(get(ENV,"MXNET_HOME",""))","lib"),
                                      joinpath(Pkg.dir("MXNet"),"deps/usr/lib")])
if isempty(MXNET_LIB)
  # touch this file, so that after the user properly build libmxnet, the precompiled
  # MXNet.ji will be re-compiled to get MXNET_LIB properly.
  touch(@__FILE__)
  error("Cannot find or load libmxnet.so. Please see the document on how to build it.")
else
  include_dependency(MXNET_LIB)
end

function __init__()
  _populate_symbol_creator_cache!()
  _populate_iter_creator_cache!()

  atexit() do
    # notify libmxnet we are shutting down
    ccall( ("MXNotifyShutdown", MXNET_LIB), Cint, () )
  end
end

function mx_get_last_error()
  msg = ccall( ("MXGetLastError", MXNET_LIB), char_p, () )
  if msg == C_NULL
    throw(MXError("Failed to get last error message"))
  end
  return unsafe_string(msg)
end

"Utility macro to call MXNet API functions"
macro mxcall(fv, argtypes, args...)
  f = eval(fv)
  args = map(esc, args)
  quote
    _mxret = ccall( ($(Meta.quot(f)), $MXNET_LIB),
                    Cint, $argtypes, $(args...) )
    if _mxret != 0
      err_msg = mx_get_last_error()
      throw(MXError(err_msg))
    end
  end
end

################################################################################
# Handle types
################################################################################
macro mx_define_handle_t(name, destructor)
  name = esc(name)
  quote
    type $name
      value :: MX_handle

      function $name(value = C_NULL)
        hdr = new(value)

        $(if destructor != :nop
          :(finalizer(hdr, delete!))
        end)

        return hdr
      end
    end

    $(if finalizer != :nop
      quote
        function delete!(h :: $name)
          if h.value != C_NULL
            @mxcall($(Meta.quot(destructor)), (MX_handle,), h.value)
            h.value = C_NULL
          end
        end
      end
    end)

    function Base.unsafe_convert(::Type{MX_handle}, obj::$name)
      obj.value
    end
    Base.convert(t::Type{MX_handle}, obj::$name) = Base.unsafe_convert(t, obj)
    Base.cconvert(t::Type{MX_handle}, obj::$name) = Base.unsafe_convert(t, obj)

    function Base.isnull(obj::$name) obj.value == C_NULL end
  end
end

@mx_define_handle_t(MX_NDArrayHandle, MXNDArrayFree)
@mx_define_handle_t(MX_FunctionHandle, nop)
@mx_define_handle_t(MX_SymbolHandle, MXSymbolFree)
@mx_define_handle_t(MX_ExecutorHandle, MXExecutorFree)
@mx_define_handle_t(MX_DataIterHandle, MXDataIterFree)
@mx_define_handle_t(MX_KVStoreHandle, MXKVStoreFree)

################################################################################
# MXNet Params
#
# MXNet API use string to pass some common parameters like the configurations
# when defining layers. Typically, it is enough to use string(obj) to get a
# recognizable representation for libmxnet. However, there is currently a
# caveat:
#
# Because Julia use column-major ordering for tensors. In order to properly
# interact with Julia Arrays, the shape will look "reversed" from the Julia
# side. For example, a typical MNIST mini-batch tensor is of shape (28,28,1,100)
# from Julia side, while the shape information for the same piece of memory
# should be interpreted as (100,1,28,28) from C/C++/Python side.
#
# Therefore, when passing parameters to libmxnet, we should reverse the shape
# parameter. For example, when the user specify a non-square kernel size for
# a convolution or pooling layer. Unfortunately, those operators are automatically
# imported, and information about the type of each parameter is somehow limited.
# One hacky way is to match the type description for the string "Shape(tuple)"
# when importing operators. But currently we simply decided to reverse **all**
# NTuple{N, Int} passed to libmxnet.
#
# TODO: find a better solution in case this cause issues in the future.
################################################################################
function dump_mx_param(val :: Any)
  string(val)
end
function dump_mx_param{N,T<:Integer}(shape :: NTuple{N, T})
  string(tuple(flipdim([shape...],1)...))
end

"""A convenient macro copied from Mocha.jl that could be used to define structs
with default values and type checks. For example
```julia
@defstruct MyStruct Any (
  field1 :: Int = 0,
  (field2 :: AbstractString = "", !isempty(field2))
)
```
where each field could be either
```julia
field_name :: field_type = default_value
```
or put within a tuple, with the second element
specifying a validation check on the field value.
In the example above, the default value for
field2 does not satisfy the assertion, this
could be used to force user to provide a
valid value when no meaningful default value
is available.

The macro will define a constructor that could accept
the keyword arguments.
"""
macro defstruct(name, fields)
  _defstruct_impl(false, name, fields)
end

"""A convenient macro to define immutable structs. The same as
`@defstruct` except that the defined type is immutable.
"""
macro defimmutable(name, fields)
  _defstruct_impl(true, name, fields)
end

"""Internal use only, this value is used to indicate a required value
is not specified.
"""
immutable __Undefined
end

function _defstruct_impl(is_immutable, name, fields)
  if isa(fields, Expr) && fields.head == :tuple
    fields = fields.args
  else
    fields = [fields]
  end
  @assert length(fields) > 0

  if isa(name, Symbol)
    name       = esc(name)
    super_name = :Any
  elseif VERSION >= v"0.5-"
    @assert(isa(name, Expr) && name.head == :(<:) && length(name.args) == 2 &&
            isa(name.args[1], Symbol) && isa(name.args[2], Symbol),
            "name must be of form 'Name <: SuperType'")

    super_name = esc(name.args[2])
    name       = esc(name.args[1])
  else
    @assert(isa(name, Expr) && name.head == :comparison &&
            length(name.args) == 3 && name.args[2] == :(<:) &&
            isa(name.args[1], Symbol) && isa(name.args[3], Symbol),
            "name must be of form 'Name <: SuperType'")

    super_name = esc(name.args[3])
    name       = esc(name.args[1])
  end

  field_defs     = Array(Expr, length(fields))        # :(field2 :: Int)
  field_names    = Array(Expr, length(fields))        # :field2
  field_defaults = Array(Expr, length(fields))        # :(field2 = 0)
  field_types    = Array(Expr, length(fields))        # Int
  field_asserts  = Array(Expr, length(fields))        # :(field2 >= 0)
  required_field = Symbol[]

  for i = 1:length(fields)
    field = fields[i]
    if field.head == :tuple
      field_asserts[i] = esc(field.args[2])
      field = field.args[1]
    end
    if field.head == :(=)
      fname             = field.args[1].args[1]
      field_defs[i]     = esc(field.args[1])
      field_names[i]    = esc(fname)
      field_types[i]    = esc(field.args[1].args[2])
      field_defaults[i] = Expr(:kw, fname, esc(field.args[2]))
    else
      # no default value provided, required field
      fname             = field.args[1]
      field_defs[i]     = esc(field)
      field_names[i]    = esc(fname)
      field_types[i]    = esc(field.args[2])
      field_defaults[i] = Expr(:kw, fname, __Undefined())
      push!(required_field, fname)
    end
  end

  # body of layer type, defining fields
  type_body = Expr(:block, field_defs...)

  # constructor
  requires = map(required_field) do fname
    :(@assert(!isa($fname, __Undefined), "value for " * string($fname) * " is required"))
  end
  converts = map(zip(field_names, field_types)) do param
    f_name, f_type = param
    :($f_name = convert($f_type, $f_name))
  end
  asserts = map(filter(i -> isdefined(field_asserts,i), 1:length(fields))) do i
    :(@assert($(field_asserts[i])))
  end
  construct = Expr(:call, name, field_names...)
  ctor_body = Expr(:block, requires..., converts..., asserts..., construct)
  ctor_def = Expr(:call, name, Expr(:parameters, field_defaults...))
  ctor = Expr(:(=), ctor_def, ctor_body)

  if is_immutable
    quote
      immutable $(name) <: $(super_name)
        $type_body
      end

      $ctor
    end
  else
    quote
      type $(name) <: $(super_name)
        $type_body
      end

      $ctor
    end
  end
end


text:
# script to run inference for the OU SDEMEM model

#using Pkg
using PyPlot
using LinearAlgebra
using KernelDensity
using DataFrames
using KernelDensity
using Distributions

# load functions
include(pwd()*"/src/SDEMEM OU process/ou_sdemem.jl")
include(pwd()*"/src/SDEMEM OU process/mcmc.jl")

N_time = 200
M_subjects = 40

nbr_particles_pf = 5
nbr_particles_cpf = 5

y,x,t_vec,dt,η,σ_ϵ,ϕ,prior_parameters_η,prior_parameters_σ_ϵ = set_up(M=M_subjects,N=N_time,seed=100)

PyPlot.figure()
PyPlot.plot(Array(x'))

PyPlot.figure()
PyPlot.plot(Array(y'))

σ_ϵ = 0.31

@time loglik_kalman, path_kalman = kalman(y, σ_ϵ, ϕ, dt, true)

@time loglik_pf, path_pf = cpf(y, σ_ϵ, ϕ, dt, randn(nbr_particles_cpf,N_time+1,M_subjects), randn(N_time,2,M_subjects), nbr_particles_pf, false, true)

@time loglik_corr_pf, path_corr_pf = cpf(y, σ_ϵ, ϕ, dt, randn(nbr_particles_cpf,N_time+1,M_subjects), randn(N_time,2,M_subjects), nbr_particles_cpf, true, true)

sum(loglik_kalman)

sum(loglik_pf)

sum(loglik_corr_pf)

subject = rand(1:M_subjects)

PyPlot.figure()
PyPlot.plot(path_kalman[subject], "r")
PyPlot.plot(y[subject,:], "g")
display(gcf())

PyPlot.figure()
PyPlot.plot(path_pf[subject][:,:]', "r")
PyPlot.plot(y[subject,:], "g")
display(gcf())

PyPlot.figure()
PyPlot.plot(path_corr_pf[subject][:,:]', "r")
PyPlot.plot(y[subject,:], "g")
display(gcf())

loglik_pf_m = zeros(M_subjects,100)
loglik_cpf_m = zeros(M_subjects,100)

for i in 1:100
    loglik_pf_m[:,i] = cpf(y, σ_ϵ, ϕ, dt, randn(nbr_particles_pf,N_time+1,M_subjects), randn(N_time,2,M_subjects), nbr_particles_pf, false)
    loglik_cpf_m[:,i] = cpf(y, σ_ϵ, ϕ, dt, randn(nbr_particles_cpf,N_time+1,M_subjects), randn(N_time,2,M_subjects), nbr_particles_cpf, true)
end

loglik_pf_m

loglik_cpf_m

loglik_pf_m = sum(loglik_pf_m, dims = 1)

loglik_cpf_m = sum(loglik_cpf_m, dims = 1)

std(loglik_pf_m)

std(loglik_cpf_m)

mean(loglik_pf_m,dims=2)

mean(loglik_cpf_m,dims=2)

#=

# old stuff test for different parameter values

loglik_pg_dist = kde(loglik_pf_m[:])
loglik_pg_dist_cpf = kde(loglik_cpf_m[:])

PyPlot.figure()
PyPlot.plot(loglik_pg_dist.x,loglik_pg_dist.density, "b")
PyPlot.plot(loglik_pg_dist_cpf.x,loglik_pg_dist_cpf.density, "--b")
PyPlot.plot((sum(loglik_kalman), sum(loglik_kalman)), (0, maximum(loglik_pg_dist.density)), "k")

# test for other parameter values
prior_dist_σ_ϵ = Gamma(prior_parameters_σ_ϵ[1],1/prior_parameters_σ_ϵ[2])

prior_dist_τ_1 = Gamma(prior_parameters_η[1,3],1/prior_parameters_η[1,4])
prior_dist_τ_2 = Gamma(prior_parameters_η[2,3],1/prior_parameters_η[2,4])
prior_dist_τ_3 = Gamma(prior_parameters_η[3,3],1/prior_parameters_η[3,4])

# check loglik est val
diff_kalman_pf = zeros(5)

PyPlot.figure()

for i in 1:5
    Random.seed!(i)

    σ_ϵ_new = rand(prior_dist_σ_ϵ)

    τ_1_new  = rand(prior_dist_τ_1)
    τ_2_new  = rand(prior_dist_τ_2)
    τ_3_new  = rand(prior_dist_τ_3)

    μ_1_new  = rand(Normal(prior_parameters_η[1,1], sqrt(1/(prior_parameters_η[1,2]*τ_1_new))))
    μ_2_new  = rand(Normal(prior_parameters_η[2,1], sqrt(1/(prior_parameters_η[2,2]*τ_2_new))))
    μ_3_new  = rand(Normal(prior_parameters_η[3,1], sqrt(1/(prior_parameters_η[3,2]*τ_3_new))))

    log_θ_1_new = μ_1_new .+ sqrt(1/τ_1_new)*randn(M_subjects)
    log_θ_2_new = μ_1_new .+ sqrt(1/τ_2_new)*randn(M_subjects)
    log_θ_3_new = μ_1_new .+ sqrt(1/τ_3_new)*randn(M_subjects)

    η_new = [μ_1_new; μ_2_new; μ_3_new; τ_1_new; τ_2_new; τ_3_new]
    ϕ_new = [log_θ_1_new log_θ_2_new log_θ_3_new]

    @time loglik_kalman, y_hat, x_hat, res_kalman = kalman(y, σ_ϵ_new, ϕ_new, dt, true)

    loglik_kalman = sum(loglik_kalman)

    loglik_pf_m = zeros(M_subjects,100)

    for i in 1:100
        loglik_pf_m[:,i] = cpf(y, σ_ϵ_new, ϕ_new, dt, randn(nbr_particles_pf,N_time+1,M_subjects), randn(N_time,2,M_subjects), nbr_particles_pf, false)
    end

    loglik_pf_m = sum(loglik_pf_m, dims = 1)

    loglik_pg_dist = kde(loglik_pf_m[:])

    loglik_cpf_m = zeros(M_subjects,100)

    for i in 1:100
        loglik_cpf_m[:,i] = cpf(y, σ_ϵ_new, ϕ_new, dt, randn(nbr_particles_cpf,N_time+1,M_subjects), randn(N_time,2,M_subjects), nbr_particles_cpf, true)
    end

    loglik_cpf_m = sum(loglik_cpf_m, dims = 1)

    loglik_pg_dist_cpf = kde(loglik_cpf_m[:])

    PyPlot.subplot(2,2,i)
    PyPlot.plot(loglik_pg_dist.x,loglik_pg_dist.density, "b")
    PyPlot.plot(loglik_pg_dist_cpf.x,loglik_pg_dist_cpf.density, "--b")
    PyPlot.plot((loglik_kalman, loglik_kalman), (0, maximum(loglik_pg_dist.density)), "k")

    diff_kalman_pf[i] = abs(loglik_kalman - mean(loglik_pf_m))
end

diff_kalman_pf

PyPlot.figure()
PyPlot.plot(diff_kalman_pf, "*")

# check loglik tracking

for i in 1:10
    PyPlot.figure()
    plot_counter = 0
    for idx in 1:M_subjects

        Random.seed!(i)
        σ_ϵ_new = rand(prior_dist_σ_ϵ)
        θ_1_new = rand(prior_dist_θ_1)
        θ_3_new = rand(prior_dist_θ_3)
        ϕ_new = rand(dist_θ_2, M_subjects)

        @time loglik_kalman, partial_loglik_sum_kalman = kalman(y, [θ_1_new; θ_3_new], θ_3_new, ϕ_new, dt, false, true)

        @time loglik_pf, w, x, partial_loglik_sum_pf = pf(y, [θ_1_new; θ_3_new], θ_3_new, ϕ_new, dt, 1000, true, true)

        #idx = 10
        partial_loglik_kalman = diff(partial_loglik_sum_kalman,dims=2)
        partial_loglik_pf = diff(partial_loglik_sum_pf,dims=2)

        plot_counter = plot_counter+1
        PyPlot.subplot(10,4,plot_counter)
        PyPlot.plot(partial_loglik_pf[idx,:], "r")
        PyPlot.plot(partial_loglik_kalman[idx,:], "b")

    end

end

# check path tracking

for i in 1:10
    PyPlot.figure()
    plot_counter = 0
    for idx in 1:M_subjects

        Random.seed!(i)
        σ_ϵ_new = rand(prior_dist_σ_ϵ)
        θ_1_new = rand(prior_dist_θ_1)
        θ_3_new = rand(prior_dist_θ_3)
        ϕ_new = rand(dist_θ_2, M_subjects)

        @time loglik_kalman, partial_loglik_sum_kalman = kalman(y, [θ_1_new; θ_3_new], θ_3_new, ϕ_new, dt, false, true)

        @time loglik_pf, w, x, partial_loglik_sum_pf = pf(y, [θ_1_new; θ_3_new], θ_3_new, ϕ_new, dt, 1000, true, true)

        plot_counter = plot_counter+1
        PyPlot.subplot(10,4,plot_counter)
        PyPlot.plot(Array(x[idx,1:10,:]'), "--")
        PyPlot.plot(Array(y[idx,:]), "b")

    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

function _check_rand_compat(s::Sampleable{Multivariate}, A::Union{AbstractVector,AbstractMatrix})
    size(A, 1) == length(s) || throw(DimensionMismatch("Output size inconsistent with sample length."))
    nothing
end

@doc doc"""
    bat_sampler(d::Distribution)

*BAT-internal, not part of stable public API.*

Tries to return a BAT-compatible sampler for Distribution d. A sampler is
BAT-compatible if it supports random number generation using an arbitrary
`AbstractRNG`:

    rand(rng::AbstractRNG, s::SamplerType)
    rand!(rng::AbstractRNG, s::SamplerType, x::AbstractArray)

If no specific method of `bat_sampler` is defined for the type of `d`, it will
default to `sampler(d)`, which may or may not return a BAT-compatible
sampler.
"""
function bat_sampler end

# ToDo/Decision: Rename? Replace with `bat_default` mechanism?
bat_sampler(d::Distribution) = Distributions.sampler(d)

@doc doc"""
    issymmetric_around_origin(d::Distribution)

*BAT-internal, not part of stable public API.*

Returns `true` (resp. `false`) if the Distribution is symmetric (resp.
non-symmetric) around the origin.
"""
function issymmetric_around_origin end

issymmetric_around_origin(d::Normal) = d.μ ≈ 0

issymmetric_around_origin(d::Gamma) = false

issymmetric_around_origin(d::Chisq) = false

issymmetric_around_origin(d::TDist) = true

issymmetric_around_origin(d::MvNormal) = iszero(d.μ)

issymmetric_around_origin(d::Distributions.GenericMvTDist) = iszero(d.μ)

const PosDefMatLike = Union{AbstractMatrix{T},AbstractPDMat{T},Cholesky{T}} where {T<:Real}

function cov2pdmat(::Type{T}, Σ::AbstractMatrix{<:Real}) where {T<:Real}
    Σ_conv = convert(Matrix{T}, Σ)
    PDMat(cholesky(Positive, Hermitian(Σ_conv)))
end

cov2pdmat(::Type{T}, Σ::PDMat{T}) where {T<:Real} = Σ
cov2pdmat(::Type{T}, Σ::AbstractPDMat) where {T<:Real} = cov2pdmat(T, Matrix(Σ))

cov2pdmat(::Type{T}, Σ::Cholesky{T}) where {T<:Real} = PDMat(Σ)
cov2pdmat(::Type{T}, Σ::Cholesky) where {T<:Real} = cov2pdmat(T, Matrix(Σ))

function get_cov end

get_cov(d::Distributions.GenericMvTDist) = d.Σ

function set_cov end

function set_cov(d::Distributions.GenericMvTDist{T,Cov}, Σ::PosDefMatLike) where {T,Cov<:PDMat{T}}
    Σ_conv = cov2pdmat(T, Σ)
    Distributions.GenericMvTDist(d.df, deepcopy(d.μ), Σ_conv)
end


text:
"""
Air properties.

* [`air_molar`](@ref)
* [`air_density`](@ref)
"""
module Air

include("docstring_style.jl")

using Canopy.Constants: R, M_w, M_d
using Canopy.Water: vapor_mole_frac

export air_molar,
    air_density

"""
Calculate the molar concentration [mol m^-3] of air from temperature [K] and
pressure [Pa].

# Examples

```jldoctest
julia> air_molar(298.15, 101325.0)
40.87405837767172
```
"""
air_molar(temp, pressure) = pressure / (R * temp)

"""
Calculate the air density [kg m^-3] from temperature [K], pressure [Pa], and
relative humidity [0--1].

# Examples

```jldoctest
julia> air_density(273.15, 101325.0, 0.0)
1.2922525730764916

julia> air_density(298.15, 101325.0, 1.0)
1.1699164154668222
```
"""
function air_density(temp, pressure, RH)
    chi_w = vapor_mole_frac(temp, pressure, RH)
    # molar mass of moist air [kg mol^-1]
    M_a = (1. - chi_w) * M_d + chi_w * M_w
    return M_a * air_molar(temp, pressure)
end

end  # module


text:
ENV["TRAIN"] = true
ENV["DEBUG"] = false
ENV["EVALUATE"] = false
include("notebook.jl")


text:
#!/usr/bin/env julia
using Lint
using LightXML

#=
    Code check prints xml checkStyle
=#

function main(args)

    xdoc = XMLDocument()
    xroot = create_root(xdoc, "checkstyle")

    for file in args
        file_root = new_child(xroot, "file")

        errs = collect(lintfile(file))
        if length(errs) > 0
            set_attribute(file_root, "name", file)
            for err in errs
                xs1 = new_child(file_root, "error")
                xml_info = Dict("line"=>err.line, "column" => "1",
                            "severity"=>"warning",
                            "message"=>err.scope*"(): "*err.variable*": "*err.message*"("*string(err.code)*")",
                            "source"=>"com.puppycrawl.tools.checkstyle.checks.AvoidEscapedUnicodeCharactersCheck")
                set_attributes(xs1, xml_info)
            end
        end
    end
    save_file(xdoc, "julia_lint-sonar.xml")

end

main(ARGS)

max_stars_repo_name: KarolisM/rnd-it-molbio-fidelity-umi

max_issues_repo_name: KarolisM/rnd-it-molbio-fidelity-umi

max_forks_repo_name: KarolisM/rnd-it-molbio-fidelity-umi


text:
using Turing: _hmc_step

include("unit_test_helper.jl")
include("gdemo.jl")

# Turing

mf = gdemo()
chn = sample(mf, HMC(0.05, 5), 2000)

println("mean of s: ", mean(chn[1000:end, :s]))
println("mean of m: ", mean(chn[1000:end, :m]))

# Plain Julia

std = ones(θ_dim)
θ = randn(θ_dim)
lj = lj_func(θ)

chn = Dict(:θ=>Vector{Vector{Float64}}(), :logϵ=>Vector{Float64}())
accept_num = 0

function dummy_print(args...)
  nothing
end

totla_num = 5000
for iter = 1:totla_num
  push!(chn[:θ], θ)
  θ, lj, is_accept, τ_valid, α = _hmc_step(θ, lj, lj_func, grad_func, 5, 0.05, std)
  accept_num += is_accept
end

@show lj
samples_s = exp.(map(x -> x[1], chn[:θ]))
samples_m = map(x -> x[2], chn[:θ])
@show mean(samples_s[1000:end])
@show mean(samples_m[1000:end])
@show std(samples_s[1000:end])
@show std(samples_m[1000:end])

@show accept_num / totla_num

# Unit tests


text:
const LeafOrNode{S, T} = Union{Leaf{T}, Node{S, T}}

struct TreeOOB{S, T}
    tree     :: LeafOrNode{S, T}
    inds     :: AbstractVector{Int}
    oobfeat  :: AbstractMatrix{S}
end

struct EnsembleOOB{S, T}
    nsamp    :: Int
    features :: AbstractMatrix{S}
    trees    :: Vector{TreeOOB{S, T}}
end

mk_rng(rng::Random.AbstractRNG) = rng
mk_rng(seed::T) where T <: Integer = Random.MersenneTwister(seed)


text:
# define new aggregator algorithms here and specify their properties

struct Direct <: AbstractAggregatorAlgorithm end
struct DirectFW <: AbstractAggregatorAlgorithm end
struct DirectCR <: AbstractAggregatorAlgorithm end
struct SortingDirect <: AbstractAggregatorAlgorithm end
struct RSSA <: AbstractAggregatorAlgorithm end
struct FRM <: AbstractAggregatorAlgorithm end
struct FRMFW <: AbstractAggregatorAlgorithm end
struct NRM <: AbstractAggregatorAlgorithm end

const JUMP_AGGREGATORS = (Direct(),DirectFW(),DirectCR(),SortingDirect(),RSSA(),FRM(),FRMFW(),NRM())

# For JumpProblem construction without an aggregator
struct NullAggregator <: AbstractAggregatorAlgorithm end

# true if aggregator requires a jump dependency graph
needs_depgraph(aggregator::AbstractAggregatorAlgorithm) = false
needs_depgraph(aggregator::DirectCR) = true
needs_depgraph(aggregator::SortingDirect) = true
needs_depgraph(aggregator::NRM) = true

# true if aggregator requires a map from solution variable to dependent jumps.
# It is implicitly assumed these aggregators also require the reverse map, from
# jumps to variables they update.
needs_vartojumps_map(aggregator::AbstractAggregatorAlgorithm) = false
needs_vartojumps_map(aggregator::RSSA) = true


text:
import Base.==

mutable struct Basis{B, V, P, T} <: AbstractBasis
    """The equations of the basis"""
    basis::B
    """Dependent (state) variables"""
    variables::V
    """Parameter variables"""
    parameter::P
    """Independent variable"""
    iv::T
    """Internal function representation of the basis field"""
    f_::Function
end

Base.show(io::IO, x::Basis) = print(io, "$(length(x.basis)) dimensional basis in ", "$(String.([v.op.name for v in x.variables]))")

@inline function Base.print(io::IO, x::Basis)
    show(io, x)
    println()
    if length(x.variables) == length(x.basis)
        for (i, bi) in enumerate(x.basis)
            println("d$(x.variables[i]) = $bi")
        end
    else
        for (i, bi) in enumerate(x.basis)
            println("f_$i = $bi")
        end
    end
end

is_independent(o::Operation) = isempty(o.args)

"""
    Basis(f, u; p, iv)

A basis over the variables `u` with parameters `p` and independent variable `iv`.
`f` can either be a Julia function which is able to use ModelingToolkit variables or
a vector of `Operation`.
It can be called with the typical DiffEq signature, meaning out of place with `f(u,p,t)`
or in place with `f(du, u, p, t)`.

# Example

```julia
using ModelingToolkit
using DataDrivenDiffEq

@parameters w[1:2] t
@variables u[1:2](t)

Ψ = Basis([u; sin.(w.*u)], u, parameters = p, iv = t)
```
"""
function Basis(basis::AbstractArray{Operation}, variables::AbstractArray{Operation};  parameters::AbstractArray =  Operation[], iv = nothing)
    @assert all(is_independent.(variables)) "Please provide independent variables for basis."

    bs = unique(basis)

    if isnothing(iv)
        @parameters t
        iv = t
    end

    vs = [ModelingToolkit.Variable(Symbol(i)) for i in variables]
    ps = [ModelingToolkit.Variable(Symbol(i)) for i in parameters]

    f_oop, f_iip = ModelingToolkit.build_function(bs, vs, ps, [iv], expression = Val{false})

    f_(u, p, t) = f_oop(u, p, t)
    f_(du, u, p, t) = f_iip(du, u, p, t)

    return Basis(bs, variables, parameters, iv, f_)
end

function Basis(basis::Function, variables::AbstractArray{Operation};  parameters::AbstractArray =  Operation[], iv = nothing)
    @assert all(is_independent.(variables)) "Please provide independent variables for basis."

    if isnothing(iv)
        @parameters t
        iv = t
    end

    try
        eqs = basis(variables, parameters, iv)
        return Basis(eqs, variables, parameters = parameters, iv = iv)
    catch e
        rethrow(e)
    end
end

function update!(basis::Basis)

    vs = [ModelingToolkit.Variable(Symbol(i))(basis.iv) for i in variables(basis)]
    ps = [ModelingToolkit.Variable(Symbol(i)) for i in parameters(basis)]

    f_oop, f_iip = ModelingToolkit.build_function(basis.basis, vs, ps, [basis.iv], expression = Val{false})

    f_(u, p, t) = f_oop(u, p, t)
    f_(du, u, p, t) = f_iip(du, u, p, t)

    basis.f_ = f_
    return
end

"""
    push!(basis, op)

    Push the operation(s) in `op` into the basis and update all internal fields accordingly.
    `op` can either be a single `Operation` or an Array of `Operation`s.
"""
function Base.push!(b::Basis, ops::AbstractArray{Operation})
    @inbounds for o in ops
        push!(b.basis, o)
    end
    unique!(b.basis)
    update!(b)
    return
end

function Base.push!(b::Basis, op₀::Operation)
    op = simplify(op₀)
    push!(b.basis, op)
    # Check for uniqueness
    unique!(b)
    update!(b)
    return
end

"""
    deleteat!(basis, inds)

    Delete the entries specified by `inds` and update the `Basis` accordingly.
"""
function Base.deleteat!(b::Basis, inds)
    deleteat!(b.basis, inds)
    update!(b)
    return
end

"""
    merge(basis_a, basis_b)

    Return a new `Basis`, which is defined via the union of both bases.
"""
function Base.merge(basis_a::Basis, basis_b::Basis)
    b =  unique(vcat(basis_a.basis, basis_b.basis))
    vs = unique(vcat(basis_a.variables, basis_b.variables))
    ps = unique(vcat(basis_a.parameter, basis_b.parameter))
    return Basis(b, vs, parameters = ps)
end

"""
    merge!(basis_a, basis_b)

    Updates the `Basis` to include the union of both bases.
"""
function Base.merge!(basis_a::Basis, basis_b::Basis)
    push!(basis_a, basis_b.basis)
    basis_a.variables = unique(vcat(basis_a.variables, basis_b.variables))
    basis_a.parameter = unique(vcat(basis_a.parameter, basis_b.parameter))
    update!(basis_a)
    return
end

Base.getindex(b::Basis, idx) = b.basis[idx]
Base.firstindex(b::Basis) = firstindex(b.basis)
Base.lastindex(b::Basis) = lastindex(b.basis)
Base.iterate(b::Basis) = iterate(b.basis)
Base.iterate(b::Basis, id) = iterate(b.basis, id)

function (==)(x::Basis, y::Basis)
    n = zeros(Bool, length(x.basis))
    @inbounds for (i, xi) in enumerate(x)
        n[i] = any(isequal.(xi, y.basis))
    end
    return all(n)
end

function count_operation(o::Expression, ops::AbstractArray)
    if isa(o, ModelingToolkit.Constant)
        return 0
    end
    k = o.op ∈ ops ? 1 : 0
    if !isempty(o.args)
        k += sum([count_operation(ai, ops) for ai in o.args])
    end
    return k
end

free_parameters(b::Basis; operations = [+]) = sum([count_operation(bi, operations) for bi in b.basis]) + length(b.basis)

(b::Basis)(u, p::DiffEqBase.NullParameters, t) = b(u, [], t)
(b::Basis)(du, u, p::DiffEqBase.NullParameters, t) = b(du, u, [], t)
(b::Basis)(u::AbstractVector,  p::AbstractArray = [], t = nothing) = b.f_(u, isempty(p) ? parameters(b) : p, isnothing(t) ? zero(eltype(u)) : t)
(b::Basis)(du::AbstractVector, u::AbstractVector, p::AbstractArray = [], t = nothing) = b.f_(du, u, isempty(p) ? parameters(b) : p, isnothing(t) ? zero(eltype(u)) : t)

function (b::Basis)(x::AbstractMatrix, p::AbstractArray = [], t::AbstractArray = [])
    isempty(t) ? nothing : @assert size(x, 2) == length(t)

    if (isempty(p) || eltype(p) <: Expression) && !isempty(parameters(b))
        pi = isempty(p) ? parameters(b) : p
        res = zeros(eltype(pi), length(b), size(x)[2])
    else
        pi = p
        res = zeros(eltype(x), length(b), size(x)[2])
    end

    @inbounds for i in 1:size(x)[2]
        res[:, i] .= b.f_(x[:, i], isempty(p) ? parameters(b) : p, isempty(t) ? zero(eltype(x)) : t[i])
    end

    return res
end

function (b::Basis)(y::AbstractMatrix, x::AbstractMatrix, p::AbstractArray = [], t::AbstractArray = [])
    @assert size(x, 2) == size(y, 2) "Provide consistent arrays."
    @assert size(y, 1) == length(b) "Provide consistent arrays."
    isempty(t) ? nothing : @assert size(x, 2) == length(t)

    @inbounds for i in 1:size(x, 2)
        b.f_(view(y, :, i), view(x, :, i), isempty(p) ? parameters(b) : p, isempty(t) ? zero(eltype(x)) : t[i])
    end

end

Base.size(b::Basis) = size(b.basis)
Base.length(b::Basis) = length(b.basis)

"""
    parameters(basis)

    Returns the parameters of the basis.
"""
ModelingToolkit.parameters(b::Basis) = b.parameter

"""
    variables(basis)

    Returns the dependent variables of the basis.
"""
variables(b::Basis) = b.variables

"""
    independent_variable(basis)

    Returns the independent_variable of the basis.
"""
ModelingToolkit.independent_variable(b::Basis) = b.iv

"""
    jacobian(basis)

    Returns a function representing the jacobian matrix / gradient of the `Basis` with respect to the
    dependent variables as a function with the common signature `f(u,p,t)` for out of place and `f(du, u, p, t)` for in place computation.
"""
function jacobian(basis::Basis)

    vs = [ModelingToolkit.Variable(Symbol(i))(independent_variable(basis)) for i in variables(basis)]
    ps = [ModelingToolkit.Variable(Symbol(i)) for i in parameters(basis)]

    j = calculate_jacobian(basis.basis, variables(basis))

    f_oop, f_iip = ModelingToolkit.build_function(expand_derivatives.(j), vs, ps, [basis.iv], expression = Val{false})

    f_(u, p, t) = f_oop(u, p, t)
    f_(du, u, p, t) = f_iip(du, u, p, t)

    return f_
end

function Base.unique!(b::Basis)
    N = length(b.basis)
    removes = Vector{Bool}()
    for i ∈ 1:N
        push!(removes, any([isequal(b.basis[i], b.basis[j]) for j in i+1:N]))
    end
    deleteat!(b, removes)
    update!(b)
end

function Base.unique(b::Basis)
    N = length(b.basis)
    returns = Vector{Bool}()
    for i ∈ 1:N
        push!(returns, any([isequal(b.basis[i], b.basis[j]) for j in i+1:N]))
    end
    returns = [!r for r in returns]
    return Basis(b.basis[returns], variables(b), parameters = parameters(b))
end

function Base.unique(b₀::AbstractVector{Operation})
    b = simplify.(b₀)
    N = length(b)
    returns = Vector{Bool}()
    for i ∈ 1:N
        push!(returns, any([isequal(b[i], b[j]) for j in i+1:N]))
    end
    returns = [!r for r in returns]
    return b[returns]
end

function Base.unique!(b::AbstractArray{Operation})
    N = length(b)
    removes = Vector{Bool}()
    for i ∈ 1:N
        push!(removes, any([isequal(b[i], b[j]) for j in i+1:N]))
    end
    deleteat!(b, removes)
end

"""
    dynamics(basis)

    Returns the internal function representing the dynamics of the `Basis`.
"""
function dynamics(b::Basis)
    return b.f_
end

"""
    ODESystem(basis)

    Converts the `Basis` into an `ODESystem` defined via `ModelingToolkit.jl`.
"""
function ModelingToolkit.ODESystem(b::Basis)
    @assert length(b) == length(variables(b))
    # Define the time
    @derivatives D'~independent_variable(b)

    vs = similar(b.variables)
    dvs = similar(b.variables)

    for (i, vi) in enumerate(b.variables)
        vs[i] = ModelingToolkit.Operation(vi.op, [independent_variable(b)])
        dvs[i] = D(vs[i])
    end
    eqs = dvs .~ b(vs, parameters(b), independent_variable(b))
    return ODESystem(eqs, independent_variable(b), variables(b), parameters(b))
end


text:
### ScaleTuner

# ScaleTuner tunes the scaling or stepsize parameters of a sampler's proposal step

immutable ScaleTuner <: AbstractTuner
    period::Int
    verbose::Bool
    targetrate::Real
    score::Function
    scoreargs::Tuple
    function ScaleTuner(period::Int,verbose::Bool,targetrate::Real,score::Function,scoreargs...)
        @assert 0 < period "Tuning period should be positive"
        @assert 0 <= targetrate <= 1 "Target acceptance rate should be between 0 and 1"
        new(period,verbose,targetrate,score,scoreargs)
    end
end

_tunerfunc(::Type{Val{:logistic}}) = (tunelogistic,7.0)
_tunerfunc(::Type{Val{:erf}}) = (tuneerf,3.0)

_tuner(::Type{Val{:scale}},period::Int,targetrate::Real,score::Symbol;verbose::Bool =true) = ScaleTuner(period,verbose,targetrate,_tunerfunc(Val{score})...)
_tuner(::Type{Val{:scale}},period::Int,targetrate::Real,score::Function,scoreargs...;verbose::Bool =true) = ScaleTuner(period,verbose,targetrate,score,scoreargs...)

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

type ScaleTunerState{T<:AbstractFloat} <: AbstractTunerState
    accepted::Vector{Int}
    proposed::Vector{Int}
    scalefactor::Vector{T}
    index::Int
    function ScaleTunerState(accepted::Vector{Int},proposed::Vector{Int},scalefactor::Vector{T},index::Int)
        @assert index > 0 "Total number of proposed MCMC samples should be non-negative"
        new(accepted,proposed,scalefactor,index)
    end
end

function _tunerstate{T<:AbstractFloat}(tuner::ScaleTuner,nburnin::Int,::Type{T})
    nsteps = ceil(Int,nburnin/tuner.period)
    ScaleTunerState{T}(zeros(Int,nsteps),zeros(Int,nsteps),zeros(T,nsteps),1)
end

function tune(tuner::ScaleTuner,state::ScaleTunerState)
    r = state.accepted[state.index]/state.proposed[state.index]
    s = tuner.score(r-tuner.targetrate,tuner.scoreargs...)
    state.scalefactor[state.index] = convert(eltype(state.scalefactor),s)
end

function show(io::IO,tuner::ScaleTuner)
    println(io,"ScaleTuner: period = $(tuner.period), verbose = $(tuner.verbose) ,targetrate = $(tuner.targetrate)")
    println(io,"            score = $(tuner.score)",isempty(tuner.scoreargs)?"":", scoreargs = $(tuner.scoreargs)")
end

function show(io::IO,state::ScaleTunerState)
    println(io,"ScaleTunerState: ")
    println(io," accepted = $(accepted(state))")
    println(io," proposed = $(proposed(state))")
    println(io," scalefactor = $(round(state.scalefactor,3))")
    println(io," acceptance rate = $(round(rate(state),3))")
    println(io," total proposed = $(total(state))")
end

function showstep(t::ScaleTuner,state::ScaleTunerState)
    if verbose(t) && state.index <= numtunesteps(state)
        a,p = current(state)
        println("  accepted/proposed = $a/$p")
        println("  acceptance rate = $(round(a/p,3))")
        println("  scalefactor = $(round(state.scalefactor[state.index],3))")
        println("  cummulative scaling = $(round((r = one(eltype(state.scalefactor)) ; for i=1:state.index r*=state.scalefactor[i] end ; r),3))")
    end
end


text:
# save current status to file
function checkpoint!(
    file :: HDF5.File,
    Λ    :: Float64,
    dΛ   :: Float64,
    m    :: Mesh,
    a    :: Action_su2
    )    :: Nothing

    # save step size
    file["dΛ/$(Λ)"] = dΛ

    # save frequency meshes
    file["σ/$(Λ)"]  = m.σ
    file["Ωs/$(Λ)"] = m.Ωs
    file["νs/$(Λ)"] = m.νs
    file["Ωt/$(Λ)"] = m.Ωt
    file["νt/$(Λ)"] = m.νt
    file["χ/$(Λ)"]  = m.χ

    # save spin length
    if haskey(file, "S") == false
        file["S"] = a.S
    end

    # save symmetry group
    if haskey(file, "symmetry") == false
        file["symmetry"] = "su2"
    end

    # save self energy
    file["a/$(Λ)/Σ"] = a.Σ

    # save vertex
    save!(file, "a/$(Λ)/Γ/spin", a.Γ[1])
    save!(file, "a/$(Λ)/Γ/dens", a.Γ[2])

    return nothing
end

# read checkpoint from file
function read_checkpoint_su2(
    file :: HDF5.File,
    Λ    :: Float64
    )    :: Tuple{Float64, Float64, Mesh, Action_su2}

    # filter out nearest available cutoff
    list    = keys(file["σ"])
    cutoffs = parse.(Float64, list)
    index   = argmin(abs.(cutoffs .- Λ))
    println("Λ was adjusted to $(cutoffs[index]).")

    # read step size
    dΛ = read(file, "dΛ/$(cutoffs[index])")

    # read frequency meshes
    σ  = read(file, "σ/$(cutoffs[index])")
    Ωs = read(file, "Ωs/$(cutoffs[index])")
    νs = read(file, "νs/$(cutoffs[index])")
    Ωt = read(file, "Ωt/$(cutoffs[index])")
    νt = read(file, "νt/$(cutoffs[index])")
    χ  = read(file, "χ/$(cutoffs[index])")
    m  = Mesh(length(σ), length(Ωs), length(νs), length(χ), σ, Ωs, νs, Ωt, νt, χ)

    # read spin length
    S = read(file, "S")

    # read self energy
    Σ = read(file, "a/$(cutoffs[index])/Σ")

    # read vertex
    Γ = Vertex[read_vertex(file, "a/$(cutoffs[index])/Γ/spin"),
               read_vertex(file, "a/$(cutoffs[index])/Γ/dens")]

    # build action
    a = Action_su2(S, Σ, Γ)

    return cutoffs[index], dΛ, m, a
end


text:
@overload
function Concepts.forward_map(distribution::Union{AbstractPoisson,Type{Val{:Poisson}}},
                              canonical_parameter::Array{T};
                              non_canonical_parameter::Union{Array{Float64},Nothing} = nothing,
                              non_canonical_map = nothing) where T<:Real
  if !isnothing(non_canonical_parameter)
    ## TODO
  end
  return exp.(canonical_parameter)
end

@overload
function Concepts.forward_map(distribution::Union{AbstractGamma,Type{Val{:Gamma}}},
                              canonical_parameter::Array{T};
                              non_canonical_parameter::Union{Array{Float64},Nothing} = nothing,
                              non_canonical_map = nothing) where T<:Real
  if !isnothing(non_canonical_parameter)
    ## TODO
  end

  canonical_parameter = 1 ./ exp.(canonical_parameter)
  return canonical_parameter
  # canonical_parameter = -1 .* abs.(canonical_parameter)
  # return -1 ./ canonical_parameter
  # return 1 ./ canonical_parameter
end

@overload
function Concepts.forward_map(distribution::Union{AbstractNegativeBinomial,Type{Val{:NegativeBinomial}}},
                              canonical_parameter::Array{T};
                              non_canonical_parameter::Union{Array{Float64},Nothing} = nothing,
                              non_canonical_map = nothing,
                              r_estimate = nothing) where T<:Real
  if !isnothing(non_canonical_parameter)
    ## TODO
  end
  # @show(r_estimate)
  return r_estimate ./ (exp.(exp.(canonical_parameter)) .- 1)
end

@overload
function Concepts.forward_map(distribution::Union{AbstractBernoulli,Type{Val{:Bernoulli}}},
                              canonical_parameter::Array{T};
                              non_canonical_parameter::Union{Array{Float64},Nothing} = nothing,
                              non_canonical_map = nothing) where T<:Real
  if !isnothing(non_canonical_parameter)
    ## TODO
  end
  ex = exp.(canonical_parameter)
  return ex ./ (1 .+ ex)
  #    return (Int.(sign.(canonical_parameter)) .+ 1) ./ 2
end

@overload
function Concepts.forward_map(distribution::Union{AbstractGaussian,Type{Val{:Gaussian}}},
                              canonical_parameter::Array{T};
                              non_canonical_parameter::Union{Array{Float64},Nothing} = nothing,
                              non_canonical_map = nothing) where T<:Real
  if !isnothing(non_canonical_parameter)
    ## TODO
  end
  return canonical_parameter
end

@overload
function Concepts.forward_map(distribution::Symbol,
                              canonical_parameter::Array{T};
                              non_canonical_parameter::Union{Array{Float64},Nothing} = nothing,
                              non_canonical_map = nothing) where T<:Real
  return Concepts.forward_map(Val{distribution},canonical_parameter,
                              non_canonical_map=non_canonical_map,
                              non_canonical_parameter=non_canonical_parameter)
end

@overload
function Concepts.predict(distribution::Union{AbstractPoisson,Type{Val{:Poisson}}},mean::Any;
                          custom_prediction_function=nothing)
  if !isnothing(custom_prediction_function)
    return -1.0
  end
  return round.(mean)
end

@overload
function Concepts.predict(distribution::Union{AbstractBernoulli,Type{Val{:Bernoulli}}},mean::Any;
                          custom_prediction_function=nothing)
  if !isnothing(custom_prediction_function)
    return -1.0
  end
  return Int.(mean .> 0.5)
end

@overload
function Concepts.predict(distribution::Union{AbstractGaussian,Type{Val{:Gaussian}}},mean::Any;
                          custom_prediction_function=nothing)
  if !isnothing(custom_prediction_function)
    return -1.0
  end
  return mean
end

@overload
function Concepts.predict(distribution::Union{AbstractGamma,Type{Val{:Gamma}}},mean::Any;
                          custom_prediction_function=nothing)
  if !isnothing(custom_prediction_function)
    return -1.0
  end
  return mean
end

@overload
function Concepts.predict(distribution::Union{AbstractNegativeBinomial,Type{Val{:NegativeBinomial}}},mean::Any;
                          custom_prediction_function=nothing)
  if !isnothing(custom_prediction_function)
    return -1.0
  end
  return round.(mean)
end

@overload
const Concepts.predict(obj::Symbol,arg1;custom_prediction_function=nothing) =
  Concepts.predict(Val{obj},arg1;custom_prediction_function=custom_prediction_function)


text:
using Parquet
using Test

function test_col_cursor(file::String)
    p = ParFile(file)
    println("loaded ", file)

    nr = nrows(p)
    cnames = colnames(p)
    for cname in cnames
        rr = 1:nr
        println("\tcolumn ", cname, " rows:", rr)
        println("\tvalue, defn, repn, next idx")
        t1 = time()
        cc = Parquet.ColCursor(p, rr, cname)
        num_read = 0
        for (v,i) in enumerate(cc)
            val,defn,repn = v
            num_read += 1
        end
        println("\t\t", isnull(val) ? nothing : get(val), ", ", defn, ", ", repn, ", ", i)
        println("\t\tread $num_read values in $(time()-t1) time")
    end
end

function test_row_cursor(file::String, map_logical_types::Bool)
    p = ParFile(file; map_logical_types=map_logical_types)

    t1 = time()
    nr = nrows(p)
    cnames = colnames(p)
    rc = RecordCursor(p)
    rec = nothing
    for i in rc
        rec = i
    end
    @info("loaded", file, count=nr, last_record=rec, time_to_read=time()-t1)
end

function test_batchedcols_cursor(file::String, map_logical_types::Bool)
    p = ParFile(file; map_logical_types=map_logical_types)

    t1 = time()
    nr = nrows(p)
    cnames = colnames(p)
    cc = BatchedColumnsCursor(p)
    batch = nothing
    for i in cc
        batch = i
    end
    @info("loaded", file, count=nr, ncols=length(propertynames(batch)), time_to_read=time()-t1)
end

function test_col_cursor_all_files()
    for encformat in ("SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            test_col_cursor(joinpath(@__DIR__, "parquet-compatibility", "parquet-testdata", "impala", "1.1.1-$encformat/$fname.impala.parquet"))
        end
    end
end

function test_row_cursor_all_files()
    for encformat in ("SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            for map_logical_types in (true,false)
                test_row_cursor(joinpath(@__DIR__, "parquet-compatibility", "parquet-testdata", "impala", "1.1.1-$encformat/$fname.impala.parquet"), map_logical_types)
            end
        end
    end
end

function test_batchedcols_cursor_all_files()
    for encformat in ("SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            for map_logical_types in (true,false)
                test_batchedcols_cursor(joinpath(@__DIR__, "parquet-compatibility", "parquet-testdata", "impala", "1.1.1-$encformat/$fname.impala.parquet"), map_logical_types)
            end
        end
    end
end

#test_col_cursor_all_files()
test_row_cursor_all_files()
test_batchedcols_cursor_all_files()


text:
# Compares SQUAREEM variants on given network structure
# Example of commaind line call:
#     julia --threads=4 --color=yes squarem.jl assets/nltcs.spn assets/nltcs.train.csv assets/nltcs.valid.csv
#using DataFrames, CSV
using DelimitedFiles # much faster than above
using Random
using SPNetworks
import SPNetworks: NLL
import SPNetworks.ParameterLearning: SQUAREM, initialize, converged, update

if length(ARGS) < 3
    println("Usage: julia --color=yes em.jl spn_filename train_filename validation_filename")
    exit()
end
spn_filename = normpath(ARGS[1])
train_filename = normpath(ARGS[2])
valid_filename = normpath(ARGS[3])
println("Threads: ", Threads.nthreads())

println("SPN: $spn_filename")
spn = SumProductNetwork(spn_filename)
@show summary(spn)
nvars = length(scope(spn))
println("Training data: $train_filename")
# Load datasets
@time tdata = readdlm(train_filename, ',', Float64) .+ 1;
@show summary(tdata)
println("Validation data: $valid_filename")
@time vdata = readdlm(valid_filename, ',', Float64) .+ 1;
@show summary(vdata)
# initialize EM learner
learner = SQUAREM(spn)
Random.seed!(3)
initialize(learner) # generate random weights for each sum node
println("It: $(learner.steps) \t NLL: $(NLL(spn, tdata)) \t held-out NLL: $(NLL(spn, vdata))")
# Running Expectation Maximization
while !converged(learner) && learner.steps < 21
    @time score, α = update(learner, tdata)
    testnll = NLL(spn, vdata)
    trainnll = NLL(spn, tdata)
    println("It: $(learner.steps) \t NLL: $trainnll \t held-out NLL: $testnll \t α: $α")
end


text:
module IntegrationWrapper

using JavaCall

const BaseAbstractUnivariateIntegrator = @jimport org.hipparchus.analysis.integration.BaseAbstractUnivariateIntegrator
const Class = @jimport java.lang.Class
const IterativeLegendreGaussIntegrator = @jimport org.hipparchus.analysis.integration.IterativeLegendreGaussIntegrator
const JString = @jimport java.lang.String
const MidPointIntegrator = @jimport org.hipparchus.analysis.integration.MidPointIntegrator
const Object = @jimport java.lang.Object
const RombergIntegrator = @jimport org.hipparchus.analysis.integration.RombergIntegrator
const SimpsonIntegrator = @jimport org.hipparchus.analysis.integration.SimpsonIntegrator
const TrapezoidIntegrator = @jimport org.hipparchus.analysis.integration.TrapezoidIntegrator
const UnivariateFunction = @jimport org.hipparchus.analysis.UnivariateFunction
const UnivariateIntegrator = @jimport org.hipparchus.analysis.integration.UnivariateIntegrator

export BaseAbstractUnivariateIntegrator
export IterativeLegendreGaussIntegrator
export MidPointIntegrator
export RombergIntegrator
export SimpsonIntegrator
export TrapezoidIntegrator
export UnivariateIntegrator
export equals
export get_absolute_accuracy
export get_class
export get_evaluations
export get_iterations
export get_maximal_iteration_count
export get_minimal_iteration_count
export get_relative_accuracy
export hash_code
export integrate
export notify
export notify_all
export to_string
export wait

include("GaussWrapper/GaussWrapper.jl")
include("base_abstract_univariate_integrator.jl")
include("iterative_legendre_gauss_integrator.jl")
include("mid_point_integrator.jl")
include("romberg_integrator.jl")
include("simpson_integrator.jl")
include("trapezoid_integrator.jl")
include("univariate_integrator.jl")

end


text:
@testset "pipe" begin

    # basic usage
    @test isa(pipe, Function)
    p = pipe(
        (v) -> next(v + 1),
        (v) -> done(v * 2) )

    @test isa(p, Pipeline)

    c = run(p, 4)

    @test isa(c, Continuation) == true
    @test isdone(c) == true
    @test c.data == (10,)

    # short-circuit
    p = pipe(
        (v) -> done(v + 1),
        (v) -> done(v * 2) )

    @test isa(p, Pipeline)

    c = run(p, 4)

    @test isa(c, Continuation) == true
    @test isdone(c) == true
    @test c.data == (5,)

    # concatenation
    p2 = pipe(v -> next(v + 1))

    p = pipe(
        p2,
        (v) -> done(v * 2) )

    @test isa(p, Pipeline)

    c = run(p, 4)

    @test isa(c, Continuation) == true
    @test isdone(c) == true
    @test c.data == (10,)

    # multiple parameters
    p = pipe(
        (v1, v2) -> next(v1 + 2, v1 + v2),
        (v1, v2) -> next(v1 * 2, v2 + 7) )

    c = run(p, 1, 1)
    @test c.data == (6, 9)
end

@testset "resume" begin
    # normal usage
    p1 = pipe(v -> done(v + 1))
    p2 = pipe(v -> done(v * 2))

    p3 = resume(p1, p2)

    @test isa(p3, Pipeline) == true

    c = run(p3, 3)

    @test isa(c, Continuation) == true
    @test c.data == (8,)

    # when not done
    p1 = pipe(v -> next(v + 1))
    p2 = pipe(v -> done(v * 2))

    p3 = resume(p1, p2)

    @test isa(p3, Pipeline) == true

    c = run(p3, 3)

    @test isa(c, Continuation) == true
    @test c.data == (4,)

    # when not done, but always set to true
    p1 = pipe(v -> next(v + 1))
    p2 = pipe(v -> done(v * 2))

    p3 = resume(p1, p2; always=true)

    @test isa(p3, Pipeline) == true

    c = run(p3, 3)

    @test isa(c, Continuation) == true
    @test c.data == (8,)

    # when done, and always set to true
    p1 = pipe(v -> done(v + 1))
    p2 = pipe(v -> done(v * 2))

    p3 = resume(p1, p2; always=true)

    @test isa(p3, Pipeline) == true

    c = run(p3, 3)

    @test isa(c, Continuation) == true
    @test c.data == (8,)

    # with a function
    p1 = pipe(v -> done(v + 1))
    p2 = v -> done(v * 2)

    p3 = resume(p1, p2)

    @test isa(p3, Pipeline) == true

    c = run(p3, 3)

    @test isa(c, Continuation) == true
    @test c.data == (8,)

    # with multiple arguments
    p1 = pipe(v -> done(v + 1))
    p2 = pipe(v -> next(v * 2))

    p3 = resume(p1, p2, v -> done(v + 9))

    @test isa(p3, Pipeline) == true

    c = run(p3, 3)

    @test isa(c, Continuation) == true
    @test c.data == (17,)

end

@testset "recover" begin
    # with no exception
    p = pipe(v -> done(v + 2))
    p = recover(p, (err, v) -> done(v * 12))
    @test isa(p, Pipeline) == true

    c = run(p, 4)
    @test isa(c, Continuation) == true

    @test isdone(c) == true
    @test c.data == (6,)

    # with exception
    p = pipe(v -> throw(UndefVarError(:sushi)))
    p = recover(p, (err, v) -> begin
        @test isa(err, UndefVarError)
        @test v == 4
        done(v * 12)
    end)
    @test isa(p, Pipeline) == true

    c = run(p, 4)
    @test isa(c, Continuation) == true

    @test isdone(c) == true
    @test c.data == (48,)
end


text:
function common_prefix(completions)
    ret = ""
    i = nexti = 1
    cc, nexti = next(completions[1], 1)
    while true
        for c in completions
            (i > length(c) || c[i] != cc) && return ret
        end
        ret *= string(cc)
        i >= length(completions[1]) && return ret
        i = nexti
        cc, nexti = next(completions[1], i)
    end
end

# Show available completions
function show_completions(s::PromptState, completions)
    colmax = maximum(map(length, completions))
    num_cols = max(div(width(terminal(s)), colmax+2), 1)
    entries_per_col, r = divrem(length(completions), num_cols)
    entries_per_col += r != 0
    # skip any lines of input after the cursor
    cmove_down(terminal(s), input_string_newlines_aftercursor(s))
    println(terminal(s))
    for row = 1:entries_per_col
        for col = 0:num_cols
            idx = row + col*entries_per_col
            if idx <= length(completions)
                cmove_col(terminal(s), (colmax+2)*col)
                print(terminal(s), completions[idx])
            end
        end
        println(terminal(s))
    end
    # make space for the prompt
    for i = 1:input_string_newlines(s)
        println(terminal(s))
    end
end

max_stars_repo_name: dadarakt/jamon

max_issues_repo_name: dadarakt/jamon

max_forks_repo_name: dadarakt/jamon


text:
@assert VERSION >= v"1.0" "Your Julia version does not satisfy the requirements for this package. You need to use Julia v1.x"
t0 = time_ns()
# cd(string(ENV["HOMEPATH"],"\\Documents\\ASync\\home\\Code\\Julia\\DecisionTrees.jl"))

using Distributed
# addprocs(22)
# Distributed.@everywhere using Revise
Distributed.@everywhere import CSV
Distributed.@everywhere import DataFrames
Distributed.@everywhere import DataFrames: DataFrame

Distributed.@everywhere using DecisionTrees  # this may take some time

tela = (time_ns() - t0) / 1e9
@show tela # precompilation (first time usage) may take some time

############################################################
# 1. Model for the French MTPL Data
############################################################

##############################
# Read the data
##############################
datafile = joinpath("data", "freMTPL2", "freMTPL2.csv")
@assert isfile(datafile) "You need to unzip the file freMTPL2.zip in the folder data/freMTPL2"
@time fullData = CSV.read(datafile, rows_for_type_detect=100000, allowmissing=:none, categorical=false);

areasSorted = sort(unique(fullData[:Area]))
AreaInteger = map(x->findall((in)([x]), areasSorted)[1], fullData[:Area])
fullData[:AreaInteger] = AreaInteger

# correct for unreasonable observations
for i = 1:size(fullData, 1)
    fullData[:ClaimNb][i] = min(4, fullData[:ClaimNb][i])
    fullData[:Exposure][i] = min(1.0, fullData[:Exposure][i])
end

# set independent variables
selected_explanatory_vars = ["Area","AreaInteger","VehPower","VehAge","DrivAge","BonusMalus","VehBrand","VehGas","Density","Region"]
# note: we refrain from using the exposure as explanatory variable, as this is not meaningful in practical applications
# where one aims to predict the claim frequency of a given policy (without knowing how long that policy will remain within a certain portfolio)

dtmtable, sett, dfprepped = prepare_dataframe_for_dtm!(fullData, keycol="IDpol", trnvalcol="trnTest", numcol="ClaimNb", denomcol="Exposure", weightcol="Exposure", independent_vars=selected_explanatory_vars);

updateSettingsMod!(sett,minWeight=-0.03,model_type="build_tree",calculatePoissonError=true)

resultingFiles, resM = dtm(dtmtable, sett)


text:
function read_model_configuration(pm::PredictionModel{BMTMKLModel}, data::Dict{String, Any})
	if haskey(data, "model_configuration")
		mc_data = data["model_configuration"]
		# general setting
		if haskey(mc_data, "alpha") && haskey(mc_data, "beta") && haskey(mc_data, "mu") && haskey(mc_data, "sigma2")
			return ModelConfiguration(pm, mc_data["alpha"], mc_data["beta"], mc_data["mu"], mc_data["sigma2"])
		else
			mc = ModelConfiguration(pm)
			# refined settings if present
			if haskey(mc_data, "bias_precision_alpha") && haskey(mc_data, "bias_precision_beta")
				mc.parameters["α_ɣ"] = mc_data["bias_precision_alpha"]
				mc.parameters["β_ɣ"] = mc_data["bias_precision_beta"]
			end
			if haskey(mc_data, "weights_precision_alpha") && haskey(mc_data, "weights_precision_beta")
				mc.parameters["α_λ"] = mc_data["weights_precision_alpha"]
				mc.parameters["β_λ"] = mc_data["weights_precision_beta"]
			end
			if haskey(mc_data, "outcome_precision_alpha") && haskey(mc_data, "outcome_precision_beta")
				mc.parameters["α_ε"] = mc_data["outcome_precision_alpha"]
				mc.parameters["β_ε"] = mc_data["outcome_precision_beta"]
			end
			if haskey(mc_data, "intermed_results_precision_alpha") && haskey(mc_data, "intermed_results_precision_beta")
				mc.parameters["α_ν"] = mc_data["intermed_results_precision_alpha"]
				mc.parameters["β_ν"] = mc_data["intermed_results_precision_beta"]
			end
			if haskey(mc_data, "kernel_weights_precision_alpha") && haskey(mc_data, "kernel_weights_precision_beta")
				mc.parameters["α_⍵"] = mc_data["kernel_weights_precision_alpha"]
				mc.parameters["β_⍵"] = mc_data["kernel_weights_precision_beta"]
			end
			if haskey(mc_data, "bias_mean") && haskey(mc_data, "bias_sigma2")
				mc.parameters["μ_b"] = mc_data["bias_mean"]
				mc.parameters["σ_b"] = mc_data["bias_sigma2"]
			end
			if haskey(mc_data, "kernel_weights_mean") && haskey(mc_data, "kernel_weights_sigma2")
				mc.parameters["μ_e"] = mc_data["kernel_weights_mean"]
				mc.parameters["σ_e"] = mc_data["kernel_weights_sigma2"]
			end
			if haskey(mc_data, "weights_mean") && haskey(mc_data, "weights_sigma2")
				mc.parameters["μ_a"] = mc_data["weights_mean"]
				mc.parameters["Σ_a"] = mc_data["weights_sigma2"]
			end
			if haskey(mc_data, "intermed_results_mean") && haskey(mc_data, "intermed_results_sigma2")
				mc.parameters["μ_g"] = mc_data["intermed_results_mean"]
				mc.parameters["Σ_g"] = mc_data["intermed_results_sigma2"]
			end
		end
	end
	# standard settings if no config present
	ModelConfiguration(pm)
end

function init!(m::PredictionModel{BMTMKLModel}, ic::InferenceConfiguration)
	(base_kernels, pathway_specific_kernels) = compute_all_kernels(m.data, collect(values(m.data.cell_lines)), subsume_pathways = ic.subsume_pathways)
	m.precomputations["base_kernels"] = base_kernels
	m.precomputations["pathway_specific_kernels"] = pathway_specific_kernels
end

function update_model_config!(pm::PredictionModel{BMTMKLModel}, mc::ModelConfiguration)
	mc.parameters["T"] = length(pm.data.drugs)
	mc.parameters["N"] = OrderedDict{Drug, Int64}()
	for drug in collect(keys(pm.data.results))
		# N[t] is the number of cell lines for which we have a result and that are in the training set
		mc.parameters["N"][drug] = length(filter((cl) -> !cl.in_test_set && haskey(pm.data.results[drug].outcome_values, cl), collect(values(pm.data.cell_lines))))
	end
	# recompute overall number of kernels (+4 accomodates for combined kernels)
	mc.parameters["K"] = length(pm.precomputations["base_kernels"]) + sum(map(length, collect(values(pm.precomputations["pathway_specific_kernels"])))) + 4
end

function post_init!(m::PredictionModel{BMTMKLModel})
    # collect outcomes statistics, only on results used for training
    for (t, drug) in enumerate(keys(m.data.results))
        vals = Float64[]
        [if !cl.in_test_set push!(vals, m.data.results[drug].outcome_values[cl]) end for cl in keys(m.data.results[drug].outcome_values)]
        m.data.results[drug].outcome_mean = mean(vals)
        m.data.results[drug].outcome_std = stdm(vals, m.data.results[drug].outcome_mean)
    end
    all_drugs = collect(keys(m.data.results)) # the tasks
    cell_lines = collect(values(m.data.cell_lines)) # all cell lines in the experiment, excluding held-out data

    # the kernels between cell lines
    kernels = DataStructures.OrderedDict{Drug, Vector{Matrix{Float64}}}()
    # the target values
    targets = DataStructures.OrderedDict{Drug, Vector{Float64}}()

    # cross kernels between training data and test data
    cross_kernels = DataStructures.OrderedDict{Drug, Vector{Matrix{Float64}}}()
    # target values for test data
    test_targets = DataStructures.OrderedDict{Drug, Vector{Float64}}()

    for (t, drug) in enumerate(all_drugs)
        #init
        kernels[drug] = Vector{Matrix{Float64}}()
        cross_kernels[drug] = Vector{Matrix{Float64}}()

        # determine cell lines available for this drug
        result_cell_lines = filter((cl) ->  !cl.in_test_set && haskey(m.data.results[drug].outcome_values, cl), cell_lines)
        test_result_cell_lines = filter((cl) -> cl.in_test_set && haskey(m.data.results[drug].outcome_values, cl), cell_lines)

        # get their indicex in original cell line array
        idx_in_cell_lines = findall((in)(result_cell_lines), cell_lines)
        test_idx_in_cell_lines = findall((in)(test_result_cell_lines), cell_lines)

        #construct the kernels
        for v in m.data.views
            push!(kernels[drug], m.precomputations["base_kernels"][v][idx_in_cell_lines, idx_in_cell_lines])
            push!(cross_kernels[drug], m.precomputations["base_kernels"][v][test_idx_in_cell_lines, idx_in_cell_lines])
            if length(m.data.pathway_information) != 0
                for pw_kernel in m.precomputations["pathway_specific_kernels"][v]
                    push!(kernels[drug], pw_kernel[idx_in_cell_lines, idx_in_cell_lines])
                    push!(cross_kernels[drug], pw_kernel[test_idx_in_cell_lines, idx_in_cell_lines])
                end
            end
        end
        # compute additional kernels
        gene_expression = m.precomputations["base_kernels"][GeneExpression][idx_in_cell_lines, idx_in_cell_lines]
        methylation = m.precomputations["base_kernels"][Methylation][idx_in_cell_lines, idx_in_cell_lines]
        cnv = m.precomputations["base_kernels"][CNV][idx_in_cell_lines, idx_in_cell_lines]

        push!(kernels[drug], gene_expression .* methylation)
        push!(kernels[drug], gene_expression .* cnv)
        push!(kernels[drug], cnv .* methylation)
        push!(kernels[drug], gene_expression .* methylation .* cnv)

        cross_gene_expression = m.precomputations["base_kernels"][GeneExpression][test_idx_in_cell_lines, idx_in_cell_lines]
        cross_methylation = m.precomputations["base_kernels"][Methylation][test_idx_in_cell_lines, idx_in_cell_lines]
        cross_cnv = m.precomputations["base_kernels"][CNV][test_idx_in_cell_lines, idx_in_cell_lines]

        push!(cross_kernels[drug], cross_gene_expression .* cross_methylation)
        push!(cross_kernels[drug], cross_gene_expression .* cross_cnv)
        push!(cross_kernels[drug], cross_cnv .* cross_methylation)
        push!(cross_kernels[drug], cross_gene_expression .* cross_methylation .* cross_cnv)

        train_outcome = map(cl -> m.data.results[drug].outcome_values[cl], result_cell_lines)
        test_outcome = map(cl -> m.data.results[drug].outcome_values[cl], test_result_cell_lines)

        #normalize the outcome data for each drug across the cell lines
        targets[drug] = (train_outcome .- m.data.results[drug].outcome_mean) ./ m.data.results[drug].outcome_std

        # test outcomes, not normalized
        test_targets[drug] = test_outcome
    end

	m.precomputations["kernels"] = kernels
	m.precomputations["cross_kernels"] = cross_kernels
	m.precomputations["targets"] = targets
	m.precomputations["test_targets"] = test_targets
	log_message("done post init")
	nothing
 end

function compute_all_kernels(experiment::Experiment, cell_lines::Vector{CellLine}, cell_lines_test::Union{Vector{CellLine}, Nothing}=nothing; subsume_pathways::Bool=true)
	log_message("in compute kernel method")
	base_kernels = OrderedDict{Type{<:ViewType}, Matrix{Float64}}()
    pathway_specific_kernels = OrderedDict{Type{<:ViewType}, Vector{Matrix{Float64}}}()
    # compute a base kernel for each view, containing all cell lines that are present in all views
    num_views = length(experiment.views)
    for v in experiment.views
        # get views of this type for all cell lines
        # v = experiment.views[v_id]
        data_views = [map_data_views(cell_lines, v)]
        if cell_lines_test != nothing
            push!(data_views, map_data_views(cell_lines_test, v))
        end
        # compute the base kernel for this view
        kct = @elapsed base_kernels[v] = compute_kernel(data_views...)
        if cell_lines_test == nothing
            log_message("computing base kernel for $v $kct seconds)")
        end
        if length(experiment.pathway_information) != 0
            if cell_lines_test == nothing
                log_message("computing pathway specific kernels...")
            end
            pathway_specific_kernels[v] = Vector{Matrix{Float64}}()

            # # pathway specific kernels method 1: compute a kernel for each pathway
            if !subsume_pathways
                num_pws = length(experiment.pathway_information)
                tt = 0
                for pathway in collect(values(experiment.pathway_information))
                    tt += @elapsed push!(pathway_specific_kernels[v], compute_kernel(data_views..., pathway=pathway))
                end

                # method 1 end
            else
                # pathway specific kernels method 2:
                # - for each cell line:
                #   - for each pathway:
                #       - compute average (expression data) or max value (otherwise) over pathway genes
                # - compute cell line similarities as usual
                tt = @elapsed begin
                    cell_line_gene_set_views = map((x) -> Vector{Vector{Float64}}(undef, length(x)), data_views)
                    for (cl_set_idx, cl_set) in enumerate(data_views)
                        for (cl_idx, dv) in enumerate(cl_set)
                            cell_line_gene_set_views[cl_set_idx][cl_idx] = zeros(length(experiment.pathway_information))
                            for (pw_idx, pathway) in enumerate(values(experiment.pathway_information))

                                (m_values, _) = prepare_kernel(dv, dv, pathway)
                                if v in [RNASeq, GeneExpression]
                                    cell_line_gene_set_views[cl_set_idx][cl_idx][pw_idx] = mean(m_values)
                                else
                                    cell_line_gene_set_views[cl_set_idx][cl_idx][pw_idx] = maximum(m_values)
                                end
                            end
                        end
                    end
                    len_views = length(cell_line_gene_set_views)
                    pw_kernel = zeros(length(cell_line_gene_set_views[1]), length(cell_line_gene_set_views[len_views]))
                    for i in 1:size(pw_kernel, 1), j in 1:size(pw_kernel, 2)
                        pw_kernel[i,j] = kernel_function(cell_line_gene_set_views[1][i], cell_line_gene_set_views[len_views][j])
                    end
                    push!(pathway_specific_kernels[v], pw_kernel)
                end
            # # method 2 end
            end
            if cell_lines_test == nothing
                log_message("computing pathway specific kernels for $v took $tt seconds)")
            end
        else
            if cell_lines_test == nothing
                log_message("Note: no pathway information provided")
            end
        end
    end
    if cell_lines_test == nothing
        log_message("done computing kernels")
    end
    (base_kernels, pathway_specific_kernels)
end

function map_data_views(cell_lines::Vector{CellLine}, view_type::Type{<:ViewType})
    ret = Vector{DataView{<:KeyType, <:ViewType}}()
    for cl in cell_lines
        if haskey(cl.views, view_type)
            push!(ret, cl.views[view_type])
        else
            if view_type == RPPA
                push!(ret, DataView{Protein, NAViewType}(cl.id))
            else
                push!(ret, DataView{Gene, NAViewType}(cl.id))
            end
        end
    end
    ret
end

get_complete_cell_lines(experiment::Experiment) = filter(cl -> length(cl.views) == length(experiment.views), collect(values(experiment.cell_lines)))

function compute_kernel(dvs::Vector{DataView}; pathway::Union{Pathway, Nothing} = nothing)
    N = length(dvs)
    k = zeros(N,N)
    prep_time = 0
    comp_time = 0
    kernel_compute_time = @elapsed begin
        Threads.@threads for i in 1:N
            for  j in i:N
                kpt, kct, k[i,j] = compute_kernel_value(dvs[i], dvs[j], pathway)
                prep_time += kpt
                comp_time += kct
                if isnan(k[i,j])
                    k[i,j] = 1e-7
                end
                if i≠j k[j,i] = k[i,j] end
            end
        end
        if !isposdef(k)
            h_k = Array(Hermitian(k))
            # @warn "kernel is not positive definite, using only upper triangular form" isposdef(h_k) sum_diff=sum(abs.(h_k - k))
            k = h_k
        end
    end
    log_message("computed kernel ($(N*(N-1)) computations), preparation time=$prep_time computation time=$comp_time")
    k
end

function compute_kernel(dvs_train::Vector{DataView}, dvs_test::Vector{DataView}; pathway::Union{Pathway, Nothing}=nothing)
    N = length(dvs_train)
    M = length(dvs_test)
    k = zeros(N,M)
    Threads.@threads for i in 1:N
        for j in 1:M
            kpt, kct, k[i,j] = compute_kernel_value(dvs_train[i], dvs_test[j], pathway)
            if isnan(k[i,j])
                k[i,j] = 1e-7
            end
        end
    end
    k
end

function compute_kernel_value(a::DataView, b::DataView, pathway::Union{Pathway, Nothing})
    kpc = @elapsed kernel_params = prepare_kernel(a,b,pathway)
    kcc = @elapsed kfv = kernel_function(kernel_params...)
    kpc, kcc, kfv
end

# general view type vs. NAViewType
function prepare_kernel(a::DataView{K, <:ValueViewType}, b::DataView{K, NAViewType}, pathway::Union{Pathway, Nothing}) where {K<:KeyType}
    keys_to_use = get_common_keys(a, b, pathway)
    (map(m -> a.measurements[m].normalized_value, keys_to_use), zeros(length(keys_to_use)))
end

# NAViewType vs. general view type
function prepare_kernel(a::DataView{K, NAViewType}, b::DataView{K, <:ValueViewType}, pathway::Union{Pathway, Nothing}) where {K<:KeyType}
    (x_prime, x) = prepare_kernel(b,a, pathway)
    (x, x_prime)
end

# RNASeqCall vs. NAViewType
function prepare_kernel(a::DataView{K, RNASeqCall}, b::DataView{K, NAViewType}, pathway::Union{Pathway{K}, Nothing}) where {K<:KeyType}
    keys_to_use = get_common_keys(a, b, pathway)
    (map(m -> a.measurements[m].normalized_value, keys_to_use), [falses(length(keys_to_use))...])
end

# general view type vs. general view type
function prepare_kernel(a::DataView{K, V}, b::DataView{K, V}, pathway::Union{Pathway, Nothing}) where {K<:KeyType, V <:ValueViewType}
    common_keys = get_common_keys(a, b, pathway)
    # hack: when no keys are shared, we have small similarity
    if length(common_keys) == 0
        x = [-2.]
        x_prime = [2.]
        return (x, x_prime  )
    end
    x = map(k -> a.measurements[k].normalized_value, common_keys)
    x_prime = map(k -> b.measurements[k].normalized_value, common_keys)
    (x, x_prime)
end

function prepare_kernel(a::DataView{K, ExomeSeq}, b::DataView{K, ExomeSeq}, pathway::Union{Pathway{K}, Nothing})  where {K<:KeyType}
    common_keys = get_common_keys(a, b, pathway)
    x = Vector{Float64}()
    x_prime = Vector{Float64}()
    factor_exact_matches = 1.
    for k in common_keys
        # check for match in protein change (this happens rarely and we want to emphasize this)
        # if there is match, make sure to include only the ExomeSeq measurement containing it
        measurements_a = a.measurements[k]
        measurements_b = b.measurements[k]

        pcs_a = map(es -> es.protein_change, measurements_a)
        pcs_b = map(es -> es.protein_change, measurements_b)

        protein_change_intersection = intersect(pcs_a, pcs_b)
        if length(protein_change_intersection) > 0
            for pc in protein_change_intersection
                push!(x, get_measurement_value(filter(m -> m.protein_change == pc, measurements_a))...)
                push!(x_prime, get_measurement_value(filter(m -> m.protein_change == pc, measurements_a))...)
                if pc != "None" factor_exact_matches += 1 end
            end
        else
            push!(x, mean(get_measurement_value(a, k)))
            push!(x_prime, mean(get_measurement_value(b, k)))
        end

    end
    if length(common_keys) == 0
        x = [-2.]
        x_prime = [2.]
    end
    (x, x_prime, Float64(length(x) * factor_exact_matches))
end

function prepare_kernel(a::DataView{K, ExomeSeq}, b::DataView{K, NAViewType}, pathway::Union{Pathway{K}, Nothing}) where {K<:KeyType}
    x = map(k -> mean(get_measurement_value(a, k)), [a.used_keys...])
    x_prime = zeros(length(x))
    (x, x_prime)
end

prepare_kernel(a::DataView{K,V}, b::DataView{K,V}, pathway::Union{Pathway, Nothing}) where {K <: KeyType, V <: NAViewType} = ([randn()],[randn()])

function get_common_keys(a::DataView, b::DataView, pathway::Union{Pathway, Nothing})
    # @info "getting common keys" pathway
    if pathway != nothing
        return extract_pathway_keys(a, b, pathway)
    end
    if (b.view_type == NAViewType) || (b.used_keys == a.used_keys)
        return collect(a.used_keys)
    elseif haskey(a.common_keys, b.cell_line_id)
        return a.common_keys[b.cell_line_id]
    elseif haskey(b.common_keys, a.cell_line_id)
        return b.common_keys[a.cell_line_id]
    end
    common_keys = collect(intersect(a.used_keys, b.used_keys))
    a.common_keys[b.cell_line_id] = common_keys
    b.common_keys[a.cell_line_id] = common_keys
    common_keys
end

function extract_pathway_keys(a::DataView, b::DataView, pathway::Pathway)
    # if b.view_type == NAViewType
    #     return intersect(pathway.genes, a.used_keys)
    # end
    # if haskey(a.common_keys, b.cell_line_id)
    #     # @info "found precomputed key intersection"
    #     return intersect(pathway.genes, a.common_keys[b.cell_line_id])
    # end
    # if haskey(b.common_keys, a.cell_line_id)
    #     # @info "found precomputed key intersection"
    #     return intersect(pathway.genes, b.common_keys[a.cell_line_id])
    # end
    common_keys = get_common_keys(a,b,nothing)
    intersect(pathway.genes, common_keys)
end

function extract_pathway_keys(a::DataView{<:KeyType, RPPA}, b::DataView{<:KeyType, RPPA}, pathway::Pathway)
    common_keys = get_common_keys(a,b,nothing)

    keys_to_use = Vector{KeyType}()
    pathway_hgnc_ids = map(k -> k.hgnc_id, pathway.genes)
    hgnc_ids_in_common_keys = map(g -> g.hgnc_id, common_keys)

    hgnc_intersection = intersect(pathway_hgnc_ids, hgnc_ids_in_common_keys)
    key_idx = findall((in)(hgnc_intersection), hgnc_ids_in_common_keys)
    common_keys[key_idx]
end

#actual kernel functions
# Gaussian
kernel_function(x::Vector{Float64}, x_prime::Vector{Float64}, l::Float64 = Float64(length(x))) = exp(-dot(x-x_prime, x-x_prime)/l)
# Jaccard
kernel_function(x::Vector{Bool}, x_prime::Vector{Bool}) = x'*x_prime/(x'*x + x_prime'*x_prime - x'*x_prime)


text:
include("LQSolver.jl")
include("subproblem.jl")
include("AbstractProgressiveHedging.jl")


text:
newref(x) = K_Ref(K_new(x))
@testset "K_new" begin
  @test (x = newref(Int64(1)); (xt(x.x), xj(x.x)) == (-KJ, 1))
  @test (x = newref(Date(2000)); (xt(x.x), xi(x.x)) == (-KD, 0))
  @test (x = newref(DateTime(2000)); (xt(x.x), xj(x.x)) == (-KP, 0))
  @test (x = newref(Dates.Microsecond(2)); (xt(x.x), xj(x.x)) == (-KN, 2000))
end
struct TestType end
@testset "ktypecode(::Type) tests" begin
  @test_throws ErrorException ktypecode(TestType)
  @test ktypecode(Bool) == KB
  @test ktypecode(UInt128) == UU
  @test ktypecode(UInt8) == KG
  @test ktypecode(Int16) == KH
  @test ktypecode(Int32) == KI
  @test ktypecode(Int64) == KJ
  @test ktypecode(Float32) == KE
  @test ktypecode(Float64) == KF
  @test ktypecode(Char) == KC
  @test ktypecode(Symbol) == KS
  @test ktypecode(Q.TimeStamp) == KP
  @test ktypecode(Q.Month) == KM
  @test ktypecode(Q.Date) == KD
  @test ktypecode(Q.DateTimeF) == KZ
  @test ktypecode(Q.TimeSpan) == KN
  @test ktypecode(Q.Minute) == KU
  @test ktypecode(Q.Second) == KV
  @test ktypecode(Q.Time) == KT
end


text:
module TestInverseCholeskyDot
using LinearAlgebra
using LazyInverses
using LazyInverses: diag_mul
using Test

# testing particularly efficient algebraic operations with inverse(Cholesky)
@testset "cholesky" begin
    n, k = 128, 128 # TODO: also include test for n > 1024 to test multi-threaded ternary dot
    element_types = (Float32, Float64, ComplexF32, ComplexF64)
    for elty in element_types
        @testset "eltype $elty" begin
            x, y = rand(elty, n), rand(elty, n)
            X, Y = rand(elty, k, n), rand(elty, n, k)
            A = randn(elty, n, n)
            A = A'A + I # for numerical stability, especially for low-precision floats
            # with Cholesky
            C = cholesky(A)
            invC = Inverse(C)
            @test dot(x, invC, x) ≈ dot(x, C\x)
            @test dot(x, invC, y) ≈ dot(x, C\y)
            @test *(X, invC, X') ≈ *(X, C\X')
            @test *(X, invC, Y) ≈ *(X, C\Y)
            @test diag_mul(X, invC, X') ≈ diag(*(X, invC, X'))
            @test diag_mul(X, invC, Y) ≈ diag(*(X, invC, Y))

            inv_A = inv(A)
            diag_inv_A = diag(inv_A)
            @test diag(invC) ≈ diag_inv_A

            # with pivoted Cholesky
            C = cholesky(A, Val(true))
            invC = Inverse(C)
            @test dot(x, invC, x) ≈ dot(x, C\x)
            @test dot(x, invC, y) ≈ dot(x, C\y)
            @test *(X, invC, X') ≈ *(X, C\X')
            @test *(X, invC, Y) ≈ *(X, C\Y)
            @test diag(invC) ≈ diag_inv_A
            @test diag_mul(X, invC, X') ≈ diag(*(X, invC, X'))
            @test diag_mul(X, invC, Y) ≈ diag(*(X, invC, Y))
        end # testset elty
    end # loop elty
end # testset cholesky
end # TestInverseCholeskyDot


text:
# Defining the size of the word vectors
println("== Defining word vectors size...")
m = 350

# Important packages
println("== Compiling useful packages...")
using JLD
using JSON

# Importing the words file
println("== Importing dictionary...")
f = open("preprocessed_data/words.txt")
raw_words = readlines(f)
close(f)

# Preprocessing words
println("== Preprocessing dictionary...")
words = Array(AbstractString, 0)
for raw_word in raw_words
	push!(words,strip(raw_word, ['\n']))
end

# Importing the sents file
println("== Importing sentences...")
f = open("preprocessed_data/sents.txt")
sents = readlines(f)
close(f)

# Checking array sizes
println("- Dictionary size: ", length(words))
println("- Number of sentences: ", length(sents))

# Initializing the count array
println("== Initializing the count array...")
count = zeros(Int32, length(words), length(words))
occurencies = zeros(Int32, length(words))

# Counting word matches (symmetric)
println("== Counting...")
w = 0
for sent in sents
	for i = 1:length(words)
        if (contains(sent, words[i]))
            occurencies[i] += 1
			for j = i:length(words)
                if (contains(sent, words[j]))
                    occurencies[j] += 1
					count[i,j] += 1
					count[j,i] += 1
				end
			end
		end
	end
	w+=1
	println(w, "/", length(sents), " sents processed")
end

# Saving the array
println("== Exporting count array...")
save("backup_matrices/counts.jld", "counts", count)
save("backup_matrices/occurencies.jld", "occurencies", occurencies)

# Checks if the count array is symmetric
isSymmetric = Symmetric(count) == count
println("== The count array is symmetric: ", isSymmetric)
if (isSymmetric == false)
	exit
end

# Saving the normalized array
println("== Exporting normalized array...")
save("backup_matrices/normalized.jld", "norm", count)

# Decomposing w/ eigendecomposition
println("== Eigendecomposition (ETA: 50min from start)...")
d, v = eig(count)

# Saving matrices
println("== Exporting eigenvalues and eigenvectors arrays...")
save("backup_matrices/eigenvalues.jld", "values", d)
save("backup_matrices/eigenvectors.jld", "vectors", v)

# Sorting eigenvalues by its absolute value (BubbleSort)
println("== BubbleSorting eigenvalues by its absolute values...")
for i = 1:length(d)
    swapped = false
    for j = 1:length(d) - i - 1
        if (abs(d[j]) < abs(d[j + 1]))
            aux = d[j]
            d[j] = d[j+1]
            d[j + 1] = aux
            swapped = true
        end
    end
    if swapped == false
        break
    end
end

# Filling a trunkated matrix (MxN) w/ most significant eigenvalues
println("== Filling the trunkated matrix (MxN)...")
d2 = zeros(m, length(d))
for i = 1:m
    d2[i, i] += sqrt(d[i])
end

# Constructing W matrix with all the word vectors
println("== Constructing W matrix with all the word vectors...")
word_vectors = d2 * v.'

# Exporting word vectors matrix
println("== Exporting word vectors matrix...")
save("backup_matrices/wordvectors.jld", "wordvectors", word_vectors)

# Generating dictionary
println("== Generating dictionary...")
dict = [words[i]=>word_vectors.'[i,:] for i = 1:length(words)]
occur = [words[i]=>occurencies[i] for i = 1:length(words)]

# Exporting dictionary as JSON
println("== Exporting dictionary as JSON...")
json_string = JSON.json(dict)
json_occur = JSON.json(occur)
open("output_json/w2v_dict.json", "w") do f
    write(f, json_string)
end
open("output_json/count_dict.json", "w") do g
    write(g, json_occur)
end

# Prints successful ending
println("== Congratulations! This was successful!")


text:
# Some code in `src/problem_depot` was modified from MathOptInterface
# which is available under an MIT license (see LICENSE).

module ProblemDepot
using BenchmarkTools, Test

using Convex
using LinearAlgebra
using LinearAlgebra: eigen, I, opnorm

randperm(d) = sortperm(rand(d))
shuffle(x) = x[randperm(length(x))]
mean(x) = sum(x) / length(x)
eye(n, T) = Matrix{T}(I, n, n)
eye(n) = Matrix{Float64}(I, n, n)

"""
    const PROBLEMS = Dict{String, Dict{String, Function}}()

A "depot" of Convex.jl problems, subdivided into categories.
Each problem is stored as a function with the signature

    f(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}

where `handle_problem!` specifies what to do with the `Problem` instance
(e.g., `solve!` it with a chosen solver), an option `test` to choose
whether or not to test the values (assuming it has been solved),
tolerances for the tests, and a numeric type in which the problem
should be specified (currently, this is not respected and all
problems are specified in `Float64` precision).

See also [`run_tests`](@ref) and [`benchmark_suite`](@ref) for helpers
to use these problems in testing or benchmarking.

### Examples

```julia
julia> PROBLEMS["affine"]["affine_diag_atom"]
affine_diag_atom (generic function with 1 method)
```
"""
const PROBLEMS = Dict{String, Dict{String, Function}}()

"""
    foreach_problem(apply::Function, [class::String],
        problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing;
        exclude::Vector{Regex} = Regex[])

Provides a convience method for iterating over problems in [`PROBLEMS`](@ref).
For each problem in [`PROBLEMS`](@ref), apply the function `apply`, which
takes two arguments: the name of the function associated to the problem,
and the function associated to the problem itself.

Optionally, pass a second argument `class` to only iterate over a class of
problems (`class` should satsify `class ∈ keys(PROBLEMS)`), and pass third
argument `problems` to only allow certain problems (specified by exact names or
regex). Use the `exclude` keyword argument to exclude problems by regex.
"""
function foreach_problem(   apply::Function,
                            problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing;
                            exclude::Vector{Regex} = Regex[])
    for class in keys(PROBLEMS)
        any(occursin.(exclude, Ref(class))) && continue
        foreach_problem(apply, class, problems; exclude=exclude)
    end
end

function foreach_problem(   apply::Function,
                            class::String,
                            problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing;
                            exclude::Vector{Regex} = Regex[])
    for (name, func) in PROBLEMS[class]
        any(occursin.(exclude, Ref(name))) && continue
        if problems !== nothing
            problems isa Vector{String} && !(name ∈ problems) && continue
            problems isa Vector{Regex} && !any(occursin.(problems, Ref(name))) && continue
        end
        apply(name, func)
    end
end

"""
    run_tests(
        handle_problem!::Function;
        problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing;
        exclude::Vector{Regex} = Regex[],
        T=Float64, atol=1e-3, rtol=0.0,
    )

Run a set of tests. `handle_problem!` should be a function that takes one
argument, a Convex.jl `Problem` and processes it (e.g. `solve!` the problem with
a specific solver).

Use `exclude` to exclude a subset of sets; automatically excludes
`r"benchmark"`. Optionally, pass a second argument `problems` to only allow certain problems
(specified by exact names or regex). The test tolerances specified by `atol` and
`rtol`. Set `T` to choose a numeric type for the problem. Currently
this is only used for choosing the type parameter of the underlying
MathOptInterface model, but not for the actual problem data.

### Examples

```julia
run_tests(exclude=[r"mip"]) do p
    solve!(p, SCSSolver(verbose=0))
end
```
"""
function run_tests( handle_problem!::Function,
                    problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing;
                    exclude::Vector{Regex} = Regex[], T=Float64, atol=1e-3, rtol=0.0)
    push!(exclude, r"benchmark")
    for class in keys(PROBLEMS)
        any(occursin.(exclude, Ref(class))) && continue
        @testset "$class" begin
            foreach_problem(class, problems; exclude=exclude) do name, problem_func
                @testset "$name" begin
                    problem_func(handle_problem!, Val(true), atol, rtol, T)
                end
            end
        end
    end
end

"""
    benchmark_suite(
        handle_problem!::Function,
        problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing;
        exclude::Vector{Regex} = Regex[],
        test = Val(false),
        T=Float64, atol=1e-3, rtol=0.0,
    )

Create a benchmark_suite of benchmarks. `handle_problem!` should be a function
that takes one argument, a Convex.jl `Problem` and processes it (e.g. `solve!`
the problem with a specific solver). Pass a second argument `problems` to specify
run benchmarks only with certain problems (specified by exact names or regex).

Use `exclude` to exclude a subset of benchmarks. Optionally, pass a second
argument `problems` to only allow certain problems (specified by exact names or
regex). Set `test=true` to also check the answers, with tolerances specified by
`atol` and `rtol`. Set `T` to choose a numeric type for the problem. Currently
this is only used for choosing the type parameter of the underlying
MathOptInterface model, but not for the actual problem data.

### Examples

```julia
benchmark_suite(exclude=[r"mip"]) do p
    solve!(p, SCSSolver(verbose=0))
end
```
"""
function benchmark_suite(handle_problem!::Function,
                         problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing;
                         exclude::Vector{Regex} = Regex[],
                         T=Float64, atol=1e-3, rtol=0.0, test = Val(false))
    group = BenchmarkGroup()
    for class in keys(PROBLEMS)
        any(occursin.(exclude, Ref(class))) && continue
        group[class] = BenchmarkGroup()
        foreach_problem(class, problems; exclude=exclude) do name, problem_func
            group[class][name] = @benchmarkable $problem_func($handle_problem!, $test, $atol, $rtol, $T)
        end
    end
    return group
end

macro add_problem(prefix, q)
    @assert prefix isa Symbol
    if q.head == :block
        f = q.args[2]
    elseif q.head == :function
        f = q
    else
        error("head $(q.head) unexpected")
    end
    name = f.args[1].args[1]
    if name isa Expr
        name = name.args[1]
    end
    return quote
        $(esc(f))
        dict = get!(PROBLEMS, String($(Base.Meta.quot(prefix))), Dict{String,Function}())
        dict[String($(Base.Meta.quot(name)))] = $(esc(name))
    end
end

include("problems/affine.jl")
include("problems/constant.jl")
include("problems/exp.jl")
include("problems/lp.jl")
include("problems/mip.jl")
include("problems/sdp_and_exp.jl")
include("problems/sdp.jl")
include("problems/socp.jl")

end


text:
using COSMO
using Random, Test, LinearAlgebra

verbosity = 0 # Change to 2 to see LOBPCG's output
rng = MersenneTwister(1)
λ = [1:6; -50.00:-0.00]
n = length(λ)
Q = Matrix(qr(randn(rng, n, n)).Q)
X = Q*Diagonal(λ)*Q'

n_ = Int(n*(n + 1)/2)

function exact_projection(A::Matrix)
    l, V = eigen(Symmetric(A))
    return Symmetric(V*Diagonal(max.(l, 0))*V')
end

function lobpcg_unit_test!(cone, X::Matrix)
    x = Main.COSMO.extract_upper_triangle(X)
    X_proj_exact = exact_projection(X)
    COSMO.project!(x, cone)
    X_proj = Main.COSMO.populate_upper_triangle(x)

    @test norm(X_proj - X_proj_exact) < 1e-5 # This is tolerance is hand-wavy
end

function lobpcg_unit_tests(X::Matrix)
    cone = Main.COSMO.PsdConeTriangleLOBPCG(n_, verbosity=verbosity, max_iter=200)
    @testset "Initial call (exact eigendecomposition)" begin
        lobpcg_unit_test!(cone, X)
        @test cone.exact_projections == 1
        @test cone.lobpcg_iterations == 0
    end

    @testset "Warm starting - already solved to tolerance" begin
        cone.iteration = 10^3
        lobpcg_unit_test!(cone, X + 1e-6*Symmetric(randn(rng, n, n)))
        @test cone.exact_projections == 1
        @test cone.lobpcg_iterations == 1
    end
    @testset "Warm starting - a few iterations required" begin
        cone.iteration = 10^8
        lobpcg_unit_test!(cone, X + 1e-6*Symmetric(randn(rng, n, n)))
        @test cone.exact_projections == 1
        @test cone.lobpcg_iterations >= 3
        @test cone.lobpcg_iterations < 10
    end

    @testset "Warm starting - expand subspace" begin
        cone.iteration = 10^8
        if cone.is_subspace_positive
            lobpcg_unit_test!(cone, X + 3.1*I)
        else
            lobpcg_unit_test!(cone, X - 3.1*I)
        end
        @test size(cone.U, 2) > 6
        @test cone.exact_projections == 1
        @test cone.lobpcg_iterations >= 3
        @test cone.lobpcg_iterations < 50
    end

    @testset "Warm starting - contract subspace" begin
        cone.iteration = 10^8
        if cone.is_subspace_positive
            lobpcg_unit_test!(cone, X - 1.1*I)
        else
            lobpcg_unit_test!(cone, X + 1.1*I)
        end
        @test size(cone.U, 2) < 6
        @test cone.exact_projections == 1
        @test cone.lobpcg_iterations >= 3
        @test cone.lobpcg_iterations < 60
    end
end

@testset "LOBPCG PSD Projection" begin
    @testset "Track Positive Eigenspace" begin
        lobpcg_unit_tests(X)
    end
    @testset "Track Negative Eigenspace" begin
        lobpcg_unit_tests(-X)
    end
end

nothing