text:
using InterpretableValidation
using Test
using Distributions
using Random
using ExprRules
using POMDPs
using POMDPSimulators
using POMDPModels
using POMDPPolicies

rng = Random.GLOBAL_RNG
x = collect(1:100.)
N = length(x)
uniform = ConstrainedTimeseriesDistribution(IID(x, Uniform(2,3)))
normal = ConstrainedTimeseriesDistribution(IID(x, Normal(0.5, 1.)))
categorical = ConstrainedTimeseriesDistribution(IID(x, Categorical([0.1, 0.3, 0.2, 0.2, 0.2])))
gp = ConstrainedTimeseriesDistribution(GaussianProcess(m = (x) -> 0, k = squared_exp_kernel(l=4), x = x))
mvts = MvTimeseriesDistribution(:x => uniform, :y => normal, :z => categorical, :gp =>gp)

# Test the contents of the default range
comparison_distribution = default_comparison_distribution(mvts)
@test comparison_distribution[:x] == Uniform(2,3)
@test comparison_distribution[:y] == Uniform(0.5-5,0.5+5)
@test comparison_distribution[:z] == Categorical(5)
@test comparison_distribution[:gp] isa Uniform
@test comparison_distribution[:gp].a < -2
@test comparison_distribution[:gp].b > 2

comparisons = default_comparisons(mvts)
@test comparisons[:x] == [Symbol(".<="), Symbol(".>=")]
@test comparisons[:z] == [Symbol(".==")]

# Sample some random comparisons
for i=1:100
    expr = sample_comparison(comparison_distribution, comparisons, Random.GLOBAL_RNG)
    @test expr.head == :call
    sym = expr.args[2]
    if sym == :z
        @test expr.args[1] == Symbol(".==")
    else
        @test expr.args[1] in [Symbol(".=="),Symbol(".<="), Symbol(".>=")]
        @test expr.args[3] >= comparison_distribution[sym].a
        @test expr.args[3] <= comparison_distribution[sym].b
    end
    @test sym in keys(comparison_distribution)
end

# Test the grammar rules and grammar sampling
g = create_stl_grammar(length(mvts), comparison_distribution, comparisons)
@test g.rules[1] == Meta.parse("R && R")
# @test g.rules[5] == Meta.parse("any_between(τ, C, C)")
# @test g.iseval[7]
# @test g.iseval[10]

# Test the loss function
set_global_grammar_params(N, comparison_distribution, comparisons, rng)
ev(t::Dict{Symbol, Array}) = sum(sum.(values(t)))
ev(rand(mvts))
rn = rand(RuleNode, g, :R, 30)
loss = loss_fn(ev, mvts)
@test loss(rn, g) isa Real

results = optimize_timed_stl((x) -> rand(), mvts, Npop=10, Niter=3, verbose = false)
@test results.loss > 0

# Sample until we have a feasible rule
max_tries, trial = 1000, 1
while true
    trial > max_tries && break
    global rn = rand(RuleNode, g, :R, 3)
    constraints = sample_constraints(get_executable(rn, g), N_pts(mvts), Random.GLOBAL_RNG)
    mvts2 = deepcopy(mvts)
    constrain_timeseries!(mvts2, constraints)
    isfeasible(mvts2) && break
end
ex = get_executable(rn, g)
constraints = sample_constraints(ex, N_pts(mvts), Random.GLOBAL_RNG)
mvts2 = deepcopy(mvts)
constrain_timeseries!(mvts2, constraints)
@test isfeasible(mvts2)

loss = loss_fn((x) -> throw(error("itsaerror")), mvts)
@test_throws ErrorException loss(rn, g)

loss = loss_fn((x) -> throw(InfeasibleConstraint("itsaerror")), mvts)
@test loss(rn, g) >= 1e9

mdp = SimpleGridWorld(tprob = 1)
hist = simulate(HistoryRecorder(), mdp, RandomPolicy(mdp), GWPos(3,3))

## Test the playback simulator
playback = PlaybackPolicy(collect(action_hist(hist)), RandomPolicy(mdp))
@test all(playback.actions .== action_hist(hist))
@test playback.backup_policy isa RandomPolicy
@test playback.i == 1

hist2 = simulate(HistoryRecorder(), mdp, playback, GWPos(3,3))
@test hist == hist2

action(playback, GWPos(3,3))

t, fn = discrete_action_mdp(mdp, 1000, use_prob = false)
@test t[:a].timeseries_distribution.distribution.p == 0.25*ones(4)
@test length(t) == 1
@test collect(keys(t))[1] == :a
@test all([t[:a].feasible[i] == [1,2,3,4] for i=1:100])
d = rand(t)
lnp = logpdf(t, d)
fn(d)

aprob = Dict(:up =>0.1, :down => 0.1, :left =>0.1, :right=>0.7)
t, fn = discrete_action_mdp(mdp, 1000, use_prob = true, action_probability = (mdp, a) -> aprob[a])
@test t[:a].timeseries_distribution.distribution.p == [0.1, 0.1, 0.1, 0.7]

h = sample_history(Meta.parse("any(a .== 2)"), t, mdp)
@test h isa SimHistory


text:
# Autogenerated wrapper script for TensorRT_jll for x86_64-linux-gnu-cuda+11.2
export libnvinfer, libnvinfer_plugin, libnvonnxparser, libnvparsers, trtexec

using CUDNN_jll
JLLWrappers.@generate_wrapper_header("TensorRT")
JLLWrappers.@declare_library_product(libnvinfer, "libnvinfer.so.8")
JLLWrappers.@declare_library_product(libnvinfer_plugin, "libnvinfer_plugin.so.8")
JLLWrappers.@declare_library_product(libnvonnxparser, "libnvonnxparser.so.8")
JLLWrappers.@declare_library_product(libnvparsers, "libnvparsers.so.8")
JLLWrappers.@declare_executable_product(trtexec)
function __init__()
    JLLWrappers.@generate_init_header(CUDNN_jll)
    JLLWrappers.@init_library_product(
        libnvinfer,
        "lib/libnvinfer.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libnvinfer_plugin,
        "lib/libnvinfer_plugin.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libnvonnxparser,
        "lib/libnvonnxparser.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libnvparsers,
        "lib/libnvparsers.so",
        nothing,
    )

    JLLWrappers.@init_executable_product(
        trtexec,
        "bin/trtexec",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# Generic Parser Helper
# =====================

"""
    EP.Named{:MyTag}(parser)

Construct an ExprParser with a type identified by `MyTag`. This is helpful if you
have multiple versions of a similar parser and would like to easily distinguish them during dispatch.

The `Named{:MyTag}` wrapper is also passed on to the parsed value.

Works with any ExprParser, also custom defined ones.

# Examples
```jldoctest
julia> using ExprParsers

julia> parser1 = EP.Named{:simple}(EP.Assignment(left = EP.anysymbol));

julia> parser2 = EP.Named{:any}(EP.Assignment());

julia> parse_expr(parser1, :(a = 4))
ExprParsers.Named{:simple,ExprParsers.Assignment_Parsed}(EP.Assignment_Parsed(left=:a, right=4))
julia> parse_expr(parser1, :(a.b = 4))
ERROR: ParseError: Expected type `Symbol`, got `a.b` of type `Expr`.
julia> parse_expr(parser2, :(a = 4))
ExprParsers.Named{:any,ExprParsers.Assignment_Parsed}(EP.Assignment_Parsed(left=:a, right=4))
julia> parse_expr(parser2, :(a.b = 4))
ExprParsers.Named{:any,ExprParsers.Assignment_Parsed}(EP.Assignment_Parsed(left=:(a.b), right=4))
```
"""
@def_structequal struct Named{Name, T} <: ExprParser
  value::T
  Named{Name}(value::T) where {Name, T} = new{Name, T}(value)
end

function parse_expr(parser::Named{Name}, expr) where Name
  parsed = parse_expr(parser.value, expr)
  Named{Name}(parsed)
end
to_expr(parsed::Named) = to_expr(parsed.value)

Base.getindex(named::Named) = named.value

@doc raw"""
    EP.Indexed(func_expecting_dict_as_only_argument)

Constructs an ExprParser where you can access dedicated subexpressions/subparsers via Dictionary lookup.
Most importantly, the shortcuts are preserved during `parse_expr()`.

Works with any ExprParser, also custom defined ones.

Concretely, here a toy example
```julia
EP.Indexed() do dict
  EP.Expr(quote
    a = $(dict[:a] = EP.Isa(Int))
    b = $(dict[:b] = EP.anysymbol)
  end)
end
```
As you can see, `EP.Indexed` is expecting a function which takes a `dict` as the only argument.
It best used with do-notation. The function then needs to return an `ExprParser`, but can do whatever it wants
in principle.
Shortcuts are now assigned by just using interpolation syntax `$(...)` and storing references to subparser
into the given `dict`. For example you see that `EP.Isa(Int)` is captured as `dict[:a]` before being used as a
subparser.

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.Indexed() do dict
         EP.Expr(quote
           a = $(dict[:a] = EP.Isa(Int))
           b = $(dict[:b] = EP.anysymbol)
         end)
       end;

julia> parser[:a]
ExprParsers.Isa{Int64}()
julia> parsed = parse_expr(parser, quote
         a = 42
         b = a
       end);

julia> parsed[:a], parsed[:b]
(42, :a)

julia> parse_expr(parser, quote
         a = 42
         b = :notasymbol
       end);
ERROR: ParseError: Expected type `Symbol`, got `:notasymbol` of type `QuoteNode`.
```
"""
@def_structequal struct Indexed{T} <: ExprParser
  _mapping::Dict
  _paths::Dict
  value::T
end
function Indexed(func_dict!)
  mapping = Dict()
  value = func_dict!(mapping)
  paths = Dict(k => _find_object(v, value) for (k,v) in mapping)
  Indexed(mapping, paths, value)
end

ProxyInterfaces.dict(parser::Indexed) = parser._mapping
ProxyInterfaces.dict(::Base.Type{Indexed{T}}) where T = T
ProxyInterfaces.@dict Indexed

function parse_expr(parser::Indexed{T}, expr) where T
  parsed = parse_expr(parser.value, expr)
  mapping = Dict(k => _get_object(parsed, nested_index) for (k, nested_index) in parser._paths)
  Indexed(mapping, parser._paths, parsed)
end
to_expr(parsed::Indexed) = to_expr(parsed.value)

function Base.setindex!(parsed::Indexed, value, key)
  nested_index = parsed._paths[key]
  # we need to update both the ._mapping as well as the .expr
  parsed._mapping[key] = value
  _set_object!(parsed.value, value, nested_index)
end

# Helper
# ------

# TODO these little helpers might be well replaced with some Lense Library
# but for now they are enough

"""
  return list of keys/symbols to access a reference or nothing if no is found
"""
function _find_object(obj, nested::Union{Vector, Tuple, Dict})
  # first check whether the object is any of the elements itself
  i = findfirst(nested) do x
    obj === x
  end
  issomething(i) && return i

  # recurse if nothing found
  # return first non-nothing value
  for (k, v) in zip(keys(nested), values(nested))
    subkeys = _find_object(obj, v)
    issomething(subkeys) && return [k; subkeys]
  end
  return nothing
end
_find_object(obj, nested) = _find_object(obj, Dict(key => getproperty(nested, key) for key in propertynames(nested)))

function _get_object(nested, nested_index, get_key)
  if length(nested_index) == 0
     # return everything for empty index
    nested
  else
    subobject = get_key(nested, nested_index[1])

    if length(nested_index) == 1
      subobject
    else
      _get_object(subobject, nested_index[2:end])
    end
  end
end
_get_object(nested::Union{Vector, Tuple, Dict}, nested_index) = _get_object(nested, nested_index, getindex)
_get_object(nested, nested_index) = _get_object(nested, nested_index, getproperty)

function _set_object!(nested, value, nested_index, get_key, set_key!)
  @assert length(nested_index) > 0 "we should never reach empty index in this recursion"
  if length(nested_index) == 1
    # use the last index to update the final object
    set_key!(nested, value, nested_index[1])
  else
    # otherwise recurse into subobject
    subobject = get_key(nested, nested_index[1])
    _set_object!(subobject, value, nested_index[2:end])
  end
  nested
end
function _set_object!(nested::Union{Vector, Tuple, Dict}, value, nested_index)
  _set_object!(nested, value, nested_index, getindex, setindex!)
  nested
end
function _set_object!(nested, value, nested_index)
  _set_object!(nested, value, nested_index, getproperty, (object, value, key) -> setproperty!(object, key, value))
  nested
end


text:
function test_padded_copy(backend::Backend, T)
  println("-- Testing CUDA PaddedCopy Utilities{$T}")
  eps = 1e-10
  width, height, channels, num = (4,5,6,7)
  pad = (2,3)

  orig = rand(T, width, height, channels, num)
  padded = round(100*rand(T, width+2*pad[1], height+2*pad[2], channels, num))

  orig_blob = make_blob(backend, orig)
  padded_blob = make_blob(backend, padded)
  dense2padded!(backend, padded_blob, orig_blob, pad)
  copy!(padded, padded_blob)

  got = padded[pad[1]+1:end-pad[1],pad[2]+1:end-pad[2],:,:]

  @test all(abs.(got - orig) .< eps)

  padded = rand(T, size(padded))
  copy!(padded_blob, padded)
  padded2dense!(backend, orig_blob, padded_blob, pad)
  copy!(orig, orig_blob)

  expected = padded[pad[1]+1:end-pad[1],pad[2]+1:end-pad[2],:,:]
  @test all(abs.(expected - orig) .< eps)
end

function test_padded_copy(backend::Backend)
  test_padded_copy(backend, Float32)
  test_padded_copy(backend, Float64)
end

test_padded_copy(backend_gpu)


text:
using Infinity
using StatsDiscretizations

ℝ_discr = RealLineDiscretizer{:open,:closed}(-5:0.1:5)
@test ℝ_discr(2) == Interval{:open,:closed}(1.9,2)
@test ℝ_discr(2.0001) == Interval{:open,:closed}(2,2.1)

@test_throws Exception ℝ_discr[0]
@test_throws Exception ℝ_discr[1000]

𝕀_discr = BoundedIntervalDiscretizer{:open, :closed}(-5:0.1:5)

@test collect(𝕀_discr) == collect(ℝ_discr[2:(end-1)])

xs = rand(1000)
discr_𝕀_xs = 𝕀_discr.(xs)
@test all( in.(xs, discr_𝕀_xs))
@test discr_𝕀_xs == 𝕀_discr.(discr_𝕀_xs)
@test 𝕀_discr(Interval{:open,:open}(0.05,0.06))  == Interval{:open,:closed}(0.0,0.1)
@test 𝕀_discr(Interval{:open,:closed}(0.099,0.1))  == Interval{:open,:closed}(0.0,0.1)

@test_throws ArgumentError 𝕀_discr(Interval{:open,:closed}(0.05,0.15))
@test_throws ArgumentError 𝕀_discr(Interval{:open,:closed}(0.05,0.15))

@test_throws Exception 𝕀_discr(-10)
@test_throws ArgumentError 𝕀_discr(+10)

discr_ℝ_xs = ℝ_discr.(xs)
@test all( in.(xs, discr_ℝ_xs))
@test discr_ℝ_xs == 𝕀_discr.(discr_ℝ_xs)
@test ℝ_discr(Interval{:open,:open}(0.05,0.06))  == Interval{:open,:closed}(0.0,0.1)
@test ℝ_discr(Interval{:open,:closed}(0.099,0.1))  == Interval{:open,:closed}(0.0,0.1)

@test_throws ArgumentError ℝ_discr(Interval{:open,:closed}(0.05,0.15))

@test ℝ_discr(-10) == Interval{:open,:closed}(-∞,-5.0)
@test ℝ_discr(10) == Interval{:open,:closed}(5.0, ∞)

finite_discretizer = FiniteSupportDiscretizer(-4:1:4)

@inferred finite_discretizer(2)
@test_throws String finite_discretizer(10)

double_fun = DiscretizedFunction(finite_discretizer, x->2x)

@test double_fun(4) == 8
@test double_fun(3) == 6

findfirst( ==(2), finite_discretizer)::Int

searchsorted(finite_discretizer, 2)


text:
using Formatting

function width{T<:Integer}(n::T)
    w = ndigits(n)
    n < 0 || return w
    return w + 1
end

function pretty{T<:Integer}(a::Array{T,2}, indent::Int=4)
    lo, hi = extrema(a)
    w = max(width(lo), width(hi))
    id = " "^indent
    fe = FormatExpr(@sprintf(" {:%dd}", w))
    s = id
    nrow = size(a)[1]
    for i in 1:nrow
        for j in a[i,:]
            s *= format(fe, j)
        end
        i != nrow || continue
        s *= "\n"*id
    end
    return s
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

function atest(model::MOI.ModelLike, config::MOIT.TestConfig{T}) where {T <: Real}
    @test config.atol == Base.rtoldefault(T)
    @test config.rtol == Base.rtoldefault(T)
    @test config.solve
    @test config.query
    @test config.duals
    @test config.infeas_certificates
end

function btest(model::MOI.ModelLike, config::MOIT.TestConfig)
    @test false # b is in exclude
end

const customtests = Dict("a" => atest,
                         "b" => btest)

MOIT.@moitestset custom

@testset "TestConfig" begin
    mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
    config = MOIT.TestConfig()
    customtest(mock, config, ["b"])
end


text:
export RandomSearch, RandomSearch_populate, RandomSearch_generation, populate, evaluate, generation

# CMA-ES implementation based on pycma: https://github.com/CMA-ES/pycma
# To install it: pip install cma

import Cambrian: populate, evaluate, generation
using PyCall
using Random

# With x random vector: new_genes = ax + b

mutable struct RandomSearch <: Cambrian.AbstractEvolution
    config::NamedTuple
    logger::CambrianLogger
    population::Array{AbstractESIndividual}
    elites::Array{AbstractESIndividual}
    fitness::Function
    gen::Int
    a::Array{Float64}
    b::Array{Float64}
end

function RandomSearch(cfg::NamedTuple, fitness::Function, a::Array{Float64}, b::Array{Float64}; T::Type=ESIndividual ,logfile=string("logs/", cfg.id, ".csv"))
    logger = CambrianLogger(logfile)
    elite = T(cfg)
    e = RandomSearch(cfg, logger, [], [], fitness, 0, a, b)
    e.config = merge((n_population=n_population = 4+3*ceil(Int, log(cfg.n_genes)), ), e.config)
    RandomSearch_populate(e, T)
    e.elites = deepcopy([e.population[i] for i in 1:cfg.n_elite])
    e
end

function RandomSearch(cfg::NamedTuple, fitness::Function, a::Float64, b::Float64; T::Type=ESIndividual ,logfile=string("logs/", cfg.id, ".csv"))
    arr_a = ones(cfg.n_genes) * a
    arr_b = ones(cfg.n_genes) * b
    RandomSearch(cfg, fitness, arr_a, arr_b; T=T, logfile=logfile)
end

function RandomSearch(cfg::NamedTuple, fitness::Function, a::Int64, b::Int64; T::Type=ESIndividual ,logfile=string("logs/", cfg.id, ".csv"))
    RandomSearch(cfg, fitness, a*1.0, b*1.0; T=T, logfile=logfile)
end

function RandomSearch(cfg::NamedTuple, fitness::Function; T::Type=ESIndividual ,logfile=string("logs/", cfg.id, ".csv"))
    RandomSearch(cfg, fitness, cfg.a, cfg.b; T=T, logfile=logfile)
end

"get next population from pycma and create Cambrian individuals"
function RandomSearch_populate(e::RandomSearch, T::Type)
    e.population=[]
    for i in 1:e.config.n_population
        genes = rand(Float64, e.config.n_genes) .* e.config.a .+ e.config.b
        indiv = T(genes, -Inf*ones(e.config.d_fitness))
        push!(e.population, indiv)
    end
    e
end

"get next population from pycma and create Cambrian individuals"
function RandomSearch_populate(e::RandomSearch)
    T = typeof(e.population[1])
    RandomSearch_populate(e, T)
end

"update pycma, called after populate and evaluate"
function RandomSearch_generation(e::RandomSearch)
    Cambrian.elites_generation(e)
end

populate(e::RandomSearch) = RandomSearch_populate(e)
evaluate(e::RandomSearch) = fitness_evaluate(e, e.fitness)
generation(e::RandomSearch) = RandomSearch_generation(e)


text:
# Example: Compartmental models in epidemiology
#=
- https://github.com/epirecipes/sir-julia/blob/master/markdown/function_map/function_map.md
- https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#Deterministic_versus_stochastic_epidemic_models
=#
using ModelingToolkit, Test

@inline function rate_to_proportion(r,t)
    1-exp(-r*t)
end;

# Independent and dependent variables and parameters
@parameters t c nsteps δt β γ
D = Difference(t; dt=0.1)
@variables S(t) I(t) R(t)

infection = rate_to_proportion(β*c*I/(S+I+R),δt)*S
recovery = rate_to_proportion(γ,δt)*I

# Equations
eqs = [D(S) ~ S-infection,
       D(I) ~ I+infection-recovery,
       D(R) ~ R+recovery]

# System
@named sys = DiscreteSystem(eqs,t,[S,I,R],[c,nsteps,δt,β,γ]; controls = [β, γ])
syss = structural_simplify(sys)
@test syss == syss

# Problem
u0 = [S => 990.0, I => 10.0, R => 0.0]
p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1, nsteps => 400]
tspan = (0.0,ModelingToolkit.value(substitute(nsteps,p))) # value function (from Symbolics) is used to convert a Num to Float64
prob_map = DiscreteProblem(sys,u0,tspan,p)

# Solution
using OrdinaryDiffEq
sol_map = solve(prob_map,FunctionMap());
@test sol_map[S] isa Vector

# Using defaults constructor
@parameters t c=10.0 nsteps=400 δt=0.1 β=0.05 γ=0.25
Diff = Difference(t; dt=0.1)
@variables S(t)=990.0 I(t)=10.0 R(t)=0.0

infection2 = rate_to_proportion(β*c*I/(S+I+R),δt)*S
recovery2 = rate_to_proportion(γ,δt)*I

eqs2 = [D(S) ~ S-infection2,
       D(I) ~ I+infection2-recovery2,
       D(R) ~ R+recovery2]

@named sys = DiscreteSystem(eqs2; controls = [β, γ])
@test ModelingToolkit.defaults(sys) != Dict()

prob_map2 = DiscreteProblem(sys,[],tspan)
sol_map2 = solve(prob_map,FunctionMap());

@test sol_map.u == sol_map2.u
@test sol_map.prob.p == sol_map2.prob.p

# Direct Implementation

function sir_map!(u_diff,u,p,t)
    (S,I,R) = u
    (β,c,γ,δt) = p
    N = S+I+R
    infection = rate_to_proportion(β*c*I/N,δt)*S
    recovery = rate_to_proportion(γ,δt)*I
    @inbounds begin
        u_diff[1] = S-infection
        u_diff[2] = I+infection-recovery
        u_diff[3] = R+recovery
    end
    nothing
end;
u0 = [990.0,10.0,0.0];
p = [0.05,10.0,0.25,0.1];
prob_map = DiscreteProblem(sir_map!,u0,tspan,p);
sol_map2 = solve(prob_map,FunctionMap());

@test Array(sol_map) ≈ Array(sol_map2)

# Delayed difference equation
@parameters t
@variables x(..) y(..) z(t)
D1 = Difference(t; dt=1.5)
D2 = Difference(t; dt=2)

@test ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(x(t-2)))
@test ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(y(t-1)))
@test !ModelingToolkit.is_delay_var(Symbolics.value(t), Symbolics.value(z))
@test_throws ErrorException ModelingToolkit.get_delay_val(Symbolics.value(t), Symbolics.arguments(Symbolics.value(x(t+2)))[1])
@test_throws ErrorException z(t)

# Equations
eqs = [
    D1(x(t)) ~ 0.4x(t) + 0.3x(t-1.5) + 0.1x(t-3),
    D2(y(t)) ~ 0.3y(t) + 0.7y(t-2) + 0.1z,
]

# System
@named sys = DiscreteSystem(eqs,t,[x(t),x(t-1.5),x(t-3),y(t),y(t-2),z],[])

eqs2, max_delay = ModelingToolkit.linearize_eqs(sys; return_max_delay=true)

@test max_delay[Symbolics.operation(Symbolics.value(x(t)))] ≈ 3
@test max_delay[Symbolics.operation(Symbolics.value(y(t)))] ≈ 2

linearized_eqs = [
    eqs
    x(t - 3.0) ~ x(t - 1.5)
    x(t - 1.5) ~ x(t)
    y(t - 2.0) ~ y(t)
]
@test all(eqs2 .== linearized_eqs)

# observed variable handling
@variables t x(t) RHS(t)
@parameters τ
@named fol = DiscreteSystem([D(x) ~ (1 - x)/τ]; observed=[RHS ~ (1 - x)/τ])
@test isequal(RHS, @nonamespace fol.RHS)
RHS2 = RHS
@unpack RHS = fol
@test isequal(RHS, RHS2)

@testset "Preface tests" begin
    @parameters t
    using OrdinaryDiffEq
    using Symbolics
    using DiffEqBase: isinplace
    using ModelingToolkit
    using SymbolicUtils.Code
    using SymbolicUtils: Sym

    c = [0]
    f =  function f(c, d::Vector{Float64}, u::Vector{Float64}, p, t::Float64, dt::Float64)
            c .= [c[1] + 1]
            d .= randn(length(u))
            nothing
        end

    dummy_identity(x, _) = x
    @register dummy_identity(x, y)

    u0 = ones(5)
    p0 = Float64[]
    syms = [Symbol(:a, i) for i in 1:5]
    syms_p = Symbol[]
    dt = 0.1
    @assert isinplace(f, 6)
    wf = let c=c, buffer = similar(u0), u=similar(u0), p=similar(p0), dt=dt
            t -> (f(c, buffer, u, p, t, dt); buffer)
    end

    num = hash(f) ⊻ length(u0) ⊻ length(p0)
    buffername = Symbol(:fmi_buffer_, num)

    Δ = DiscreteUpdate(t; dt=dt)
    us = map(s->(@variables $s(t))[1], syms)
    ps = map(s->(@variables $s(t))[1], syms_p)
    buffer, = @variables $buffername[1:length(u0)]
    dummy_var = Sym{Any}(:_) # this is safe because _ cannot be a rvalue in Julia

    ss = Iterators.flatten((us, ps))
    vv = Iterators.flatten((u0, p0))
    defs = Dict{Any, Any}(s=>v for (s, v) in zip(ss, vv))

    preface = [
               Assignment(dummy_var, SetArray(true, term(getfield, wf, Meta.quot(:u)), us))
               Assignment(dummy_var, SetArray(true, term(getfield, wf, Meta.quot(:p)), ps))
               Assignment(buffer, term(wf, t))
              ]
    eqs = map(1:length(us)) do i
        Δ(us[i]) ~ dummy_identity(buffer[i], us[i])
    end

    @named sys = DiscreteSystem(eqs, t, us, ps; defaults=defs, preface=preface)
    prob = DiscreteProblem(sys, [], (0.0, 1.0))
    sol = solve(prob, FunctionMap(); dt=dt)
    @test c[1]+1 == length(sol)
end

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
using InvertibleNetworks, LinearAlgebra, Test, Random

Random.seed!(11)

# Input
nx = 16
ny = 16
nz = 16
n_in = 2
n_hidden = 4
batchsize = 2
maxiter = 2

# Observed data
nxrec = 10
nyrec = 12
nt = 8
d = randn(Float32, nt, nxrec*nyrec, 1, batchsize)

# Modeling/imaging operator
J = 1 .+ rand(Float32, nt*nxrec*nyrec, nx*ny*nz)

# Link function
Ψ(η) = identity(η)

# Unrolled loop
L = NetworkLoop3D(n_in, n_hidden, maxiter, Ψ; type="HINT")

# Initializations
η = 10*randn(Float32, nx, ny, nz, 1, batchsize)
s = 10*randn(Float32, nx, ny, nz, n_in-1, batchsize)

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

# Test invertibility
η_, s_ = L.forward(η, s, d, J)
ηInv, sInv = L.inverse(η_, s_, d, J)
@test isapprox(norm(ηInv - η)/norm(η), 0f0, atol=1e-5)
@test isapprox(norm(sInv - s)/norm(s), 0f0, atol=1e-5)

# Test invertibility
η_, s_ = L.forward(η, s, d, J)
ηInv, sInv = L.backward(0f0.*η_, 0f0.*s_, η_, s_, d, J)[3:4]
@test isapprox(norm(ηInv - η)/norm(η), 0f0, atol=1e-5)
@test isapprox(norm(sInv - s)/norm(sInv), 0f0, atol=1e-5)

η_, s_ = L.inverse(η, s, d, J)
ηInv, sInv = L.forward(η_, s_, d, J)
@test isapprox(norm(ηInv - η)/norm(η), 0f0, atol=1e-5)
@test isapprox(norm(sInv - s)/norm(sInv), 0f0, atol=1e-5)

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

# Initializations
η = randn(Float32, nx, ny, nz, 1, batchsize)
s = randn(Float32, nx, ny, nz, n_in-1, batchsize)
η0 = randn(Float32, nx, ny, nz, 1, batchsize)
s0 = randn(Float32, nx, ny, nz, n_in-1, batchsize)
Δη = η - η0
Δs = s - s0

# Observed data
η_, s_ = L.forward(η, s, d, J)   # only need η

function loss(L, η0, s0, d, η, s)
    η_, s_ = L.forward(η0, s0, d, J)  # reshape
    Δη = η_ - η
    Δs = s_ - s    # no "observed" s, so Δs=0
    f = .5f0*norm(Δη)^2 + .5f0*norm(Δs)^2
    Δη_, Δs_ = L.backward(Δη, Δs, η_, s_, d, J)[1:2]
    return f, Δη_, Δs_, L.L[1].C.v1.grad, L.L[1].CL[1].RB.W1.grad
end

# Gradient test for input
f0, gη, gs = loss(L, η0, s0, d, η_, s_)[1:3]
h = 0.1f0
maxiter = 6
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test loop unrolling\n")
for j=1:maxiter
    f = loss(L, η0 + h*Δη, s0 + h*Δs, d, η_, s_)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(Δη, gη) - h*dot(Δs, gs))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for weights
L0 = NetworkLoop3D(n_in, n_hidden, maxiter, Ψ; type="HINT")
L_ini = deepcopy(L0)
dv = L.L[1].C.v1.data - L0.L[1].C.v1.data   # just test for 2 parameters
dW = L.L[1].CL[1].RB.W1.data - L0.L[1].CL[1].RB.W1.data
f0, gη, gs, gv, gW = loss(L0, η, s, d, η_, s_)
h = 0.05f0
maxiter = 5
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test loop unrolling\n")
for j=1:maxiter
    L0.L[1].C.v1.data = L_ini.L[1].C.v1.data + h*dv
    L0.L[1].CL[1].RB.W1.data = L_ini.L[1].CL[1].RB.W1.data + h*dW
    f = loss(L0, η, s, d, η_, s_)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dv, gv) - h*dot(dW, gW))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
#### This script simulates example data in the /data/ folder.
#### Phenotypes are stored in 6th column of the `.fam` files and also
#### stored separatedly in `phenotypes.txt`. Non-genetic covariates include
#### an intercept and sex, which is stored separately in `covariates.txt`.
#### The true β is stored in `normal_true_beta.txt`.

using Revise
using MendelIHT
using SnpArrays
using Random
using GLM
using DelimitedFiles

################################
######## Gaussian data #########
################################
n = 1000            # number of samples
p = 10000           # number of SNPs
k = 10              # 8 causal SNPs and 2 causal covariates (intercept + sex)
d = Normal          # Gaussian (continuous) phenotypes
l = IdentityLink()  # canonical link function

# set random seed
Random.seed!(1111)

# simulate `.bed` file with no missing data
x = simulate_random_snparray("normal.bed", n, p)
xla = SnpLinAlg{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true)

# nongenetic covariate: first column is the intercept, second column is sex: 0 = male 1 = female
z = ones(n, 2)
z[:, 2] .= rand(0:1, n)
writedlm("covariates.txt", z, ',') # save covariates
standardize!(@view(z[:, 2:end])) # standardize covariates

# randomly set genetic predictors where causal βᵢ ~ N(0, 1)
true_b = zeros(p)
true_b[1:k-2] = randn(k-2)
shuffle!(true_b)

# find correct position of genetic predictors
correct_position = findall(!iszero, true_b)

# save true SNP's position and effect size
open("normal_true_beta.txt", "w") do io
    println(io, "snpID,effectsize")
    for pos in correct_position
        println(io, "snp$pos,", true_b[pos])
    end
end

# define effect size of non-genetic predictors: intercept & sex
true_c = [1.0; 1.5]

# simulate phenotype using genetic and nongenetic predictors
# note for d = Normal, l = IdentityLink(), `prob = xla * true_b .+ z * true_c`
prob = GLM.linkinv.(l, xla * true_b .+ z * true_c)
y = [rand(d(i)) for i in prob]
y = Float64.(y); # turn y into floating point numbers
writedlm("phenotypes.txt", y) # save phenotypes

# create `.bim` and `.bam` files using phenotype
make_bim_fam_files(x, y, "normal")

####### IHT TEST RUNS #########

# without sex as covariate (but intercept automatically included)
result = iht("normal", 9)
[true_b[correct_position] result.beta[correct_position]] # compare IHT's result with answer

# include sex as covariate
result = iht("normal", "covariates.txt", 10)
[true_b[correct_position] result.beta[correct_position]] # compare IHT's result with answer

# phenotypes inputted as separate file
result = iht("phenotypes.txt", "normal", "covariates.txt", 10)
[true_b[correct_position] result.beta[correct_position]] # compare IHT's result with answer

##### CROSS VALIDATION TESTS ######

# without sex as covariate (but intercept automatically included)
mses = cross_validate("normal", 1:20)
argmin(mses)

# include sex as covariate
mses = cross_validate("normal", "covariates.txt", 1:20)
argmin(mses)

# phenotypes inputted as separate file
mses = cross_validate("phenotypes.txt", "normal", "covariates.txt", 1:20)
argmin(mses)


text:
############################################################
# judiLRWF #############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017
# Updated: March 2021, Mathias Louboutin (mlouboutin3@gatech.edu)

export judiLRWF, judiLRWFexception

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

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

mutable struct judiLRWFexception <: Exception
    msg :: String
end

############################################################
## Constructor
"""
    judiLRWF(info, geometry)
Low-rank wavefield operator which injects a wavelet q at every point of the subsurface. \\
`info` is an `Info` structure and `wavelet` is a cell array containing the wavelet(s).
Examples
========
`F` is a modeling operator of type `judiModeling` and `w` is a weighting matrix of type `judiWeights`:
    Pr = judiProjection(info, rec_geometry)
    Pw = judiLRWF(info, q.data)
    dobs = Pr*F*Pw'*w
    dw = Pw*F'*Pr'*dobs
"""
function judiLRWF(info::Info, data::Array{T, N}; DDT::DataType=Float32, RDT::DataType=DDT) where {T, N}
    (DDT == Float32 && RDT == Float32 && T == Float32) || throw(judiProjectionException("Domain and range types not supported"))
    m = info.n * info.nsrc
    n = info.n * sum(info.nt)
    wavelet = Array{Array{DDT, N}, 1}(undef, info.nsrc)
    for j=1:info.nsrc
        wavelet[j] = data
    end
    return judiLRWF{Float32,Float32}("restriction operator",m,n,info,wavelet)
end

function judiLRWF(info::Info, wavelet::Array{Array{T, N}, 1}; DDT::DataType=Float32, RDT::DataType=DDT) where {T, N}
    (DDT == Float32 && RDT == Float32 && T == Float32) || throw(judiProjectionException("Domain and range types not supported"))
    m = info.n * info.nsrc
    n = info.n * sum(info.nt)
    return judiLRWF{Float32,Float32}("restriction operator",m,n,info,wavelet)
end

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

# conj(judiProjection)
conj(A::judiLRWF{DDT,RDT}) where {DDT,RDT} =
    judiLRWF{DDT,RDT}("conj("*A.name*")",A.m,A.n,A.info,A.wavelet)

# transpose(judiProjection)
transpose(A::judiLRWF{DDT,RDT}) where {DDT,RDT} =
    judiLRWF{DDT,RDT}("injection operator",A.n,A.m,A.info,A.wavelet)

adjoint(A::judiLRWF{DDT,RDT}) where {DDT,RDT} =
    judiLRWF{DDT,RDT}("injection operator",A.n,A.m,A.info,A.wavelet)

############################################################
## overloaded Base *(...judiProjection...)

# *(judiLRWF, judiWeights)
function *(A::judiLRWF{ADDT,ARDT}, v::judiWeights{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v,1) || throw(judiLRWFexception("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiLRWF,judiWeights):",A.name,typeof(A),vDT]," / "))
    V = judiExtendedSource(A.info,A.wavelet,v.weights)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiLRWF,judiWeights):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(judiLRWF, vec)
function *(A::judiLRWF{ADDT,ARDT}, v::AbstractVector{vDT}) where {ADDT,ARDT,vDT}
    A.n == size(v, 1) || throw(judiLRWFexception("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiLRWF,AbstractVector):",A.name,typeof(A),vDT]," / "))
    V = judiExtendedSource(A.info,A.wavelet, v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiLRWF,AbstractVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

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

# Subsample Modeling operator
function subsample(P::judiLRWF{ADDT,ARDT}, srcnum) where {ADDT,ARDT}
    info = Info(P.info.n, length(srcnum), P.info.nt[srcnum])
    return judiLRWF(info, P.wavelet[srcnum];DDT=ADDT,RDT=ARDT)
end

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


text:
function makeprimalstrictlyfeasible!{T<:Number}(sdp::SparseSDP{T}, epsilon::Number)
    for (ci, v) in rhs(sdp)
        v = zero(T)
        blockmatrix = cons(sdp)[ci]
        for (_, block) in blockmatrix
            for (indices, entry) in block
                i, j = indices
                if i == j
                    v += entry
                end
            end
        end
        setrhs!(sdp, ci, rhs(sdp, ci) + epsilon*v)
    end
end

function makedualstrictlyfeasible!{T<:Number}(sdp::SparseSDP{T}, epsilon::Number)
    for (_, block) in obj(sdp)
        for (indices, value) in block
            i, j = indices
            if i == j
                block[i, j] = block[i,j] - epsilon
            end
        end
    end
end


text:
"""
    reformulate_disjunction(m::Model, disj...; bin_var, reformulation, param)

Reformulate disjunction on a JuMP model.

    reformulate_disjunction(disj, bin_var, reformulation, param)

Reformulate disjunction.
"""
function reformulate_disjunction(m::Model, disj...; bin_var, reformulation, param)
    #check disj
    disj = [check_constraint!(m, constr) for constr in disj]#check_disjunction!(m, disj)
    #get original variable refs and variable names
    vars = setdiff(all_variables(m), m[bin_var])
    var_names = unique(Symbol.([split("$var","[")[1] for var in vars]))
    if !in(:gdp_variable_refs, keys(object_dictionary(m)))
        @expression(m, gdp_variable_refs, vars)
    end
    if !in(:gdp_variable_names, keys(object_dictionary(m)))
        @expression(m, gdp_variable_names, var_names)
    end
    #run reformulation
    if reformulation == :hull
        disaggregate_variables(m, disj, bin_var)
        sum_disaggregated_variables(m, disj, bin_var)
    end
    reformulate_disjunction(disj, bin_var, reformulation, param)

    #show new constraints as a Dict
    new_constraints = Dict{Symbol,Any}(
        Symbol(bin_var,"[$i]") => disj[i] for i in eachindex(disj)
    )
    # new_constraints[Symbol(bin_var,"_XOR")] = constraint_by_name(m, "XOR(disj_$bin_var)")
    if reformulation == :hull
        for var in m[:gdp_variable_refs]
            agg_con_name = "$(var)_$(bin_var)_aggregation"
            new_constraints[Symbol(agg_con_name)] = constraint_by_name(m, agg_con_name)
        end
    end
    return new_constraints

    #remove model.optimize_hook ?

    # return m[bin_var]
end
function reformulate_disjunction(disj, bin_var, reformulation, param)
    for (i,constr) in enumerate(disj)
        reformulate_constraint(constr, bin_var, reformulation, param, i)
    end
end

"""
    reformulate_constraint(constr::Tuple, bin_var, reformulation, param, i)

Reformulate a Tuple of constraints.

    reformulate_constraint(constr::AbstractArray{<:ConstraintRef}, bin_var, reformulation, param, i, j = missing)

Reformulate a block of constraints.

    reformulate_constraint(constr::ConstraintRef, bin_var, reformulation, param, i, j = missing, k = missing)

Reformulate a constraint.
"""
function reformulate_constraint(constr::Tuple, bin_var, reformulation, param, i)
    for (j,constr_j) in enumerate(constr)
        reformulate_constraint(constr_j, bin_var, reformulation, param, i, j)
    end
end
function reformulate_constraint(constr::AbstractArray{<:ConstraintRef}, bin_var, reformulation, param, i, j = missing)
    for k in eachindex(constr)
        reformulate_constraint(constr[k], bin_var, reformulation, param, i, j, k)
    end
end
function reformulate_constraint(constr::ConstraintRef, bin_var, reformulation, param, i, j = missing, k = missing)
    if reformulation == :big_m
        big_m_reformulation!(constr, bin_var, param, i, j, k)
    elseif reformulation == :hull
        hull_reformulation!(constr, bin_var, param, i, j, k)
    end
end
reformulate_constraint(args...) = nothing


text:
# =========================================================================================
"""
    legendre{T <: AbstractFloat}(
        maxn::Int,
        x   ::T
    )

Precomputes the Legendre polynomials ``Pₙ(x)`` for ``n = 0, 1, ..., \\texttt{maxn}-1`` and
returns a generic function to access the values.

# Return type
`(generic function)`
```julia
(n::Int) -> Pₙ(x)   # return type: T
```

# Example

```jldoctest; setup = :(using NESSie.TestModel: legendre)
julia> p = legendre(3, 0.5);

julia> [p(n) for n in 0:2]    # [P₀(0.5), P₁(0.5), P₂(0.5)]
3-element Array{Float64,1}:
  1.0
  0.5
 -0.125
```
"""
function legendre(maxn::Int, x::T) where T <: AbstractFloat
    maxn == 0 && return n::Int -> T[][n+1]
    maxn == 1 && return n::Int -> [one(T)][n+1]

    P = Array{T}(undef, maxn)
    P[1] = one(T) # P₀(x) = 0
    P[2] = x      # P₁(x) = x

    # Use Bonnet's recursion formula for remaining values:
    # Pₙ₊₁(x) = 1/(n+1) ⋅ [(2n + 1)⋅x⋅Pₙ(x) - n⋅Pₙ₋₁}(x)]
    @inbounds for n in 1:maxn-2
        P[n + 2] = (n + 1) \ ((2n + 1) * x * P[n + 1] - n * P[n])
    end
    n::Int -> P[n+1]
end

# =========================================================================================
"""
    spherical_besseli{T <: AbstractFloat}(
        maxn::Int,
        r   ::T
    )

Precomputes the modified spherical Bessel function of the first kind ``iₙ(r)`` for
``n=-1, 0, ..., \\texttt{maxn}`` and returns a generic function to access the values.
``iₙ(r)`` is defined as

```math
iₙ(r) = \\sqrt{\\frac{π}{2r}} I_{n+0.5}(r),
```

where ``I_ν`` is the modified Bessel function of the first kind
[[Xie16]](@ref Bibliography).

# Return type
`(generic function)`
```julia
(n::Int) -> iₙ(r)   # return type: T
```
"""
function spherical_besseli(maxn::Int, r::T) where T <: AbstractFloat
    i = √(π / 2r) * [besseli(n + one(T)/2, r) for n in -1:maxn]
    n::Int -> i[n+2]
end

# =========================================================================================
"""
    spherical_besselk{T <: AbstractFloat}(
        maxn::Int,
        r   ::T
    )

Precomputes the modified spherical Bessel function of the second kind ``kₙ(r)`` for
``n=-1, 0, ..., \\texttt{maxn}`` and returns a generic function to access the values.
``kₙ(r)`` is defined as

```math
kₙ(r) = \\sqrt{\\frac{π}{2r}} K_{n+0.5}(r),
```

where ``K_ν`` is the modified Bessel function of the second kind
[[Xie16]](@ref Bibliography).

# Return type
`(generic function)`
```julia
(n::Int) -> kₙ(r)   # return type: T
```
"""
function spherical_besselk(maxn::Int, r::T) where T <: AbstractFloat
    k = √(π / 2r) * [besselk(n + one(T)/2, r) for n in -1:maxn]
    n::Int -> k[n+2]
end


text:
mutable struct ApplicationReference
    jurisdiction::String
    doc_number::String
    kind::String
    date::Date
    ApplicationReference() = new("", "", "", Date(9999))
end

mutable struct EarliestClaim
    date::Date
    EarliestClaim() = new(Date(9999))
end

mutable struct PriorityClaim
    earliest_claim::EarliestClaim
    PriorityClaim() = new(EarliestClaim())
end

mutable struct ExtractedName
    value::String
    ExtractedName() = new("")
end

mutable struct Applicant
    extracted_name::ExtractedName
    Applicant() = new(ExtractedName())
end

mutable struct Inventor
    extracted_name::ExtractedName
    Inventor() = new(ExtractedName())
end

mutable struct Party
    applicants::Array{Applicant, 1}
    inventors::Array{Inventor, 1}
    Party() = new(Applicant[], Inventor[])
end

mutable struct ClassificationsIpcr
    classifications::Array{Classification, 1}
    ClassificationsIpcr() = new(Classification[])
end

mutable struct ClassificationsCpc
    classifications::Array{Classification, 1}
    ClassificationsCpc() = new(Classification[])
end

mutable struct ClassificationsNational
    classifications::Array{Classification, 1}
    ClassificationsNational() = new(Classification[])
end

mutable struct DocumentId
    jurisdiction::String
    doc_number::String
    kind::String
    date::Date
    DocumentId() = new("", "", "", Date(9999))
end

mutable struct PatentCitation
    document_id::DocumentId
    lens_id::String
    PatentCitation() = new(DocumentId(), "")
end

mutable struct Citation
    patcit::Union{Nothing, PatentCitation}
    nplcit::Union{Nothing, NPLCitation}
    Citation() = new()
end

mutable struct Cited
    citations::Array{Citation, 1}
    npl_count::Int64
    patent_count::Int64
    Cited() = new(Citation[], 0, 0)
end

mutable struct CitingDoc
    document_id::DocumentId
    lens_id::String
    CitingDoc() = new(DocumentId(), "")
end

mutable struct CitedBy
    patents::Array{CitingDoc, 1}
    patent_count::Int64
    CitedBy() = new(CitingDoc[], 0)
end

mutable struct Biblio
    application_reference::ApplicationReference
    priority_claims::PriorityClaim
    invention_title::Array{Title, 1}
    parties::Party
    classifications_ipcr::ClassificationsIpcr
    classifications_cpc::ClassificationsCpc
    classifications_national::ClassificationsNational
    references_cited::Cited
    cited_by::CitedBy
    Biblio() = new(ApplicationReference(), PriorityClaim(), Title[], Party(),
                    ClassificationsIpcr(), ClassificationsCpc(), ClassificationsNational(), Cited(), CitedBy())
end

mutable struct Member
    document_id::DocumentId
    lens_id::String
    Member() = new(DocumentId(), "")
end

mutable struct SimpleFamily
    members::Array{Member, 1}
    size::Int64
    SimpleFamily() = new(Member[], 0)
end

mutable struct ExtendedFamily
    members::Array{Member, 1}
    size::Int64
    ExtendedFamily() = new(Member[], 0)
end

mutable struct LensFamily
    simple_family::SimpleFamily
    extended_family::ExtendedFamily
    LensFamily() = new(SimpleFamily(), ExtendedFamily())
end

mutable struct LegalStatus
    calculation_log::Vector{String}
    patent_status::String
    LegalStatus() = new(String[], "")
end

mutable struct ClaimText
    claim_text::Array{String, 1}
    ClaimText() = new(String[])
end

mutable struct Claim
    claims::Array{ClaimText, 1}
    lang::String
    Claim() = new(ClaimText[], "")
end

mutable struct LensApplication
    lens_id::String
    jurisdiction::String
    doc_number::String
    kind::String
    date_published::Date
    biblio::Biblio
    families::LensFamily
    legal_status::LegalStatus
    abstract::Array{Abstract, 1}
    claims::Array{Claim, 1}
    publication_type::String
    LensApplication() = new("", "", "", "", Date(9999), Biblio(), LensFamily(), LegalStatus(), Abstract[], Claim[], "")
end

StructTypes.StructType(::Type{ApplicationReference}) = StructTypes.Mutable()
StructTypes.StructType(::Type{EarliestClaim}) = StructTypes.Mutable()
StructTypes.StructType(::Type{PriorityClaim}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Title}) = StructTypes.Mutable()
StructTypes.StructType(::Type{ExtractedName}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Applicant}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Inventor}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Party}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Classification}) = StructTypes.Mutable()
StructTypes.StructType(::Type{ClassificationsIpcr}) = StructTypes.Mutable()
StructTypes.StructType(::Type{ClassificationsCpc}) = StructTypes.Mutable()
StructTypes.StructType(::Type{ClassificationsNational}) = StructTypes.Mutable()
StructTypes.StructType(::Type{DocumentId}) = StructTypes.Mutable()
StructTypes.StructType(::Type{PatentCitation}) = StructTypes.Mutable()
StructTypes.StructType(::Type{NPLCitation}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Citation}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Cited}) = StructTypes.Mutable()
StructTypes.StructType(::Type{CitingDoc}) = StructTypes.Mutable()
StructTypes.StructType(::Type{CitedBy}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Biblio}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Member}) = StructTypes.Mutable()
StructTypes.StructType(::Type{SimpleFamily}) = StructTypes.Mutable()
StructTypes.StructType(::Type{ExtendedFamily}) = StructTypes.Mutable()
StructTypes.StructType(::Type{LensFamily}) = StructTypes.Mutable()
StructTypes.StructType(::Type{LegalStatus}) = StructTypes.Mutable()
StructTypes.StructType(::Type{ClaimText}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Claim}) = StructTypes.Mutable()
StructTypes.StructType(::Type{Abstract}) = StructTypes.Mutable()
StructTypes.StructType(::Type{LensApplication}) = StructTypes.Mutable()

lensid(a::LensApplication) = a.lens_id
jurisdiction(a::LensApplication) = a.jurisdiction
docnr(a::LensApplication) = a.doc_number
kind(a::LensApplication) = a.kind
date(a::LensApplication) = a.date_published
status(a::LensApplication) = a.legal_status
type(a::LensApplication) = a.publication_type

title(a::LensApplication) = a.biblio.invention_title
abstract(a::LensApplication) = a.abstract

inventors(a::LensApplication) = [inv.extracted_name.value for inv in a.biblio.parties.inventors]
applicants(a::LensApplication) = [app.extracted_name.value for app in a.biblio.parties.applicants]

siblings(a::LensApplication) = a.families.simple_family.members
siblings_extended(a::LensApplication) = a.families.extended_family.members
family_size_simple(a::LensApplication) = a.families.simple_family.size
family_size_extended(a::LensApplication) = a.families.extended_family.size

classification(a::LensApplication) = a.biblio.classifications_cpc.classifications

_is_patcit(x) = isdefined(x, :patcit)
_is_nplcit(x) = isdefined(x, :nplcit)

function cites(a::LensApplication)
    cit = a.biblio.references_cited.citations
    idx = findall(_is_patcit, cit)
    [c.patcit for c in cit[idx]]
end

function cites_npl(a::LensApplication)
    cit = a.biblio.references_cited.citations
    idx = findall(_is_nplcit, cit)
    [c.nplcit for c in cit[idx]]
end

cites_count(a::LensApplication) = a.biblio.references_cited.patent_count
cites_count_npl(a::LensApplication) = a.biblio.references_cited.npl_count

citedby(a::LensApplication) = a.biblio.cited_by.patents
citedby_count(a::LensApplication) = a.biblio.cited_by.patent_count


text:
#using MetidaReports
using Test
using DataFrames, CSV, CategoricalArrays

path     = dirname(@__FILE__)
io       = IOBuffer();

@testset "  Test" begin
    ds  = CSV.File(path*"/csv/rds1.csv") |> DataFrame

    be = MetidaReports.bioquivalence(ds; variable = nothing, subject = :subject,
    period = :period, formulation = :treatment,
    sequence = :sequence, design = Symbol("2x2x4"))

end


text:
module cemgym

# this is a port of the example found at
#	http://rl-gym-doc.s3-website-us-west-2.amazonaws.com/mlss/lab1.html#starter

using Reinforce
using OpenAIGym
# using Distributions
using Transformations
using StochasticOptimization
using LearningStrategies

using MLPlots; gr(size=(500,500))

# ----------------------------------------------------------------

"""
Wraps a LearnBase.Transformation which converts an input vector to action values.
For discrete actions, it chooses the action which produces the highest value.
For continuous (interval) actions, it squashes actions to [0,1].
"""
struct TransformPolicy{T} <: AbstractPolicy
	trans::T
end

Transformations.params(tp::TransformPolicy) = params(tp.trans)

get_action(A::DiscreteSet, ŷ::AbstractVector) = A[indmax(ŷ)]
function get_action(A::IntervalSet, ŷ)
    val = 0.5 * (1.0 + clamp(first(ŷ), -1.0, 1.0))
    val * (A.hi - A.lo) + A.lo
end

# # discrete: our action is the action which maximizes the affine transform
# function Reinforce.action(π::TransformPolicy, r, s, A::DiscreteSet)
# 	A[indmax(transform!(π.trans, s))]
# end

# most actionsets pass through to get_action
Reinforce.action(π::TransformPolicy, r, s, A) = get_action(A)

# get an action where each part of the TupleSet is associated with part of the
# transformation output
function Reinforce.action(π::TransformPolicy, r, s, A::TupleSet)
    ŷ = transform!(π.trans, s)
    a = []
    i = 0
    for Aᵢ in A
        nᵢ = length(Aᵢ)
        push!(a, get_action(Aᵢ, view(ŷ, i+1:i+nᵢ)))
        i += nᵢ
    end
    # @show a
    a
end

# # continuous: return the transform value, squashed to [0,1]
# # TODO: remove the squashing when a "Affine + Sigmoid" transformation is available
# function Reinforce.action(π::TransformPolicy, r, s, A::IntervalSet)
# 	Transformations.sigmoid(transform(π.trans, s)[1]) * (A.amax-A.amin) + A.amin
# end

# ----------------------------------------------------------------

# function myplot(t, hist_min, hist_mean, hist_max, anim=nothing)
# 	(env,i,sars) -> if mod1(t,3)==1 && mod1(i,10)==1
# 		plot(
# 			plot(hist_mean, c=:black, fill=((hist_min,hist_max), 0.2), title="Progress", leg=false),
# 			plot(env, title = "Episode: $t  Iter: $i")
# 		)
# 		if anim == nothing
# 			gui()
# 		else
# 			frame(anim)
# 		end
# 	else
# 		return
# 	end
# end

# ----------------------------------------------------------------

# initialize a policy, do the learning, then return the policy
function do_cem_test(sublearners...; env = GymEnv("SoccerEmptyGoal-v0"),
                       maxsteps = 200, # max number of steps in one episode
                       maxiter = 1000, # max learning iterations
					   # noise_max = 1.0,
					   # noise_steps = 20,
					   # noise_func = t -> max(noise_max - t/noise_steps, 0.0),
					   kw...)
    s = state(env)
    A = actions(env,s)

    nin = length(s)
    nout = length(A)
    @show nin, nout

    t = nnet(nin, nout, [2], :softplus)
    @show t

    policy = TransformPolicy(t)
    cem = CrossEntropyMethod(;maxsteps=maxsteps, kw...)

    tp = TracePlot(2, layout=@layout([a;b{0.2h}]))
    tracer = IterFunction((policy,i) -> begin
        mod1(i,10)==1 || return

        #run one episode
        R,N = episode!(env, policy, stepfunc = OpenAIGym.render)
        @show R, N
        push!(tp, i, [R,N])
        gui(tp.plt)
    end)

    learner = strategy(cem, tracer, sublearners...; maxiter=maxiter, kw...)
    learn!(policy, learner, repeated(env))

    @show policy cem
end

end #module


text:
"""
    Parameterizations

Useful parameterizations that map various constrained parameters to
the real line or ``R^n``.  The naming convention follows:

```julia
p = XXX_param(x, args...)
x = XXX_value(p, args...)
log_jacobian = XXX_logjac(x, p, args...)
```

where `log_jacobian` gives ``\\log \\left| \\frac{\\mathrm{d}
x}{\\mathrm{d} p} \\right|``, which is the approprate factor to multiply
a density in `x` to produce a density in `p`.  In other words,
`log_jacobian` is the factor that should be introduced into the prior
to properly account for the reparameterization.

Scalar parameterizations also have vector-based versions.
"""
module Parameterizations

"""
    bounded_param(x, low, high)

Variables that are bounded between `low` and `high`.
"""
bounded_param(x, low, high) = log(x-low) - log(high-x)

function bounded_value(p, low, high)
    if p > 0
        ep = exp(-p)
        (high + ep*low)/(one(p) + ep)
    else
        ep = exp(p)
        (high*ep + low)/(one(p) + ep)
    end
end

bounded_value(x::Vector, low, high) = [bounded_value(x[i], low, high) for i in eachindex(x)]

bounded_logjac(x, p, low, high) = log((high-x)*(x-low)/(high-low))
function bounded_logjac(x::Vector, p::Vector, low, high)
    lj = 0.0
    for i in eachindex(x)
        lj += bounded_logjac(x[i], p[i], low, high)
    end
    lj
end

"""
    positive_param(x)

Variable that is constrained to be positive.
"""
function positive_param(x)
    log(x)
end
positive_param(x::Vector) = [positive_param(x[i]) for i in eachindex(x)]

function positive_value(p)
    exp(p)
end
positive_value(p::Vector) = [positive_value(p[i]) for i in eachindex(p)]

function positive_logjac(x, p)
    p
end
function positive_logjac(x::Vector, p::Vector)
    lj = zero(x[1])
    for i in eachindex(x)
        lj += positive_logjac(x[i], p[i])
    end
    lj
end

"""
    increasing_params(x)

Variables that are constrained to be increasing: `x[1] < x[2] < ...`.
"""
function increasing_params(x)
    n = length(x)
    for i in 2:n
        @assert(x[i] > x[i-1], "values must be increasing")
    end

    p = zeros(x)
    p[1] = x[1]
    for i in 2:n
        p[i] = log(x[i]-x[i-1])
    end

    p
end

function increasing_values(p)
    n = length(p)

    x = zeros(p)
    x[1] = p[1]
    for i in 2:n
        x[i] = x[i-1] + exp(p[i])
    end

    x
end

function increasing_logjac(x, p)
    n = length(p)

    lj = 0.0
    for i in 2:n
        lj += p[i]
    end

    lj
end

"""
    increasing_positive_params(x)

Variables that are both positive and increasing.

"""
function increasing_positive_params(x)
    @assert all(x[2:end] .> x[1:end-1])

    n = size(x, 1)
    p = zeros(n)

    p[1] = log(x[1])
    for i in 2:n
        p[i] = log(x[i]-x[i-1])
    end

    p
end

function increasing_positive_values(p)
    n = size(p, 1)
    x = zeros(n)

    x[1] = exp(p[1])
    for i in 2:n
        x[i] = x[i-1] + exp(p[i])
    end

    x
end

function increasing_positive_logjac(x, p)
    sum(p)
end

"""
    simplex_params(x)

Variables that are constrained to be positive and sum to 1.  The
resulting parameterization is reduced in dimension by 1.
"""
function simplex_params(x)
    @assert(abs(sum(x) - 1) < 1e-8, "values must sum to one")

    n = length(x)

    p = zeros(n-1)
    remaining = one(p[1])
    for i in 1:n-1
        p[i] = bounded_param(x[i], zero(x[i]), remaining)
        remaining = remaining - x[i]
    end

    p
end

function simplex_values(p)
    n = length(p)

    x = zeros(n+1)
    remaining = one(p[1])
    for i in 1:n
        x[i] = bounded_value(p[i], zero(p[i]), remaining)
        remaining = remaining - x[i]
    end

    x[n+1] = remaining

    x
end

function simplex_logjac(x, p)
    lj = 0.0
    remaining = 1.0

    for i in 1:length(p)
        lj += bounded_logjac(x[i], p[i], zero(x[i]), remaining)
        remaining = remaining - x[i]
    end

    lj
end

"""
    unit_disk_param(z)

Two-dimensional variable `z` that is constrained to live in the unit
disk: ``|z| < 1``.
"""
function unit_disk_param(z)
    x = z[1]
    y = z[2]

    r2 = x*x + y*y

    a = 1.0 - r2

    Float64[x/a, y/a]
end

function unit_disk_value(p)
    u = p[1]
    v = p[2]

    R2 = u*u + v*v
    R = sqrt(R2)

    r = 0.5*(sqrt(4.0*R2 + 1.0) - 1.0)/R

    x = u*r/R
    y = v*r/R

    Float64[x, y]
end

function unit_disk_logjac(z, p)
    x = z[1]
    y = z[2]

    r2 = x*x + y*y

    3.0*log1p(-r2) - log1p(r2)
end

end


text:
using ExactConvolution
using Test

function naive_array2bigint(arr, num_bits)
	@assert all(>=(0), arr)
	ans = BigInt(0)
	for i ∈ reverse(arr)
		ans <<= num_bits
		ans |= i
	end
	ans
end

function random_input(;len=300, bits=24, max_bits=100)
	l = rand(1:len)
	arr = rand(0:2^bits-1, l)
	num_bits = rand(bits+1:max_bits)
	while num_bits%5!=0
		num_bits+=1
	end
	return arr,num_bits
end

@testset "Array->BigInt" begin
    function test_input(arr, num_bits)
		r1 = ExactConvolution.array2bigint(arr, num_bits)
		r2 = naive_array2bigint(arr, num_bits)
		@test r1==r2
	end
    for _ ∈ 1:100000
		test_input(random_input(len=3, bits=3, max_bits=4)...)
	end

	for _ ∈ 1:100000
		test_input(random_input(len=5, bits=10, max_bits=20)...)
	end

	for _ ∈ 1:1000
		test_input(random_input()...)
	end
end

@testset "BigInt->Array" begin
    function run_bigint2array(arr, num_bits)
        bi = naive_array2bigint(arr, num_bits)
        bigint2array(BigInt, bi, length(arr), num_bits)
    end

    function test_input(arr, num_bits)
		arr2 = run_bigint2array(arr, num_bits)
		@test arr == arr2
	end

	for _ ∈ 1:1000
		test_input(random_input(len=5, bits=3, max_bits=6)...)
	end

	for _ ∈ 1:1000
		test_input(random_input()...)
	end
end

function naive_conv(a, b)
    n = length(a) + length(b) - 1
    arr = zeros(BigInt, n)
    for (i, v1) ∈ enumerate(a), (j, v2) ∈ enumerate(b)
        arr[i+j-1] += big(v1) * v2
    end
    arr
end

function random_inputs_conv(;len=300, bits=20, has_neg=true)
	rand_array() = begin
		l = rand(1:len)
		arr = rand(0:(2^bits-1), l)
		rand(false:has_neg) && (arr .-= (2^(bits-1)))
		arr
	end
	rand_array(),rand_array()
end

@testset "Exact Convolutions" begin
    function test_input(arr1, arr2)
		r1 = exact_conv(arr1,arr2)
		r2 = naive_conv(arr1,arr2)
		@test r1==r2
	end

	for _ ∈ 1:100000
		test_input(random_inputs_conv(len=5, bits=5, has_neg=false)...)
	end

	for _ ∈ 1:100000
		test_input(random_inputs_conv(len=5, bits=5)...)
	end

	for _ ∈ 1:1000
		test_input(random_inputs_conv()...)
	end
end


text:
"
solves AL * C = C * A in-place
"
function uniform_leftorth!(AL,CR, A; tol = Defaults.tolgauge, maxiter = Defaults.maxiter)
    #(_,CR[end]) = leftorth!(CR[end], alg = TensorKit.QRpos());
    normalize!(CR[end]);

    iteration = 1;
    delta = 2*tol;

    while iteration < maxiter && delta > tol

        if iteration>10 #when qr starts to fail, start using eigs - should be kw arg

            alg = Arnoldi(krylovdim = 30, tol = max(delta*delta,tol/10),maxiter=maxiter)

            (vals,vecs) = @closure eigsolve(CR[end], 1, :LM,alg) do x
                transfer_left(x,A,AL)
            end

            (_,CR[end]) = leftorth!(vecs[1],alg=TensorKit.QRpos())
        end

        cold = CR[end]

        for loc in 1:length(AL)
            AL[loc] = _transpose_front(CR[mod1(loc-1,end)]*_transpose_tail(A[loc]))
            AL[loc], CR[loc] = leftorth!(AL[loc], alg = QRpos())
            normalize!(CR[loc])
        end

        #update delta
        if domain(cold) == domain(CR[end]) && codomain(cold) == codomain(CR[end])
            delta = norm(cold-CR[end])
        end

        iteration += 1
    end

    delta>tol && @warn "leftorth failed to converge $(delta)"

    AL,CR
end

"
solves C * AR = C * A in-place
"
function uniform_rightorth!(AR,CR,A; tol = Defaults.tolgauge, maxiter = Defaults.maxiter)
    #(CR[end],_) = rightorth!(CR[end], alg = TensorKit.LQpos());
    normalize!(CR[end])

    iteration = 1;
    delta = 2*tol;
    while iteration<maxiter && delta>tol

        if iteration>10#when qr starts to fail, start using eigs
            alg = Arnoldi(krylovdim = 30, tol = max(delta*delta,tol/10),maxiter=maxiter)
            #Projection of the current guess onto its largest self consistent eigenvector + isolation of the unitary part

            (vals,vecs) = @closure eigsolve(CR[end], 1, :LM,alg) do x
                transfer_right(x,A,AR)
            end
            (CR[end],_) = rightorth!(vecs[1],alg=TensorKit.LQpos())
        end

        cold = CR[end]
        for loc in length(AR):-1:1
            AR[loc] = A[loc]*CR[loc]

            CR[mod1(loc-1,end)], temp = rightorth!(_transpose_tail(AR[loc]), alg=LQpos())
            AR[loc] = _transpose_front(temp)
            normalize!(CR[mod1(loc-1,end)])
        end

        #update counters and delta
        if domain(cold) == domain(CR[end]) && codomain(cold) == codomain(CR[end])
            delta = norm(cold-CR[end])
        end

        iteration += 1
  end

  delta>tol && @warn "rightorth failed to converge $(delta)"

  return AR, CR
end


text:
"""
    _build_univariate_milp_relaxation!(m,x,y,function_data,pre_base_name)

Return a MILPRelaxation object with constraint and RHS information of the MILP
formulation of the polyhedral relaxation.
"""
function _build_univariate_milp_relaxation!(
    m::JuMP.Model,
    x::JuMP.VariableRef,
    y::JuMP.VariableRef,
    univariate_function_data::UnivariateFunctionData,
    pre_base_name::AbstractString,
)::FormulationInfo
    sec_vs, tan_vs = _collect_vertices(univariate_function_data)
    formulation_info = FormulationInfo()

    # create variables
    num_vars = length(univariate_function_data.partition) - 1

    delta_1 =
        formulation_info.variables[:delta_1] = @variable(
            m,
            [1:num_vars],
            lower_bound = 0.0,
            upper_bound = 1.0,
            base_name = pre_base_name * "delta_1"
        )
    delta_2 =
        formulation_info.variables[:delta_2] = @variable(
            m,
            [1:num_vars],
            lower_bound = 0.0,
            upper_bound = 1.0,
            base_name = pre_base_name * "delta_2"
        )
    z =
        formulation_info.variables[:z] = @variable(
            m,
            [1:num_vars],
            binary = true,
            base_name = pre_base_name * "z"
        )

    # add x constraints
    formulation_info.constraints[:x] = @constraint(
        m,
        x ==
        sec_vs[1][1] + sum(
            delta_1[i] * (tan_vs[i][1] - sec_vs[i][1]) +
            delta_2[i] * (sec_vs[i+1][1] - sec_vs[i][1]) for i in 1:num_vars
        )
    )

    # add y constraints
    formulation_info.constraints[:y] = @constraint(
        m,
        y ==
        sec_vs[1][2] + sum(
            delta_1[i] * (tan_vs[i][2] - sec_vs[i][2]) +
            delta_2[i] * (sec_vs[i+1][2] - sec_vs[i][2]) for i in 1:num_vars
        )
    )

    # add first delta constraint
    formulation_info.constraints[:first_delta] =
        @constraint(m, delta_1[1] + delta_2[1] <= 1)

    # add linking constraints between delta_1, delta_2 and z
    formulation_info.constraints[:below_z] =
        @constraint(m, [i = 2:num_vars], delta_1[i] + delta_2[i] <= z[i-1])
    formulation_info.constraints[:above_z] =
        @constraint(m, [i = 2:num_vars], z[i-1] <= delta_2[i-1])

    return formulation_info
end


text:
using LinearAlgebra,Statistics,MultivariateStats,Roots,HypothesisTests,RCall,Dierckx,Plots
import Base: length, push!, deleteat!

"ePPR Log Options"
Base.@kwdef mutable struct ePPRLog
    debug = false
    plot = false
    io = nothing
    dir = nothing
end
function (log::ePPRLog)(msg;file="ePPRLog.txt",close=false)
    if isnothing(log.io)
        if isnothing(log.dir)
            io=stdout
        else
            isdir(log.dir) || mkpath(log.dir)
            io = open(joinpath(log.dir,file),"a")
            log.io = io
        end
    else
        io=log.io
    end
    println(io,msg)
    flush(io)
    close && !isnothing(log.io) && Base.close(log.io)
end
function (log::ePPRLog)(msg::Plots.Plot;file="Model.png")
    if isnothing(log.dir)
        display(msg)
    else
        isdir(log.dir) || mkpath(log.dir)
        savefig(msg,joinpath(log.dir,file))
    end
end

"""
extended Projection Pursuit Regression Model

```math
\\hat{y}_i=\\bar{y}+\\sum_{d=0}^D\\sum_{m=1}^{M_d}\\beta_{m,d}\\phi_{m,d}(\\alpha_{m,d}^Tx_{i-d}), where

\\frac{1}{n}\\sum_{i=1}^n\\phi_{m,d}(\\alpha_{m,d}^Tx_{i-d})=0, \\frac{1}{n}\\sum_{i=1}^n\\phi_{m,d}^2(\\alpha_{m,d}^Tx_{i-d})=1
```
"""
Base.@kwdef mutable struct ePPRModel
    "𝑦̄"
    ymean::Float64 = 0.0
    "vector of β for each term"
    beta::Vector{Float64} = []
    "vector of Φ for each term"
    phi::Vector = []
    "vector of α for each term"
    alpha::Vector{Vector{Float64}} = []
    "vector of (temporal, spatial) index for each term"
    index::Vector = []
    "vector of ``\\phi_{m,d}(\\alpha_{m,d}^TX_{-d})`` for each term"
    phivalues::Vector{Vector{Float64}} = []
    "vector of trustregionsize for each term"
    trustregionsize::Vector{Float64} = []
    "γ"
    residuals::Vector{Float64} = []
end
"Number of terms in the model"
length(m::ePPRModel) = length(m.beta)
function deleteat!(model::ePPRModel,i::Int)
    deleteat!(model.beta,i)
    deleteat!(model.phi,i)
    deleteat!(model.alpha,i)
    deleteat!(model.index,i)
    deleteat!(model.phivalues,i)
    deleteat!(model.trustregionsize,i)
    return model
end
function push!(model::ePPRModel,β::Float64,Φ,α::Vector{Float64},index,Φvs::Vector{Float64},trustregionsize::Float64=1.0)
    push!(model.beta,β)
    push!(model.phi,Φ)
    push!(model.alpha,α)
    push!(model.index,index)
    push!(model.phivalues,Φvs)
    push!(model.trustregionsize,trustregionsize)
    return model
end
function getterm(model::ePPRModel,i::Int)
    return model.beta[i],model.phi[i],model.alpha[i],model.index[i],model.phivalues[i],model.trustregionsize[i]
end
function setterm!(model::ePPRModel,i::Int,β::Float64,Φ,α::Vector{Float64},index,Φvs::Vector{Float64},trustregionsize::Float64=1.0)
    model.beta[i]=β
    model.phi[i]=Φ
    model.alpha[i]=α
    model.index[i]=index
    model.phivalues[i]=Φvs
    model.trustregionsize[i]=trustregionsize
    return model
end
function clean!(model::ePPRModel)
    model.phivalues=[]
    model.trustregionsize=[]
    model.residuals=[]
    return model
end

"Cross Validation Parameters for ePPR"
Base.@kwdef mutable struct ePPRCrossValidation
    "Percentage of sample for training"
    trainpercent::Float64 = 0.9
    "N-Fold of training sample"
    trainfold::Int = 5
    "N-Fold of one `trainfold` sample for testing training"
    traintestfold::Int = 15
    "each combination of train and traintest fold"
    trainsets = []
    "N-Fold of testing sample for testing model"
    testfold::Int = 15
    "`testfold` of testing sample"
    tests = []
    "which trainset for training"
    trainsetindex::Int = 1
    "significent level to accept null hypothesis that current model is no better than reference model"
    h0level::Float64 = 0.01
    "significent level to accept alternative hypothesis that current model is better than reference model"
    h1level::Float64 = 0.01
    "Correlation between response and model prediction on testsets of current trainset"
    modeltraintestcor = []
    "Correlation between response and model prediction on testsets"
    modeltestcor = []
    "Correlation between response and model prediction"
    modelcor = []
end

"Hyper Parameters for ePPR"
Base.@kwdef mutable struct ePPRHyperParams
    """memory size to pool for nonlinear time interaction, ndelay=1 for linear time interaction.
    only first delay terms in `nft` is used for nonlinear time interaction."""
    ndelay::Int = 1
    "number of forward terms for each delay. [3, 2, 1] means 3 spatial terms for delay 0, 2 for delay 1, 1 for delay 2"
    nft::Vector{Int} = [3,3,3]
    "penalization parameter λ"
    lambda::Float64 = 8
    "Φ Spline degree of freedom"
    phidf::Int = 5
    "minimum number of backward terms"
    mnbt::Int = 1
    "whether to fit all spatial terms before moving to next temporal delay"
    spatialtermfirst::Bool = true
    "α priori for penalization"
    alphapenaltyoperator = []
    "`(lossₒ-lossₙ)/lossₒ`, forward converge rate threshold to decide a saturated iteration"
    forwardconvergerate::Float64 = 0.01
    "`(lossₒ-lossₙ)/lossₒ`, refit converge rate threshold to decide a saturated iteration"
    refitconvergerate::Float64 = 0.001
    "number of consecutive saturated iterations to decide a new term"
    nsaturatediteration::Int = 2
    "maximum number of iterations to fit a new term"
    newtermmaxiteration::Int = 50
    "initial size of trust region"
    trustregioninitsize::Float64 = 1
    "maximum size of trust region"
    trustregionmaxsize::Float64 = 100
    "η of trust region"
    trustregioneta::Float64 = 0.2
    "maximum iterations of trust region"
    trustregionmaxiteration::Int = 50
    "dimension of image"
    imagesize = ()
    "row vector of blank image"
    blankimage = []
    "drop term index between backward models"
    droptermindex = []
    "ePPR Cross Validation"
    cv::ePPRCrossValidation = ePPRCrossValidation()
    "Valid Image Region"
    xindex::Vector{Int} = Int[]
    "Index iⱼ where image sequence breaks between x[iⱼ-1,:] and x[iⱼ,:]"
    xbreak::Vector{Int} = Int[]
    "maximum iterations of hyperparameter search"
    hypermaxiteration::Int = 15
    "number of consecutive no model iterations to stop hyperparameter search"
    nhypernomodeliteration::Int = 3
    "number of consecutive chance iterations to stop hyperparameter search"
    nhyperchanceiteration::Int = 3
    "number of consecutive saturated iterations to stop hyperparameter search"
    nhypersaturatediteration::Int = 2
    "scale factor for λ search"
    lambdascale::Float64 = 0.5
end
function ePPRHyperParams(nrow::Int,ncol::Int;xindex::Vector{Int}=Int[],ndelay::Int=1,nft::Vector{Int}=[3,3,3],lambda::Real=32,lambdascale::Real=0.5,blankcolor::Real=127)
    hp = ePPRHyperParams(;imagesize=(nrow,ncol),xindex,ndelay,nft,lambda,lambdascale)
    hp.blankimage = fill(blankcolor,1,nrow*ncol)
    hp.alphapenaltyoperator = laplacian2dmatrix(nrow,ncol)
    return hp
end

"Image sequences at delay `d`"
function delayx(x,d,hp,xi=[])
    if d<=0
        tx = x
    else
        tx = [repeat(hp.blankimage,outer=(d,1));x[1:end-d,:]]
        if !isempty(hp.xbreak)
            for bi in hp.xbreak
                tx[bi-d:bi-1,:]=repeat(hp.blankimage,outer=(d,1))
            end
        end
    end
    isempty(xi) ? tx : tx[xi,:]
end

"Pool images in delay window"
function delaywindowpool!(x::Matrix,hp::ePPRHyperParams,log::ePPRLog)
    if isempty(hp.xindex)
        vx = x
    else
        vx = x[:,hp.xindex]
        hp.blankimage = hp.blankimage[:,hp.xindex]
        hp.alphapenaltyoperator=hp.alphapenaltyoperator[hp.xindex,hp.xindex]
    end
    hp.ndelay<=1 && return vx

    log.debug && log("Nonlinear Time Interaction, pool x[i-$(hp.ndelay-1):i, :] together ...")
    dwpx=vx
    for d in 1:hp.ndelay-1
        dwpx = [dwpx delayx(vx,d,hp)]
    end
    hp.alphapenaltyoperator = delaywindowpooloperator(hp.alphapenaltyoperator,hp.ndelay)
    return dwpx
end

function delaywindowpooloperator(operator::Matrix,ndelay::Int=1)
    ndelay<=1 && return operator

    nr,nc=size(operator)
    dwpo = zeros(ndelay*nr,ndelay*nc)
    for d in 0:ndelay-1
        dwpo[(1:nr).+d*nr, (1:nc).+d*nc] = operator
    end
    return dwpo
end

"Model prediction on training data"
(m::ePPRModel)() = m.ymean .+ mapreduce((b,ps)->b*ps, (i,j)->i.+j, m.beta, m.phivalues)
"""
Model prediction on data `x`

1. x: matrix with one image per row
2. hp: hyper parameters
3. xi: image indices on which predictions are made
"""
function (m::ePPRModel)(x::Matrix,hp::ePPRHyperParams,xi=[])
    ti=map(i->i.t,m.index);ut=unique(ti);uti=[findall(ti.==t) for t in ut]
    ŷ = m.ymean
    for i in eachindex(ut)
        tx = delayx(x,ut[i],hp,xi)
        for j in uti[i]
            ŷ = ŷ .+ m.beta[j] * m.phi[j](tx*m.alpha[j])
        end
    end
    return ŷ
end

"""
Data Partition for cross validation

1. cv: cross validation
2. n: number of samples
3. log: log options
"""
function cvpartition!(cv::ePPRCrossValidation,n::Int,log::ePPRLog)
    ntrain = cv.trainpercent*n
    ntrainfold = ntrain/cv.trainfold
    ntraintestfold = Int(floor(ntrainfold/cv.traintestfold))
    ntrainfold = ntraintestfold*cv.traintestfold
    ntrain = ntrainfold*cv.trainfold
    trainsets=[]
    for tf in 0:cv.trainfold-1
        traintest = [tf*ntrainfold .+ (1:ntraintestfold) .+ ttf*ntraintestfold for ttf in 0:cv.traintestfold-1]
        train = setdiff(1:ntrain,tf*ntrainfold .+ (1:ntrainfold))
        push!(trainsets,(train=train,traintest=traintest))
    end
    cv.trainsetindex=cv.trainfold
    ntestfold = Int(floor((n-ntrain)/cv.testfold))
    tests = [ntrain .+ (1:ntestfold) .+ tf*ntestfold for tf in 0:cv.testfold-1]
    log.debug && log("Cross Validation Data Partition: n = $n, ntrain = $ntrain in $(cv.trainfold)-fold, ntrainfold = $ntrainfold in $(cv.traintestfold)-fold, ntest = $(ntestfold*cv.testfold) in $(cv.testfold)-fold")
    cv.trainsets=trainsets;cv.tests=tests
    return cv
end

"Choose best model among models by cross validation"
function cvmodel(models::Vector{ePPRModel},x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog)
    log.debug && log("ePPR Models Cross Validation ...")
    train = hp.cv.trainsets[hp.cv.trainsetindex].train;traintest = hp.cv.trainsets[hp.cv.trainsetindex].traintest;test=hp.cv.tests
    # response and model predication
    traintestpredications = map(m->map(i->m(x,hp,i),traintest),models)
    traintestys = map(i->y[i],traintest)
    # correlation between response and predication
    traintestcors = map(mps->cor.(traintestys,mps),traintestpredications)
    hp.cv.modeltraintestcor=[];hp.cv.modeltestcor=[] # clear model test correlation
    log.plot && log(plotcor(models,traintestcors),file="Models_Goodness (λ=$(hp.lambda)).png")
    # find the model no worse than models with more terms and better than models with less terms
    mi=0;nmodel=length(models)
    for rm in 1:nmodel
        moretermp = [pvalue(SignedRankTest(traintestcors[rm],traintestcors[m]),tail=:left) for m in rm+1:nmodel]
        if rm==1 && (1==nmodel || all(moretermp .> hp.cv.h0level))
            mi=rm;break
        end
        if rm==nmodel || all(moretermp .> hp.cv.h0level)
            lesstermp = [pvalue(SignedRankTest(traintestcors[m],traintestcors[rm]),tail=:left) for m in 1:rm-1]
            if all(lesstermp .< hp.cv.h1level)
                mi=rm;break
            end
        end
    end
    if mi==0
        log.debug && log("No model not worse than models with more terms and better than models with less terms.")
        return nothing
    end
    model = deepcopy(models[mi])
    log.debug && log("$(mi)th model with $(length(model)) terms is chosen.")

    # find drop terms that do not improve model predication
    droptermp = [pvalue(SignedRankTest(traintestcors[m-1],traintestcors[m]),tail=:left) for m in 2:nmodel]
    notimprove = findall(droptermp .> hp.cv.h0level)
    # find drop term models with chance level(zero correlation) predication
    modelp = [pvalue(SignedRankTest(traintestcors[m]),tail=:both) for m in 2:nmodel]
    notpredictive = findall(modelp .> hp.cv.h0level)

    poorterm = hp.droptermindex[union(notimprove,notpredictive)]
    # spurious terms in the selected model
    spuriousterm = findall(in(poorterm),model.index)
    if !isempty(spuriousterm)
        log.debug && log("Model drop spurious term: $(model.index[spuriousterm]).")
        foreach(i->deleteat!(model,i),sort(spuriousterm,rev=true))
    end
    length(model)==0 && return nothing
    model = eppr(model,x[train,:],y[train],hp,log)
    hp.cv.modeltraintestcor = map(i->cor(y[i],model(x,hp,i)),traintest)
    hp.cv.modeltestcor = map(i->cor(y[i],model(x,hp,i)),test)
    return model
end

"extended Projection Pursuit Regression with cross validation"
function epprcv(x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog=ePPRLog())
    n = length(y);n !=size(x,1) && error("Length of x and y does not match!")
    cvpartition!(hp.cv,n,log)
    px = delaywindowpool!(x,hp,log)
    log.debug && log("Choose $(hp.cv.trainsetindex)th trainset.")
    train = hp.cv.trainsets[hp.cv.trainsetindex].train
    models = eppr(px[train,:],y[train],hp,log)
    model = cvmodel(models,px,y,hp,log)
    log.debug && log("Cross Validated ePPR Done.",close=true)
    return model,models
end

"extended Projection Pursuit Regression with cross validation and hyper parameter search"
function epprhypercv(x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog=ePPRLog())
    n = length(y);n !=size(x,1) && error("Length of x and y does not match!")
    cvpartition!(hp.cv,n,log)
    px = delaywindowpool!(x,hp,log)
    log.debug && log("Choose $(hp.cv.trainsetindex)th trainset.")
    train = hp.cv.trainsets[hp.cv.trainsetindex].train

    hi=0;hypermodel=[];hypermodels=[];λs=[];modelcor=[];saturatediteration=0;chanceiteration=0;nomodeliteration=0
    for i in 1:hp.hypermaxiteration
        log.debug && log("HyperParameter Search: λ = $(hp.lambda) ...")
        models = eppr(px[train,:],y[train],hp,log)
        model = cvmodel(models,px,y,hp,log)
        if isnothing(model)
            nomodeliteration+=1
            if nomodeliteration>=hp.nhypernomodeliteration
                hi=-2;break
            else
                hp.lambda *= hp.lambdascale
            end
        else
            nomodeliteration=0
            push!(hypermodel,model);push!(hypermodels,models);push!(λs,hp.lambda);push!(modelcor,[hp.cv.modeltraintestcor;hp.cv.modeltestcor])
            chancep = pvalue(SignedRankTest(modelcor[end]),tail=:both)
            if chancep > hp.cv.h0level
                chanceiteration+=1
                if chanceiteration>=hp.nhyperchanceiteration
                    hi=-1;break
                else
                    hp.lambda *= hp.lambdascale;continue
                end
            else
                chanceiteration=0
            end

            if length(modelcor)==1
                hp.lambda *=hp.lambdascale
            else
                improvep = pvalue(SignedRankTest(modelcor[end],modelcor[end-1]),tail=:right)
                if improvep < hp.cv.h1level
                    saturatediteration=0
                    hp.lambda *= hp.lambdascale
                else
                    impairp = pvalue(SignedRankTest(modelcor[end],modelcor[end-1]),tail=:left)
                    if impairp < hp.cv.h1level
                        hi=length(modelcor)-1;break
                    else
                        saturatediteration+=1
                        if saturatediteration>=hp.nhypersaturatediteration
                            hi = hp.lambdascale > 1 ? length(modelcor) : length(modelcor)-hp.nhypersaturatediteration;break
                        else
                            hp.lambda *= hp.lambdascale
                        end
                    end
                end
            end
        end
    end
    log.plot && !isempty(modelcor) && log(plotcor(λs,modelcor,xlabel="λ"),file="λ_Model_Goodness.png")
    if hi==-2
        log.debug && log("No valid λ and model.",close=true)
        return nothing,[]
    elseif hi==-1
        log.debug && log("No predictive λ and model.",close=true)
        return nothing,[]
    elseif hi==0
        if isempty(modelcor)
            log.debug && log("No valid λ and model.",close=true)
            return nothing,[]
        else
            hi=argmax(median.(modelcor))
        end
    end
    log.debug && log("HyperParameter search done with best λ = $(λs[hi]).",close=true)
    hp.lambda = λs[hi];hp.cv.modelcor = modelcor[hi]
    return hypermodel[hi],hypermodels[hi]
end

"""
extended Projection Pursuit Regression

Minimizing ``f=\\sum_{i=1}^N(y_i-\\hat{y}(x_i))^2+\\lambda\\sum_{d=0}^D\\sum_{m=1}^{M_d}\\Vert{L\\alpha_{m,d}}\\Vert^2``

1. x: matrix with one image per row
2. y: vector of response for each image
3. hp: hyper parameters
4. log: log options
"""
function eppr(x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog=ePPRLog())
    model = forwardstepwise(x,y,hp,log)
    model = refitmodel!(model,x,y,hp,log)
    models = backwardstepwise(model,x,y,hp,log)
end
function eppr(model::ePPRModel,x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog)
    refitmodelbeta!(forwardstepwise(model,x,y,hp,log),y,log)
end

"""
Gradually add term into the model based on terms of an existing model
"""
function forwardstepwise(m::ePPRModel,x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog)
    ti = map(i->i.t,m.index);si = map(i->i.s,m.index)
    ut = sort(unique(ti));uts = Dict(t=>si[ti.==t] for t in ut);utsn=Dict(t=>length(uts[t]) for t in ut)
    log.debug && log("ePPR Model Forward Stepwise ...")
    ym = mean(y);model = ePPRModel(ymean=ym,residuals=y.-ym)
    if hp.spatialtermfirst
        for t in ut
            tx = delayx(x,t,hp)
            for s in 1:utsn[t]
                log.debug && log("Fit Model (Temporal-$t, Spatial-$s) New Term ...")
                α = normalize(m.alpha[findfirst(==((t=t,s=uts[t][s])),m.index)], 2)
                β,Φ,α,Φvs = fitnewterm(tx,model.residuals,α,hp.phidf)
                model.residuals .-= β*Φvs
                push!(model,β,Φ,α,(t=t,s=s),Φvs)
            end
        end
    else
        for s in 1:maximum(values(utsn)), t in ut
            s>utsn[t] && continue
            log.debug && log("Fit Model (Temporal-$t, Spatial-$s) New Term ...")
            tx = delayx(x,t,hp)
            α = normalize(m.alpha[findfirst(==((t=t,s=uts[t][s])),m.index)], 2)
            β,Φ,α,Φvs = fitnewterm(tx,model.residuals,α,hp.phidf)
            model.residuals .-= β*Φvs
            push!(model,β,Φ,α,(t=t,s=s),Φvs)
        end
    end
    return model
end

"""
Gradually add term into the model
"""
function forwardstepwise(x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog)
    log.debug && log("ePPR Forward Stepwise ...")
    if hp.ndelay>1
        hp.nft=hp.nft[1:1]
    end
    ym = mean(y);model = ePPRModel(ymean=ym,residuals=y.-ym)
    if hp.spatialtermfirst
        for t in 0:length(hp.nft)-1
            tx = delayx(x,t,hp)
            for s in 1:hp.nft[t+1]
                log.debug && log("Fit (Temporal-$t, Spatial-$s) New Term ...")
                α = getinitialalpha(tx,model.residuals,log)
                β,Φ,α,Φvs,trustregionsize = fitnewterm(tx,model.residuals,α,hp,log)
                model.residuals .-= β*Φvs
                push!(model,β,Φ,α,(t=t,s=s),Φvs,trustregionsize)
            end
        end
    else
        for s in 1:maximum(hp.nft), t in 0:length(hp.nft)-1
            s>hp.nft[t+1] && continue
            log.debug && log("Fit (Temporal-$t, Spatial-$s) New Term ...")
            tx = delayx(x,t,hp)
            α = getinitialalpha(tx,model.residuals,log)
            β,Φ,α,Φvs,trustregionsize = fitnewterm(tx,model.residuals,α,hp,log)
            model.residuals .-= β*Φvs
            push!(model,β,Φ,α,(t=t,s=s),Φvs,trustregionsize)
        end
    end
    return model
end

"""
Estimate ridge parameter

1. X: design matrix without intercept
2. y: response vector

return:

- LW: L-W estimate of the ridge parameter for unscaled `X`
"""
function ridgeparam(X,y)
    n,p = size(X)
    β = llsq(X,y,bias=false)
    ŷ = X * β
    σ² = sum((y .- ŷ).^2) / (n - p)

    LW = (p - 2) * σ² * n/sum(ŷ.^2)
    Xscale² = vec(sum((X.^2)./n,dims=1))
    LW *= Xscale²
    return (;LW)
end

"""
A simplified version of `lm.ridge` in R package `MASS`.

1. X: design matrix without intercept
2. y: response vector
3. λ: ridge parameter for internal scaled `X`

return:

- coef: coefficients for unscaled `X`
- HKB: HKB estimate of the ridge parameter for internal scaled `X`
- LW: L-W estimate of the ridge parameter for internal scaled `X`
"""
function ridge_R(X,y,λ=0.0)
    n,p = size(X)
    Xscale = sqrt.(sum((X.^2)./n,dims=1))
    _X = X./Xscale
    _Xsvd = svd(_X)
    rhs = _Xsvd.U'*y
    lscoef = _Xsvd.V * (rhs./_Xsvd.S)
    lsfit = _X*lscoef
    σ² = sum((y.-lsfit).^2)/(n - p)
    HKB = (p - 2) * σ²/sum(lscoef.^2)
    LW = (p - 2) * σ² * n/sum(lsfit.^2)
    div = _Xsvd.S.^2 .+ λ
    a = (_Xsvd.S .* rhs)./div
    _coef = _Xsvd.V * a
    coef = vec(_coef./Xscale')
    return (;coef,LW,HKB)
end

function getinitialalpha(x::Matrix,r::Vector,log::ePPRLog)
    log.debug && log("Get Initial α ...")
    α = nothing
    try
        # Faster Cholesky decomposition requires x positive definite
        λ = ridgeparam(x,r).LW
        α = ridge(x,r,λ,bias=false)
    catch
        # Slower SVD decomposition
        log.debug && log("Design Matrix not positive definite, switch to SVD method ...")
        λ = ridge_R(x,r).LW
        α = ridge_R(x,r,λ).coef
    end
    α.-=mean(α);normalize!(α,2);α
end

function refitmodelbeta!(model::ePPRModel,y::Vector,log::ePPRLog)
    log.debug && log("Refit Model βs ...")
    ml = length(model);n=length(y)
    x = Matrix{Float64}(undef,n,ml)
    for i in 1:ml
        x[:,i] = model.phivalues[i]
    end
    β = llsq(x, y .- model.ymean, bias=false)
    if log.debug
        log("Old βs: $(model.beta)")
        log("New βs: $β")
    end
    model.beta = β
    model.residuals = y .- x*β
    return model
end

function refitmodel!(model::ePPRModel,x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog)
    log.debug && log("ePPR Model Refit ...")
    model = refitmodelbeta!(model,y,log)
    for i in 1:length(model)
        oldloss = lossfun(model,y,hp)
        oldβ,oldΦ,oldα,index,oldΦvs,oldtrustregionsize = getterm(model,i)
        model.residuals .+= oldβ*oldΦvs

        t = index.t;s=index.s
        tx = delayx(x,t,hp)
        log.debug && log("Refit (Temporal-$t, Spatial-$s) New Term ...")
        β,Φ,α,Φvs,trustregionsize = fitnewterm(tx,model.residuals,oldα,hp,log,convergerate=hp.refitconvergerate,trustregionsize=oldtrustregionsize)
        setterm!(model,i,β,Φ,α,index,Φvs,trustregionsize)
        newloss = lossfun(model,y,hp)
        if newloss > oldloss
            log.debug && log("Model Loss increased from $oldloss to $newloss. Discard the new term, keep the old one.")
            setterm!(model,i,oldβ,oldΦ,oldα,index,oldΦvs,oldtrustregionsize)
            model.residuals .-= oldβ*oldΦvs
        else
            model.residuals .-= β*Φvs
        end
    end
    return model
end

function backwardstepwise(model::ePPRModel,x::Matrix,y::Vector,hp::ePPRHyperParams,log::ePPRLog)
    log.debug && log("ePPR Backward Stepwise ...")
    log.plot && log(plotmodel(model,hp),file="Model_$(length(model)) (λ=$(hp.lambda)).png")
    models=[deepcopy(model)];hp.droptermindex=[]
    for i in length(model)-1:-1:hp.mnbt
        model,dropindex = dropleastimportantterm!(model,log)
        pushfirst!(hp.droptermindex,dropindex)
        model = refitmodel!(model,x,y,hp,log)
        log.plot && log(plotmodel(model,hp),file="Model_$(length(model)) (λ=$(hp.lambda)).png")
        pushfirst!(models,deepcopy(model))
    end
    return models
end

dropleastimportantterm!(model::ePPRModel,log::ePPRLog)=dropterm!(model,argmin(abs.(model.beta)),log)
function dropterm!(model::ePPRModel,i::Int,log::ePPRLog)
    dropindex = model.index[i]
    β=model.beta[i]
    log.debug && log("Drop Term: (temporal-$(dropindex.t), spatial-$(dropindex.s)) with β: $(β).")
    deleteat!(model,i)
    return model,dropindex
end

lossfun(g::Vector) = 0.5*norm(g,2)^2
"""
Loss function for a model term

``f(α) = sum((r.-Φ(x*α)).^2) + λ*norm(hp.alphapenaltyoperator*α,2)^2``
"""
lossfun(r::Vector,x::Matrix,α::Vector,Φ,hp::ePPRHyperParams) = lossfun([r.-Φ(x*α);sqrt(hp.lambda)*hp.alphapenaltyoperator*α])
"Loss function for model terms"
function lossfun(model::ePPRModel,y::Vector,hp::ePPRHyperParams)
    modelloss = lossfun(y.-model())
    penaltyloss = 0.5*hp.lambda*sum(norm.([hp.alphapenaltyoperator].*model.alpha,2).^2)
    return modelloss + penaltyloss
end

(phi::RObject)(x) = rcopy(R"predict($phi, x=$x)$y")
"Alternately fit `Φ` and `α` using smooth spline and newton trust region"
function fitnewterm(x::Matrix,r::Vector,α::Vector,hp::ePPRHyperParams,log::ePPRLog;convergerate::Float64=hp.forwardconvergerate,trustregionsize::Float64=hp.trustregioninitsize)
    saturatediteration = 0;xa=phi=Φvs=nothing
    for i in 1:hp.newtermmaxiteration
        xa = x*α
        phi = R"smooth.spline(y=$r, x=$xa, df=$(hp.phidf), spar=NULL, cv=FALSE)"
        Φvs = phi(xa)
        gt = r.-Φvs;gp = sqrt(hp.lambda)*hp.alphapenaltyoperator*α;g = [gt;gp]
        log.debug && log("New Term $(i)th iteration. TermLoss: $(lossfun(gt)), PenaltyLoss: $(lossfun(gp)).")
        f(a) = lossfun(r,x,a,phi,hp)
        # Loss before trust region
        lossₒ = lossfun(g)
        Φ′vs = rcopy(R"predict($phi, x=$xa, deriv=1)$y")
        gg = [-Φ′vs.*x;sqrt(hp.lambda)*hp.alphapenaltyoperator]'
        f′ = gg*g
        f″ = gg*gg'
        # α and Loss after trust region
        success,α,lossₙ,trustregionsize = newtontrustregion(f,α,lossₒ,f′,f″,trustregionsize,hp.trustregionmaxsize,hp.trustregioneta,hp.trustregionmaxiteration,log)
        if !success
            log.debug && log("NewtonTrustRegion failed, New Term use old α.")
            break
        end
        log.debug && lossₙ > lossₒ && log("New Term $(i)th iteration. Loss increased from $(lossₒ) to $(lossₙ).")
        cr = (lossₒ-lossₙ)/lossₒ
        if lossₙ < lossₒ && cr < convergerate
            saturatediteration+=1
            if saturatediteration >= hp.nsaturatediteration
                log.debug && log("New Term converged in $i iterations with (lossₒ-lossₙ)/lossₒ = $(cr).")
                break
            end
        else
            saturatediteration=0
        end
        log.debug && i==hp.newtermmaxiteration && log("New Term does not converge in $i iterations.")
    end
    β = std(Φvs)

    si = sortperm(xa)
    Φ = Spline1D(xa[si], Φvs[si], bc="nearest", s=0.5)
    # xknots = range(extrema(xa)...,length=21)[2:end-1]
    # Φ = Spline1D(xa[si], Φvs[si], xknots, bc="nearest")

    return β,Φ,α,Φvs/β,trustregionsize
end

function fitnewterm(x::Matrix,r::Vector,α::Vector,phidf::Int)
    xa = x*α
    phi = R"smooth.spline(y=$r, x=$xa, df=$phidf, spar=NULL, cv=FALSE)"
    Φvs = phi(xa)
    β = std(Φvs)

    si = sortperm(xa)
    Φ = Spline1D(xa[si], Φvs[si], bc="nearest", s=0.5)
    # xknots = range(extrema(xa)...,length=21)[2:end-1]
    # Φ = Spline1D(xa[si], Φvs[si], xknots, bc="nearest")

    return β,Φ,α,Φvs/β
end

"""
Update `α` only once [^1]

subproblem: Min mᵢ(p) = fᵢ + gᵢᵀp + 0.5pᵀBᵢp , ∥p∥ ⩽ rᵢ

Theorem 4.1
a: (Bᵢ + λI)pˢ = -gᵢ , λ ⩾ 0
b: λ(rᵢ - ∥pˢ∥) = 0
c: Bᵢ + λI positive definite

λ = 0 => ∥pˢ∥ ⩽ rᵢ, Bᵢpˢ = -gᵢ, Bᵢ positive definite
∥pˢ∥ = rᵢ => λ ⩾ 0, p(λ) = -(Bᵢ + λI)⁻¹gᵢ, Bᵢ + λI positive definite

[^1]

"Numerical Optimization, Nocedal and Wright, 2006"
"""
function newtontrustregion(f::Function,x₀::Vector,f₀::Float64,g₀::Vector,H₀::Matrix,r::Float64,rmax::Float64,η::Float64,maxiteration::Int,log::ePPRLog)
    eh = eigen(Symmetric(H₀))
    posdef = isposdef(eh)
    qᵀg = eh.vectors'*g₀
    if posdef
        pˢ = -eh.vectors*(qᵀg./eh.values)
        pˢₙ = norm(pˢ,2)
    end
    ehminvalue = minimum(eh.values)
    λe = eh.values.-ehminvalue
    ehminidx = λe .== 0
    C1 = sum((qᵀg./λe)[.!ehminidx].^2)
    C2 = sum(qᵀg[ehminidx].^2)
    C3 = sum(qᵀg.^2)

    for i in 1:maxiteration
        log.debug && log("NewtonTrustRegion $(i)th iteration, r = $(r)")
        # try for solution when λ = 0
        if posdef && pˢₙ <= r
            pᵢ = pˢ
            islambdazero = true
        else
            islambdazero = false
            # easy or hard-easy cases
            if C2 > 0 || C1 >= r^2
                iseasy = true
                ishard = C2==0

                λdn = sqrt(C2)/r
                λup = sqrt(C3)/r
                function ϕ(λ)
                    if λ==0
                        if C2 > 0
                            return -1/r
                        else
                            return sqrt(1/C1) - 1/r
                        end
                    end
                    return 1/norm(qᵀg./(λe.+λ),2) - 1/r
                end
                if ϕ(λup) <= 0
                    λ = λup
                elseif ϕ(λdn) >= 0
                    λ = λdn
                else
                    λ = fzero(ϕ,λdn,λup)
                end
                pᵢ = -eh.vectors*(qᵀg./(λe.+λ))
            else
                # hard-hard case
                iseasy = false
                ishard = true
                w = qᵀg./λe
                w[ehminidx]=0
                τ = sqrt(r^2-C1)
                𝑧 = eh.vectors[:,1]
                pᵢ = -eh.vectors*w + τ*𝑧
            end
        end
        # ρ: ratio of actual change versus predicted change
        xᵢ = x₀ + pᵢ
        fᵢ = f(xᵢ)
        ρ = (fᵢ - f₀) / (pᵢ'*(g₀ + H₀*pᵢ/2))
        # update trust region size
        if ρ < 0.25
            r /= 4
        elseif ρ > 0.75 && !islambdazero
            r = min(2*r,rmax)
        end
        if log.debug
            log("                                 ρ = $ρ")
            if islambdazero
                steptype="λ = 0"
            else
                if ishard
                    steptype=iseasy ? "hard-easy" : "hard-hard"
                else
                    steptype="easy"
                end
            end
            log("                                 step is $steptype")
        end
        # accept solution only once
        if ρ > η
            return true,xᵢ,fᵢ,r
        end
    end
    log.debug && log("NewtonTrustRegion does not converge in $maxiteration iterations.")
    return false,x₀,f₀,r
end

"""
2D Laplacian Filter in Matrix Form
"""
function laplacian2dmatrix(nrow::Int,ncol::Int)
    center = [-1 -1 -1;
              -1  8 -1;
              -1 -1 -1]
    firstrow=[-1  5 -1;
              -1 -1 -1;
               0  0  0]
    lastrow= [ 0  0  0;
              -1 -1 -1;
              -1  5 -1]
    firstcol=[-1 -1  0;
               5 -1  0;
              -1 -1  0]
    lastcol= [ 0 -1 -1;
               0 -1  5;
               0 -1 -1]
    topleft= [ 3 -1  0;
              -1 -1  0;
               0  0  0]
    topright=[ 0 -1  3;
               0 -1 -1;
               0  0  0]
    downleft=[ 0  0  0;
              -1 -1  0;
               3 -1  0]
    downright=[0  0  0;
               0 -1 -1;
               0 -1  3]
    lli = LinearIndices((nrow,ncol))
    lm = zeros(nrow*ncol,nrow*ncol)
    # fill center
    for r in 2:nrow-1, c in 2:ncol-1
        f=zeros(nrow,ncol)
        f[r-1:r+1,c-1:c+1]=center
        lm[lli[r,c],:]=vec(f)
    end
    for c in 2:ncol-1
        # fill first row
        r=1;f=zeros(nrow,ncol)
        f[r:r+2,c-1:c+1]=firstrow
        lm[lli[r,c],:]=vec(f)
        # fill last row
        r=nrow;f=zeros(nrow,ncol)
        f[r-2:r,c-1:c+1]=lastrow
        lm[lli[r,c],:]=vec(f)
    end
    for r in 2:nrow-1
        # fill first col
        c=1;f=zeros(nrow,ncol)
        f[r-1:r+1,c:c+2]=firstcol
        lm[lli[r,c],:]=vec(f)
        # fill last col
        c=ncol;f=zeros(nrow,ncol)
        f[r-1:r+1,c-2:c]=lastcol
        lm[lli[r,c],:]=vec(f)
    end
    # fill top left
    r=1;c=1;f=zeros(nrow,ncol)
    f[r:r+2,c:c+2]=topleft
    lm[lli[r,c],:]=vec(f)
    # fill top right
    r=1;c=ncol;f=zeros(nrow,ncol)
    f[r:r+2,c-2:c]=topright
    lm[lli[r,c],:]=vec(f)
    # fill down left
    r=nrow;c=1;f=zeros(nrow,ncol)
    f[r-2:r,c:c+2]=downleft
    lm[lli[r,c],:]=vec(f)
    # fill down right
    r=nrow;c=ncol;f=zeros(nrow,ncol)
    f[r-2:r,c-2:c]=downright
    lm[lli[r,c],:]=vec(f)
    return lm
end

## Visualization
"Plot ePPR model"
function plotmodel(model::ePPRModel,hp::ePPRHyperParams;ti = map(i->i.t+1,model.index),si = map(i->i.s,model.index),tmax=maximum(ti),smax=maximum(si),
                   color=:coolwarm,linkclim=true,xlim=200,colorbar=:none,width=200,height=200)
    plot(plotalpha(model,hp;ti,si,tmax,smax,color,linkclim,colorbar,width,height),plotphi(model,hp;ti,si,tmax,smax,xlim,width,height),
    layout=grid(2,1,heights=[hp.ndelay,1]/(hp.ndelay+1)),left_margin=4Plots.mm,size=(smax*width,tmax*height+tmax*hp.ndelay*height),link=:none)
end

"Plot ePPR α for each term"
function plotalpha(model::ePPRModel,hp::ePPRHyperParams;ti = map(i->i.t+1,model.index),si = map(i->i.s,model.index),tmax=maximum(ti),smax=maximum(si),
                   color=:coolwarm,linkclim=true,colorbar=:none,width=200,height=200)
    utsmax=Dict(t=>maximum(si[ti.==t]) for t in unique(ti));inpx=prod(hp.imagesize);xnpx=length(hp.xindex);αlim=0
    p = plot(layout=(tmax,smax),size=(smax*width,tmax*hp.ndelay*height),yflip=true,framestyle=:none,margin=-1Plots.mm)
    for i in 1:length(model)
        t=ti[i];s=si[i];iα=model.alpha[i]
        α = mapfoldl(d->begin
                            if xnpx==0
                                da = reshape(iα[(1:inpx).+d*inpx],hp.imagesize)
                            else
                                da = zeros(hp.imagesize)
                                da[hp.xindex]=iα[(1:xnpx).+d*xnpx]
                            end
                            da
                        end,(a1,a2)->[a1;a2],0:hp.ndelay-1)
        if linkclim
            αlim=max(αlim,maximum(abs.(α)))
        end
        heatmap!(p[t,s],α,color=color,ratio=:equal,colorbar=colorbar)
    end
    for t in 1:tmax
        ylabel!(p[t,1],"Delay $(hp.ndelay>1 ? "$(hp.ndelay-1) - 0" : t-1)")
    end
    linkclim && plot!(p,clims=(-αlim,αlim))
    p
end

"Plot ePPR Φ function for each term"
function plotphi(model::ePPRModel,hp::ePPRHyperParams;ti = map(i->i.t+1,model.index),si = map(i->i.s,model.index),tmax=maximum(ti),smax=maximum(si),
                 xlim=200,width=200,height=200)
    p = plot(layout=(tmax,smax),size=(smax*width,tmax*height),leg=false,grid=false,framestyle=:none,xtick=[-xlim,0,xlim],margin=-1Plots.mm,titlefontsize=10)
    for i in 1:length(model)
        t=ti[i];s=si[i]
        vline!(p[t,s],[0],linewidth=0.5,color=:grey80);hline!(p[t,s],[0],linewidth=0.5,color=:grey80)
        plot!(p[t,s],x->model.phi[i](x),-xlim,xlim,color=:dodgerblue,linewidth=3,link=:all,framestyle=:axes,title="β=$(round(model.beta[i],digits=3))")
        if t < tmax
            xaxis!(p[t,s],xformatter=_->"")
        end
        if s > 1
            yaxis!(p[t,s],yformatter=_->"")
        end
    end
    for t in 1:tmax
        ylabel!(p[t,1],"Delay $(hp.ndelay>1 ? "$(hp.ndelay-1) - 0" : t-1)")
    end
    p
end

"Plot correlation between Response and ePPR prediction"
plotcor(models::Vector{ePPRModel},mcors)=plotcor(length.(models),mcors,xlabel="Number of Terms")
function plotcor(ms,mcors;xlabel="Models")
    scatter(ms,hcat(mcors...)',ylabel="Pearson Correlation",xlabel=xlabel,leg=false,xtick=ms)
end