text:
module LRMoE

import Base: size, length, convert, show, getindex, rand, vec, inv, expm1, abs, log1p
import Base: isnan, isinf
import Base: sum, maximum, minimum, ceil, floor, extrema, +, -, *, ==
import Base: convert, copy, findfirst, summary
import Base.Math: @horner
import Base: π

using StatsFuns
import StatsFuns: log1mexp, log1pexp, logsumexp
import StatsFuns: sqrt2, invsqrt2π

using Statistics
import Statistics: quantile, mean, var, median

using Distributions
import Distributions: pdf, cdf, ccdf, logpdf, logcdf, logccdf, quantile
import Distributions: rand, AbstractRNG
import Distributions: mean, var, skewness, kurtosis
import Distributions: UnivariateDistribution, DiscreteUnivariateDistribution, ContinuousUnivariateDistribution
import Distributions: @distr_support, RecursiveProbabilityEvaluator
import Distributions: Bernoulli, Multinomial
import Distributions: Binomial, Poisson
import Distributions: Gamma, InverseGaussian, LogNormal, Normal, Weibull

using InvertedIndices
import InvertedIndices: Not

using LinearAlgebra
import LinearAlgebra: I, Cholesky

using SpecialFunctions
import SpecialFunctions: erf, loggamma, gamma_inc, gamma, beta_inc

using QuadGK
import QuadGK: quadgk

using Optim
import Optim: optimize, minimizer

using Clustering
import Clustering: kmeans, assignments, counts

using HypothesisTests
import HypothesisTests: ExactOneSampleKSTest, pvalue, ksstats

using Logging

using Roots
import Roots: find_zero, Order2

export
    # generic types
    convert,
    copy,
    ZeroInflation,
    ExpertSupport,
    AnyExpert,
    DiscreteExpert,
    ContinuousExpert,
    RealDiscreteExpert,
    RealContinuousExpert,
    NonNegDiscreteExpert,
    NonNegContinuousExpert,
    ZIDiscreteExpert,
    ZIContinuousExpert,
    NonZIDiscreteExpert,
    NonZIContinuousExpert,

    # model related
    summary,
    LRMoEModel,
    LRMoESTD,
    LRMoEFittingResult,
    LRMoESTDFit,

    # loglikelihood functions
    pdf, logpdf,
    cdf, logcdf,
    rowlogsumexp,

    expert_ll_pos,
    expert_tn_pos,
    expert_tn_bar_pos,
    expert_ll,
    expert_tn,
    expert_tn_bar,
    expert_ll_pos_list,
    expert_tn_pos_list,
    expert_tn_bar_pos_list,
    expert_ll_list,
    expert_tn_list,
    expert_tn_bar_list,

    loglik_aggre_dim,
    loglik_aggre_gate_dim,
    loglik_aggre_gate_dim_comp,
    loglik_np,
    loglik_exact,

    penalty_α,
    penalty_params,
    penalize,

    # EM related
    nan2num,
    EM_E_z_obs,
    EM_E_z_lat,
    EM_E_k,
    EM_M_α,
    EM_M_dQdα,
    EM_M_dQ2dα2,
    EM_E_z_zero_obs_update,
    EM_E_z_zero_lat_update,
    EM_E_z_zero_obs,
    EM_E_z_zero_lat,

    # init
    cmm_init,
    cmm_init_exact,

    # gating
    LogitGating,

    # experts
    Burr,
    GammaCount,

    params,
    p_zero,
    params_init,
    exposurize_expert,
    exposurize_model,

    BurrExpert, ZIBurrExpert,
    GammaExpert, ZIGammaExpert,
    InverseGaussianExpert, ZIInverseGaussianExpert,
    LogNormalExpert, ZILogNormalExpert,
    WeibullExpert, ZIWeibullExpert,

    BinomialExpert, ZIBinomialExpert,
    GammaCountExpert, ZIGammaCountExpert,
    NegativeBinomialExpert, ZINegativeBinomialExpert,
    PoissonExpert, ZIPoissonExpert,

    # fitting
    fit_main,
    fit_exact,
    fit_LRMoE,

    # simulation
    sim_expert,
    sim_logit_gating,
    sim_dataset,

    # prediction
    mean, var, quantile,
    predict_class_prior,
    predict_class_posterior,
    predict_mean_prior,
    predict_mean_posterior,
    predict_var_prior,
    predict_var_posterior,
    predict_limit_prior,
    predict_limit_posterior,
    predict_excess_prior,
    predict_excess_posterior,
    predict_VaRCTE_prior,
    predict_VaRCTE_posterior

### source files

include("utils.jl")
include("AICBIC.jl")
include("modelstruct.jl")

include("gating.jl")
include("expert.jl")
include("loglik.jl")
include("penalty.jl")

include("paramsinit.jl")
include("fit.jl")

include("simulation.jl")
include("predict.jl")

# include("experts/ll/expert_ll_pos.jl")

"""
A Julia package for the Logit-Reduced Mixture-of-Experts (LRMoE) model.
"""
LRMoE

end # module

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
abstract type Camera end

"""
    open!(camera::Camera)

Open camera.
"""
open!(camera::Camera) = error("No implementation for $(typeof(camera))")

"""
    close!(camera::Camera)

Close camera.
"""
close!(camera::Camera) = error("No implementation for $(typeof(camera))")

"""
    isopen(camera::Camera)

Return if the camera is open.
"""
isopen(camera::Camera) = error("No implementation for $(typeof(camera))")

"""
    isrunning(camera::Camera)

Return if the camera is running.
"""
isrunning(camera::Camera) = error("No implementation for $(typeof(camera))")

"""
    start!(camera::Camera)

Start camera, i.e. start image acquisition.
"""
start!(camera::Camera) = error("No implementation for $(typeof(camera))")

"""
    stop!(camera::Camera)

Stop camera, i.e. stop image acquisition.
"""
stop!(camera::Camera) = error("No implementation for $(typeof(camera))")

"""
    take!(camera::Camera)

Take an image, i.e. an [`AbstractAcquiredImage`](@ref). Blocks until an image is available.
"""
take!(camera::Camera)::AbstractAcquiredImage = error("No implementation for $(typeof(camera))")

"""
    trigger!(camera::Camera)

Trigger image acquisition.
"""
trigger!(camera::Camera) = error("No implementation for $(typeof(camera))")

ext: jl


text:
module CausalCore

using Distributions
using Statistics
export ExogenousVariable, EndogenousVariable, intervene, randomsample, prob

"""
Struct for an exogenous variable.
"""
struct ExogenousVariable
    name
    distribution
end

"""
Struct for an endogenous variable.
"""
struct EndogenousVariable
    operator
    parents
end

"""
Evalutes an endogenous variable given values of parent variables.
"""
function (u::EndogenousVariable)(u_)
    evaluated_vars = map(parent -> (parent)(u_), u.parents)
    return reduce(u.operator, evaluated_vars)
end

"""
Finds exogenous variable value from list of exogenous variable values.
"""
function (u::ExogenousVariable)(u_)
    for name in keys(u_)
        if u.name == name
            return getindex(u_, name)
        end
    end
end

"""
Randomly samples value from a variable.
"""
function randomsample(u)
    parent_values = Dict{Any,Any}()
    randomsample_helper(u, parent_values)
end

function randomsample_helper(u::EndogenousVariable, dict::Dict{Any, Any})
    if haskey(dict, u.name)
      getindex(dict, u.name)
    else
      sample_values = map(parent -> randomsample_helper(parent, dict), u.parents)
      val = u.operator(sample_values...)
      push!(dict, u.name => val)
      val
    end
end

function randomsample_helper(u::ExogenousVariable, dict::Dict{Any, Any})
    if haskey(dict, u.name)
        getindex(dict, u.name)
    else
        val = rand(u.distribution)
        push!(dict, u.name => val)
        val
    end
end

function randomsample_helper(u::Array, dict::Dict{Any, Any})
  sample_values = []
  for variable in u
    sample_value = randomsample_helper(variable, dict)
    push!(sample_values, sample_value)
  end
  sample_values
end

"""
Compute probability that u == val.
"""
function prob(u, val)
  NSAMPS = 100000
  total = 0
  for i=1:NSAMPS
    if randomsample(u) == val
      total += 1
    end
  end
  total/NSAMPS
end

"""
Constructs new endogenous variable based on an intervention of an existing variable.
@param u1: endogenous variable that undergoes intervention.
@param u2: exogenous variable that takes deterministic value under intervention.
@param val: value of exogenous variable under intervention.
"""
function intervene(u1::EndogenousVariable, u2::ExogenousVariable, val)
    new_parents = map(parent -> intervene(parent, u2, val), u1.parents)
    EndogenousVariable(u1.operator, new_parents)
end

"""
Constructs new, deterministic exogenous variable based on an intervention.
@param u1: exogenous variable that undergoes intervention.
@param u2: exogenous variable that takes deterministic value under intervention.
@param val: value of exogenous variable under intervention.
"""
function intervene(u1::ExogenousVariable, u2::ExogenousVariable, val)
    if (u1.name == u2.name)
        if val == 0
            ExogenousVariable(u1.name, Bernoulli(0))
        else
            ExogenousVariable(u1.name, Bernoulli(1))
        end
    else
        u1
    end
end

end


text:
############################################################
## joMatrix - extra functions

# elements(jo)
#elements(A::joAbstractDAparallelToggleOperator) = throw(joAbstractDAparallelToggleOperatorException("elements: pointless operation for joDAdistribute/joDAgather operations"))

# hasinverse(jo)

# issquare(jo)

# istall(jo)

# iswide(jo)

# iscomplex(jo)

# islinear(jo)

# isadjoint(jo)


text:
const Adjlist{T} = Vector{Vector{T}}

function Adjlist{T}(n::Integer) where {T}
    result = Adjlist{T}(undef, n)
    @inbounds for i in 1:n
        result[i] = T[]
    end
    return result
end

function deepcopy_adjlist(adjlist::Adjlist{T}) where {T}
    n = length(adjlist)
    result = Vector{Vector{T}}(undef, n)
    @inbounds for i in OneTo(n)
        list = adjlist[i]
        n_list = length(list)
        result_list = Vector{T}(undef, n_list)
        for j in OneTo(n_list)
            result_list[j] = deepcopy(list[j])
        end
        result[i] = result_list
    end
    return result
end

const EdgeValContainer{T} = Union{Adjlist{T},
                                  Tuple{Vararg{Adjlist}},
                                  NamedTuple{S, <: Tuple{Vararg{Adjlist}}} where S
                                 }

create_edgeval_list(nv, E_VAL::Type{Nothing}) = nothing
create_edgeval_list(nv, E_VAL::Type) = Adjlist{E_VAL}(nv)
create_edgeval_list(nv, E_VAL::Type{<:Tuple}) = Tuple(Adjlist{T}(nv) for T in E_VAL.parameters)
create_edgeval_list(nv, E_VAL::Type{<:NamedTuple}) = NamedTuple{Tuple(E_VAL.names)}(Adjlist{T}(nv) for T in E_VAL.types)

edgevals_container_type(::Val{Nothing}) where {E_VAL <: Type} = Adjlist{Nothing}
edgevals_container_type(::Val{E_VAL}) where {E_VAL <: Type} = Adjlist{E_VAL}

@generated function edgevals_container_type(::Val{E_VAL}) where {E_VAL <:Tuple}
    R = Tuple{( Adjlist{T} for T in E_VAL.types )...}
    return :($R)
end

@generated function edgevals_container_type(::Val{E_VAL}) where {E_VAL <:NamedTuple}
    R = NamedTuple{ Tuple(E_VAL.names), Tuple{( Adjlist{T} for T in E_VAL.types )...}}
    return :($R)
end

# TODO maybe move somewhere else
function set_value_for_index!(adjlist::Adjlist,
                              s::Integer,
                              index::Integer,
                              value)
    @inbounds adjlist[s][index] = value
    return nothing
end

@generated function set_value_for_index!(tup_adjlist::TupleOrNamedTuple, s::Integer, index::Integer, value::T) where {T <: TupleOrNamedTuple}
    len = length(T.types)
    exprs = Expr[]
    for i in OneTo(len)
        e = :(@inbounds tup_adjlist[$i][s][index] = value[$i])
        push!(exprs, e)
    end
    return Expr(:block, exprs...)
end

#= TODO remove
function set_value_for_index!(tup_adjlist::TupleOrNamedTuple, s::Integer, index::Integer, value::TupleOrNamedTuple)
    @inbounds for i in eachindex(value)
        tup_adjlist[i][s][index] = value[i]
    end
    return nothing
end
=#

function set_value_for_index!(tup_adjlist::TupleOrNamedTuple, s::Integer, index::Integer, key, value)
    @inbounds tup_adjlist[key][s][index] = value
    return nothing
end

function insert_value_for_index!(adjlist::Adjlist,
                                 s::Integer,
                                 index::Integer,
                                 value)
    @inbounds insert!(adjlist[s], index, value)
    return nothing
end

@generated function insert_value_for_index!(tup_adjlist::TupleOrNamedTuple,
                                 s::Integer,
                                 index::Integer,
                                 value::T) where {T <: TupleOrNamedTuple}
    len = length(T.types)
    exprs = Expr[]
    for i in OneTo(len)
        e = :(@inbounds insert!(tup_adjlist[$i][s], index, value[$i]))
        push!(exprs, e)
    end
    return Expr(:block, exprs...)
end

#= TODO remove
function insert_value_for_index!(tup_adjlist::TupleOrNamedTuple,
                                 s::Integer,
                                 index::Integer,
                                 value)
    @inbounds for i in eachindex(value)
        insert!(tup_adjlist[i][s], index, value[i])
    end
    return nothing
end
=#

function delete_value_for_index!(adjlist::Adjlist,
                                 s::Integer,
                                 index::Integer)
    @inbounds deleteat!(adjlist[s], index)
    return nothing
end

function delete_value_for_index!(tup_adjlist::TupleOrNamedTuple,
                                 s::Integer,
                                 index::Integer)
    @inbounds for adjlist in tup_adjlist
        deleteat!(adjlist[s], index)
    end
    return nothing
end

function value_for_index(adjlist::Adjlist, E_VAL, s::Integer, index::Integer)
    @inbounds return adjlist[s][index]
end

# TODO E_VAL could probably be extracted without being explicitly passed to the function
function value_for_index(tup_adjlist::TupleOrNamedTuple, E_VAL::Type, s::Integer, index::Integer)
    @inbounds return E_VAL( adjlist[s][index] for adjlist in tup_adjlist )
end

function value_for_index(tup_adjlist::TupleOrNamedTuple, E_VAL, s::Integer, index::Integer, key)
    adjlist = tup_adjlist[key]
    @inbounds return adjlist[s][index]
end

#=
# TODO binary search with hint
function searchsortedfirst_hint(v::AbstractVector, x, hint::Int)
    lo = 0
    hi = length(v) + 1
    m = hint
    @inbounds while lo < hi -1
        #=
        if v[m] < x
            lo = m
        else
            hi = m
        end
        =#
        is_less = v[m] < x
        lo = ifelse(is_less, m, lo)
        hi = ifelse(is_less, hi, m)
        m = (lo+hi)>>>1
    end
    return hi
end
=#


text:
# AllDifferent.
# Nice syntax:    @constraint(m, alldifferent(x, y, z))
# Default syntax: @constraint(m, [x, y, z] in AllDifferent(3))
function JuMP.parse_one_operator_constraint(_error::Function, ::Bool, ::Val{:alldifferent}, arg...)
  # All variables are explicit in the call. For instance: @constraint(m, alldifferent(x, y)).
  if length(arg) > 1 && typeof(arg) <: Tuple && all(typeof(arg).parameters .== Symbol)
    set = CP.AllDifferent(length(arg))
    func = Expr(:vect, arg...)

    variable, parse_code = JuMP._MA.rewrite(func)
    build_call = JuMP._build_call(_error, false, variable, set)
    return parse_code, build_call
  end

  # The number of variables is not statically known, generate some code to get it.
  # All the generated code will be run in JuMP's scope, hence the set_ variable.
  set_ = CP.AllDifferent
  variable, parse_code = JuMP._MA.rewrite(arg[1])
  build_call = :(build_constraint($_error, $variable, ($set_)(length($variable))))
  return parse_code, build_call
end

# Domain.
# Nice syntax:    @constraint(m, x in [1, 2, 3]), like Membership
# Default syntax: @constraint(m, x in Domain(Set(1, 2, 3)))
function JuMP.build_constraint(_error::Function, func::AbstractJuMPScalar, values::Vector{T}) where {T <: Real}
  return ScalarConstraint(func, CP.Domain(Set(values)))
end

# Membership.
# Nice syntax:    @constraint(m, x in [y, z]), like Domain
# Default syntax: @constraint(m, [x, y, z] in Membership(3))
# Also works with mixed variables and constants: @constraint(m, x in [y, z, 1])
function JuMP.build_constraint(_error::Function, func::AbstractJuMPScalar, values::Vector{<:JuMP.AbstractJuMPScalar})
  return VectorConstraint(vcat([func], values), CP.Membership(length(values)))
end

# DifferentFrom.
# Nice syntax:    @constraint(m, x != y)
# Default syntax: @constraint(m, x - y in DifferentFrom(0.0))
JuMP.sense_to_set(_error::Function, ::Val{:(!=)}) = CP.DifferentFrom(0.0)

# Strictly.
JuMP.sense_to_set(_error::Function, ::Val{:(<)}) = CP.Strictly(MOI.LessThan(0.0))
JuMP.sense_to_set(_error::Function, ::Val{:(>)}) = CP.Strictly(MOI.GreaterThan(0.0))

# Element.
# Nicer syntax:   @constraint(m, y == array[x])
# Nicer syntax:   @constraint(m, y == element(array, x))
# Nice syntax:    @constraint(m, element(y, array, x))
# Default syntax: @constraint(m, [y, x] in Element(array, 2))
function JuMP.parse_one_operator_constraint(_error::Function, ::Bool, ::Val{:element}, args...)
  if length(args) != 3
    error("element() constraints must have three operands: the destination, the array, the index.")
  end

  variable, parse_code_variable = JuMP._MA.rewrite(args[1])
  array = esc(args[2])
  index, parse_code_index = JuMP._MA.rewrite(args[3])
  parse_code = :($parse_code_variable; $parse_code_index)

  # TODO: Likely limitation, when passing a variable as array, modifying the array in the user code will not change the array in the constraint. Problematic or not?
  set_ = CP.Element
  build_call = :(build_constraint($_error, [$variable, $index], ($set_)($array, 2)))
  return parse_code, build_call
end

function JuMP._rewrite_expr(_error::Function, ::Val{:call}, op::Val{:element}, args...)
  # Switch to getindex to retrieve the right expression.
  if length(args) != 2
    error("element() expressions must have two operands: the array, the index.")
  end

  return :(), :(), quote; $(args[1])[$(args[2])]; end
end

function Base.getindex(v::Vector{T}, i::VariableRef) where T
  # Always returns a JuMP expression.
  m = owner_model(i)
  x = @variable(m)
  @constraint(m, element(x, v, i))
  return x
end

# Sort.
# Nicer syntax:   @constraint(m, [y1, y2] == sort([x1, x2])) TODO
# Nice syntax:    @constraint(m, sort([x1, x2], [y1, y2]))
# Default syntax: @constraint(m, [x1, x2, y1, y2] in Sort(2))
function JuMP.parse_one_operator_constraint(_error::Function, ::Bool, ::Val{:sort}, arg...)
  if length(arg) != 2
    error("sort() constraints must have two operands: the array to sort, its elements in sorted order.")
  end

  # Parse the inputs.
  original, parse_code_original = JuMP._MA.rewrite(arg[1])
  destination, parse_code_destination = JuMP._MA.rewrite(arg[2])
  parse_code = :($parse_code_original; $parse_code_destination)

  # Check whether the arrays have all the same size.
  check_code = quote
    if length($original) != length($destination)
      error("Expected to have arrays of the same size, but the array to sort has size $(length($original)) " *
            "and the output sorted array has size $(length($destination)).")
    end
  end
  parse_code = :($parse_code; $check_code)

  # Generate the constraint.
  set_ = CP.Sort
  build_call = :(build_constraint($_error, vcat($original, $destination), ($set_)(length($original))))
  return parse_code, build_call
end

# SortPermutation.
# Nicer syntax:   @constraint(m, [y1, y2] == sortpermutation([x1, x2])) TODO
# Nice syntax:    @constraint(m, sortpermutation([x1, x2], [z1, z2]))
# Nice syntax:    @constraint(m, sortpermutation([x1, x2], [y1, y2], [z1, z2]))
# Default syntax: @constraint(m, [x1, x2, y1, y2, z1, z2] in SortPermutation(2))
# I.e. the sorted values are not available with the nicer syntax, and should be retrieved through Element.
function JuMP.parse_one_operator_constraint(_error::Function, ::Bool, ::Val{:sortpermutation}, arg...)
  if length(arg) != 2 && length(arg) != 3
    error("sortpermutation() constraints must have two or three operands: the array to sort, optionally its elements in sorted order, the sorting permutation.")
  end

  # Parse the inputs.
  original, parse_code_original = JuMP._MA.rewrite(arg[1])
  array1, parse_code_array1 = JuMP._MA.rewrite(arg[2])
  parse_code = :($parse_code_original; $parse_code_array1)
  if length(arg) == 3
    array2, parse_code_array2 = JuMP._MA.rewrite(arg[3])
    parse_code = :($parse_code; $parse_code_array2)
  else
    array2 = array1
    array1 = gensym()
    create_array1 = quote
      $array1 = VariableRef[]
      for i in 1:length($original)
        push!($array1,
              add_variable(owner_model(first($original)),
                           build_variable($_error,
                                          VariableInfo(false, NaN, false, NaN, false, NaN, false, NaN, false, false)),
                           ""))
      end
    end
    parse_code = :($parse_code; $create_array1)
  end

  # Check whether the arrays have all the same size.
  if length(arg) == 2
    check_code = quote
      if length($original) != length($array1)
        error("Expected to have arrays of the same size, but the array to sort has size $(length($original)) " *
              "and the output array for the permutation has size $(length($array1)).")
      end
    end
  else
    check_code = quote
      if length($original) != length($array1)
        error("Expected to have arrays of the same size, but the array to sort has size $(length($original)) " *
              "and the output sorted array has size $(length($array1)).")
      end
      if length($original) != length($array2)
        error("Expected to have arrays of the same size, but the array to sort has size $(length($original)) " *
              "and the output array for the permutation has size $(length($array2)).")
      end
    end
  end
  parse_code = :($parse_code; $check_code)

  # Generate the constraint.
  set_ = CP.SortPermutation
  build_call = :(build_constraint($_error, vcat($original, $array1, $array2), ($set_)(length($original))))
  return parse_code, build_call
end

# BinPacking.
# Nicer syntax:   @constraint(m, [load1, load2, assigned1, assigned2] == binpacking([size1, size2])) TODO
# Nice syntax:    @constraint(m, binpacking([size1, size2], [load1, load2], [assigned1, assigned2]))
# Default syntax: @constraint(m, [load1, load2, assigned1, assigned2, size1, size2] in BinPacking(2, 2))

# CapacitatedBinPacking.
# Nicer syntax:   @constraint(m, [load1, load2, assigned1, assigned2] == binpacking([size1, size2], [capa1, capa2])) TODO
# Nice syntax:    @constraint(m, binpacking([size1, size2], [load1, load2], [assigned1, assigned2], [capa1, capa2]))
# Default syntax: @constraint(m, [load1, load2, assigned1, assigned2, size1, size2, capa1, capa2] in CapacitatedBinPacking(2, 2))
function JuMP.parse_one_operator_constraint(_error::Function, ::Bool, ::Val{:binpacking}, arg...)
  if length(arg) != 3 && length(arg) != 4
    error("binpacking() constraints must have three or four operands: the bin loads, the assignment for each item, the item sizes, and optionnally the bin capacities.")
  end

  # Parse the inputs.
  load, parse_code_load = JuMP._MA.rewrite(arg[1])
  assign, parse_code_assign = JuMP._MA.rewrite(arg[2])
  size, parse_code_size = JuMP._MA.rewrite(arg[3])
  parse_code = :($parse_code_load; $parse_code_assign; $parse_code_size)

  if length(arg) == 4
    capa, parse_code_capa = JuMP._MA.rewrite(arg[4])
    parse_code = :($parse_code; $parse_code_capa)
  end

  # Check whether the arrays have all the same size.
  check_code = quote
    if length($assign) != length($size)
      error("Expected to have arrays of the same size, but the item-assignment array has size $(length($assign)) " *
            "and the item-size array has size $(length($size)).")
    end
  end
  parse_code = :($parse_code; $check_code)

  if length(arg) == 4
    check_code = quote
      if length($load) != length($capa)
        error("Expected to have arrays of the same size, but the bin-load array has size $(length($load)) " *
              "and the bin-capacity array has size $(length($capa)).")
      end
    end
    parse_code = :($parse_code; $check_code)
  end

  # Generate the constraint.
  if length(arg) == 3
    set_ = CP.BinPacking
    build_call = :(build_constraint($_error, vcat($load, $assign, $size), ($set_)(length($load), length($assign))))
  elseif length(arg) == 4
    set_ = CP.CapacitatedBinPacking
    build_call = :(build_constraint($_error, vcat($load, $assign, $size, $capa), ($set_)(length($load), length($assign))))
  end
  return parse_code, build_call
end

function JuMP.parse_ternary_constraint(_error::Function, vectorized::Bool, a, b::Val, c, d::Val, e)
  # Hacky as hell...
  b_unboxed = typeof(b).parameters[1]
  d_unboxed = typeof(d).parameters[1]
  if applicable(JuMP.parse_one_operator_constraint, _error, vectorized, Val(a), b_unboxed, c, d_unboxed, e)
    parse_code, build_call = JuMP.parse_one_operator_constraint(_error, vectorized, Val(a), b_unboxed, c, d_unboxed, e)
    return parse_code, :(), :(), build_call
  else
    # If there is no corresponding method, return the standard error from JuMP.
    JuMP.parse_ternary_constraint(_error, nothing)
  end
end

# Count.
# Nicer syntax:   @constraint(m, y == count([x1, x2, x3] .== 1.0)) TODO
# Nice syntax:    @constraint(m, y == count(1.0, x1, x2, x3)) TODO
# Default syntax: @constraint(m, [y, x1, x2, x3] in Count(1.0, 3))

# CountDistinct.
# Nice syntax:    @constraint(m, y == countdistinct(x1, x2, x3)) TODO
# Default syntax: @constraint(m, [y, x1, x2, x3] in CountDistinct(3))

# ReificationSet.

function _build_reified_constraint(
  _error::Function, variable::AbstractVariableRef,
  constraint::ScalarConstraint, ::Type{MOI.IndicatorSet{A}}) where A
  # TODO.
  set = MOI.IndicatorSet{A}(moi_set(constraint))
  return VectorConstraint([variable, jump_function(constraint)], set)
end

function JuMP.parse_constraint(_error::Function, ::Val{:(:=)}, arg...)
  if length(arg) != 2
    error("A reification constraint must have the following form: variable := { constraint }")
  end
  lhs, rhs = arg

  # Parse the left-hand variable.
  if typeof(lhs) != Symbol
    error("A reification constraint must have the following form: variable := { constraint }. The left-hand expression is not a single variable, but rather `$(arg[1])`.")
  end

  # Parse the right-hand constraint.
  if !Meta.isexpr(rhs, :braces) || length(rhs.args) != 1
      _error("Invalid right-hand side `$(rhs)` of reification constraint. Expected constraint surrounded by `{` and `}`.")
  end

  rhs_con = rhs.args[1]
  rhs_vectorized, rhs_parse_code, rhs_build_call = parse_constraint(_error, Val(rhs_con.head), rhs_con.args...)
  if rhs_vectorized
      _error("Reified constraint cannot be vectorized.")
  end

  # Build the reified constraint.
  set_ = CP.ReificationSet
  reified_build_call = quote
    return VectorConstraint(vcat([$(esc(lhs))], jump_function($rhs_build_call)), ($set_)(moi_set($rhs_build_call)))
  end

  return rhs_parse_code, reified_build_call
end


text:
#do some pre-filtering to reduce the processing time
function prefiltercrsp!(crsp::DataFrame, validpermnos::Vector{Int})::Nothing
  #rename columns to lower case

  crsp.tokeep = trues(size(crsp,1))

  crspgrp::GroupedDataFrame = groupby(crsp, :permno)
  @mpar for i ∈ 1:length(crspgrp)
    scrsp::SubDataFrame = crspgrp[i]

    if scrsp.permno[1] ∉ validpermnos
      scrsp.tokeep .= false
    end
  end

  filter!(r::DataFrameRow->r.tokeep, crsp)
  select!(crsp, Not([:tokeep]))

  return nothing

end

@inline normalizecrspnames!(crsp::DataFrame) = rename!(crsp, (s::Symbol->Symbol(lowercase(string(s)))).(names(crsp)))

#helper to handle default case
#prepcrsp(param::AbstractDict; args...) = prepcrsp(x->x, param; args...)::DataFrame

function prepcrsp(op!::Function = (noop(x)=x);
  crsptype = :not_set,
  refreshcrsp=true,
  incsvstream::Function = IN_CSV_STREAM,
  inbinstream::Function = IN_BIN_STREAM,
  outbinstream::Function = OUT_BIN_STREAM,
  csvextension::String = CSV_EXTENSION,
  datapath::String = PARAM[:datapath],
  yearrange::Union{AbstractRange, Nothing} = nothing,
  validpermnos::Union{Nothing, Vector{Int}} = nothing,
  prefix::String = "",
  crspcolumns::Vector{Symbol} = error("crsp columns are required"))::DataFrame

  crsptype ∈ [:monthly, :daily] || error("Set crsptype to :monthly or :daily instead of $crsptype")

  local crspname
  local crspparts

  if crsptype == :monthly
    crspname = PARAM[:crspmname]
    crspparts = PARAM[:crspmparts]
  elseif crsptype==:daily
    crspname = PARAM[:crspdname]
    crspparts = PARAM[:crspdparts]
  else
    @assert false
  end

  local p::String = "$datapath\\$(prefix)$crspname"
  local crsp::DataFrame
  local crspdateformat::DateFormat = PARAM[:wrdsdateformat]

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshcrsp
    #do the below so we can pre-filter

    pin::String = "$datapath\\$(crspparts[1])"
    crsp = CSV.read(incsvstream("$pin.$csvextension"), threaded=PARAM[:csvthreaded]) |> DataFrame
    normalizecrspnames!(crsp)
    (!isnothing(validpermnos)) && (prefiltercrsp!(crsp, validpermnos))

    if length(crspparts) > 1
      #try to do this is in a space efficient way
      for crsppart ∈ crspparts[2:end]
        pin = "$datapath\\$crsppart"
        local addedcrsp::DataFrame = CSV.read(
          incsvstream("$pin.$CSV_EXTENSION"), threaded=PARAM[:csvthreaded]) |> DataFrame
        normalizecrspnames!(addedcrsp)
         (!isnothing(validpermnos)) && (prefiltercrsp!(addedcrsp, validpermnos))
        append!(crsp, addedcrsp)
      end
    end

    #filter by year first
    crsp.date = (s->parsecrsp(Date, s, crspdateformat)).(crsp.date)
    (!isnothing(yearrange)) && filter!(r::DataFrameRow->year(r.date) ∈ yearrange, crsp) #restrict to desired years

    crsp = preprocesscrsp(crsp, crsptype)
    op!(crsp)
    select!(crsp, crspcolumns)

    outbinstream("$p.$BIN_EXTENSION", crsp)

  else
    crsp = inbinstream("$p.$BIN_EXTENSION")
  end

  @info "CRSP data loaded and/or saved into file $(prefix)$crspname.$BIN_EXTENSION"

  return crsp
end

function computecrspreturns(crsp::DataFrame)
  crspparsedmissing::Vector{Float64} = PARAM[:crspparsedmissing]
  retlisted::Vector{MFloat64} = (s::MString->parsecrsp(Float64, s, crspparsedmissing)).(crsp.ret)
  dlret::Vector{MFloat64} = (s::MString->parsecrsp(Float64, s, crspparsedmissing)).(crsp.dlret)

  retlisted = (r->
    !ismissing(r) && (r > -1.0) ? r : missing).(retlisted)
  dlret .= (r->
    !ismissing(r) && (r > -1.0) ? r : missing).(dlret)

  ret::Vector{MFloat64} = ((retlisted,dlret)->
    ismissing(retlisted) && ismissing(dlret) ? missing :
      ismissing(dlret) ? retlisted :
      ismissing(retlisted) ? dlret :
      (1+retlisted)*(1+dlret) - 1.
    ).(retlisted,dlret)

  return ret
end

function computecrspprices(crsp::DataFrame, crsptype::Symbol)
  #Now work on the price
  #in this case, use the mean of the two values
  #price cannot be zero
  crsp.prc .= (
    prc->coalesce(iszero(prc),true) ? missing : prc).(crsp.prc)
  crsp.dlprc .= (
    dlprc->coalesce(iszero(dlprc),true) ? missing : dlprc).(crsp.dlprc)

  @inline ismissingorzero(p::Float64) = iszero(p) ? true : false
  @inline ismissingorzero(::Missing) = true

  price = Vector{MFloat64}(undef, size(crsp,1))
  @mpar for i ∈ 1:length(price)
    if ismissingorzero(crsp.dlprc[i])
      price[i] = crsp.prc[i]
    elseif ismissingorzero(crsp.prc[i])
      price[i] = crsp.dlprc[i]
    else
      price[i] = (crsp.prc[i]+crsp.dlprc[i])/2
    end
  end

  #add in the alternative price if nothing else is available
  #(data available monthly only)
  if crsptype == :monthly
    @mpar for i ∈ 1:length(price)
      if ismissing(price[i]) && (!ismissingorzero(crsp.altprc[i]))
        price[i] = crsp.altprc[i]
      end
    end
  end

  #negative values computed via bid-ask spread method
  price .= (abs).(price)

  return price
end

function computecrspmktcap(crsp::DataFrame)

  #adjust the units and compute the market cap
  mktcap::Vector{MInt} = Vector{MInt}((round).(crsp.price .* crsp.shrout .* 1000))
  mktcap .= (i::MInt->coalesce(iszero(i),true) ? missing : i).(mktcap)

  return mktcap
end

function preprocesscrsp(crsp::DataFrame, crsptype::Symbol)::DataFrame
  crspdateformat::DateFormat = PARAM[:wrdsdateformat]

  local scrsps::GroupedDataFrame

  #NOTE: SHRCD==10,11 : US Common Shares
  #NOTE:  :EXCHCD ∈ [1,2,3] : NYSE, Amex, NASDAQ
  deleterows!(crsp, (!).(completecases(crsp[:,[:shrcd, :exchcd]])))
  local shrcdvalid::Set{Int} = Set{Int}((10,11))
  local exchcdvalid::Set{Int} = Set{Int}((1,2,3))
  filter!(r::DataFrameRow->
    (r.shrcd ∈ shrcdvalid) && (r[:exchcd] ∈ exchcdvalid), crsp)

  #Doing this before the dedup eliminates rows without returns,
  #prices and shares outstanding first. Implicitly requires that
  #these fields exist
  crsp.ret = computecrspreturns(crsp)
  select!(crsp, Not([:dlret]))
  crsp.lret = (log).(1.0 .+ crsp.ret)

  crsp.price = computecrspprices(crsp, crsptype)
  select!(crsp, Not([:prc, :dlprc]))
  crsp.mktcap = computecrspmktcap(crsp)

  #now dedup crsp
  dedupctr::Int = 0
  crsp.keep = trues(size(crsp,1))

  scrsps = groupby(crsp, [:permno, :date])
  @mpar for i ∈ 1:length(scrsps)
    local scrsp::SubDataFrame = scrsps[i]
    if size(scrsp,1) > 1
      #first try to pick records by number missing
      nmissing::Vector{Int} = zeros(Int, size(scrsp,1))
      for f ∈ usedcrspfields
        nmissing .+= (ismissing).(scrsp[!,f])
      end
      minmissing::Int = minimum(nmissing)
      scrsp.keep .= nmissing .== minmissing

      #keep the top record iff it is identical to the others
      if sum(scrsp.keep) > 1
        sscrsp = view(scrsp, scrsp.keep, :)
        nchecks::Int = size(sscrsp,1)
        keep1::Bool = true
        for i ∈ 2:nchecks
          for f ∈ usedcrspfields
            #checks if the given used field is equal to the top row
            (coalesce(sscrsp[i,f],-99.0) ≠ coalesce(sscrsp[1,f], -99.0)) && (keep1 = false)
          end
        end

        sscrsp.keep .= false
        sscrsp.keep[1] = keep1
      end
    end
  end

  deleterows!(crsp, (!).(crsp.keep))

  return crsp
end

#parse a value where the type is known: Generic case
function parsecrsp(::Type{T}, s::String, crspparsedmissing::AbstractVector{Float64}) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  (!ismissing(v)) && (v ∈ crspparsedmissing) && (v=missing)

  return v
end

#the date case
parsecrsp(::Type{Date}, s::String, crspdateformat::DateFormat) = Dates.Date(s, crspdateformat)

#helper methods and special cases
parsecrsp(::Type{<:Any}, v::Missing, flotsam) = missing
parsecrsp(::Type{Date}, i::Int, crspdateformat::DateFormat) = parsecrsp(Date, "$i", crspdateformat)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using RiceBPH.Model
using RiceBPH.Model.Agents
using InteractiveUtils

isbph(x) = true
food(model) = count(≥(0.5), model.food)
adata = [(isbph, count)]
mdata = [food]

function run_with_seed(seed)
    crop = Model.gencrop_3x3()
    model = Model.init_model(crop, 200, :corner, 0.2f0; seed=seed)
    return adf, mdf = run!(
        model, Model.agent_step!, Model.model_step!, 2880; adata=adata, mdata=mdata
    )
end

@testset "Basic" begin
    crop = Model.gencrop_3x3()
    @test true
    model = Model.init_model(crop, 200, :corner, 0.2f0; seed=1)
    @test true
    @code_warntype run!(model, Model.agent_step!, Model.model_step!, 1; adata=adata, mdata=mdata)
    @test true
end

# Before adding type parameter
# BenchmarkTools.Trial:
#   memory estimate:  20.73 MiB
#   allocs estimate:  82035
#   --------------
#   minimum time:     309.505 ms (0.00% GC)
#   median time:      313.172 ms (0.00% GC)
#   mean time:        312.776 ms (0.46% GC)
#   maximum time:     316.794 ms (1.13% GC)
#   --------------
#   samples:          16
display(@benchmark run_with_seed(1))
@testset "Benchmark" begin
    @test true
end

@testset "Reproducible" begin
    for seed in rand(UInt8, 3)
        adf1, mdf1 = run_with_seed(seed)
        adf2, mdf2 = run_with_seed(seed)
        for name in names(adf1)
            @test all(adf1[!, name] .=== adf2[!, name])
        end
        for name in names(mdf1)
            @test all(mdf1[!, name] .=== mdf2[!, name])
        end
    end
end

@testset "Generate Video" begin
    Model.video(joinpath("..", "src", "maps", "020-1x2.csv"), 200, :corner, 0.075; seed=9)
    @test true
end

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
using SketchyCGAL
using LinearAlgebra, SparseArrays
using BSON
using Plots

include("../utils.jl")

## Timing Figure
#load data
trials_size = BSON.load(joinpath(@__DIR__, "output/trials_size.bson"))
num_trials = length(trials_size)

trial_names = Vector{String}(undef, num_trials)
ns = zeros(num_trials)
num_nonzeros = zeros(num_trials)
time_1000_iter = zeros(num_trials)
storage = zeros(num_trials)

for (ind, trial) in enumerate(trials_size)
    trial_names[ind] = trial[1]
    ns[ind] = trial[2]["dim"].n
    num_nonzeros[ind] = trial[2]["dim"].nonzeros
    time_1000_iter[ind] = trial[2]["stats"].value.log.time_sec[end] / 2
    storage[ind] = trial[2]["stats"].bytes
end

size_plt_time_edge = plot(
    num_nonzeros,
    time_1000_iter,
    seriestype = :scatter,
    title="Iteration Time vs. Graph Edge Count",
    ylabel="Time per 1,000 Iterations (s)",
    xlabel="Number of Graph Edges",
    legend=false,
    dpi=300,
    yscale=:log,
    xscale=:log,
    markersize=6,
    xlims=(10^3.5, 10^5)
    # xticks=[10^3, 10^4, 10^5, 10^6]
)

size_plt_time_vert = plot(
    ns,
    time_1000_iter,
    seriestype = :scatter,
    title="Iteration Time vs. Graph Vertex Count",
    ylabel="Time per 1,000 Iterations (s)",
    xlabel="Number of Graph Vertices (side dimension)",
    legend=false,
    dpi=300,
    yscale=:log,
    xscale=:log,
    xlims=(10^2.5, 10^4.5),
    markersize=6
    # xticks=[10^3, 10^4, 10^5, 10^6]
)

size_plt_storage_edge = plot(
    num_nonzeros,
    storage/1e6,
    seriestype = :scatter,
    title="Storage vs. Graph Edge Count",
    ylabel="Total Storage (MB)",
    xlabel="Number of Graph Edges",
    legend=false,
    dpi=300,
    yscale=:log,
    xscale=:log,
    markersize=6,
    xlims=(10^3.5, 10^5),
    yticks=([16, 32, 64, 128], ["16MB", "32MB", "64MB", "128MB"]),
    ylim=(16, 64)
)

size_plt_storage_vert = plot(
    ns,
    storage/1e6,
    seriestype = :scatter,
    title="Storage vs. Graph Vertex Count",
    ylabel="Total Storage (MB)",
    xlabel="Number of Graph Vertices (side dimension)",
    legend=false,
    dpi=300,
    yscale=:log,
    xscale=:log,
    xlims=(10^2.5, 10^4.5),
    markersize=6,
    yticks=([16, 32, 64, 128], ["16MB", "32MB", "64MB", "128MB"]),
    ylim=(16, 64)
    # xticks=[10^3, 10^4, 10^5, 10^6]
)

savefig(size_plt_time_edge, joinpath(@__DIR__, "figures/size_time_edge"))
savefig(size_plt_time_vert, joinpath(@__DIR__, "figures/size_time_vert"))
savefig(size_plt_storage_edge, joinpath(@__DIR__, "figures/size_storage_edge"))
savefig(size_plt_storage_vert, joinpath(@__DIR__, "figures/size_storage_vert"))

## Setup data
G = graph_from_file(joinpath(dirname(@__DIR__), "data/gset/G72"))
n = size(G, 1)
C = -0.25*(Diagonal(G*ones(n)) - G)
b = ones(n)
# Scaling variables -- so trace is bounded by 1
scale_C = 1 / norm(C)
scale_X = 1 / n
# pstar = -7744.42783 * scale_X * scale_C #G67
pstar = -7808.53427 * scale_X * scale_C #G72

## R parameter figure
trials_R = BSON.load(joinpath(@__DIR__, "output/trials_R_G72.bson"))
Rs = [10, 25, 50, 100]

T = 10_000
obj_vals_iter = zeros(T, length(Rs) + 1)
p_infeas_iter = zeros(T, length(Rs) + 1)

# trials_R["iter"][10].value.log

for (ind, R) in enumerate(Rs)
    obj_vals_iter[:,ind] .= trials_R["iter"][R].value.log.obj_val_Xhat
    p_infeas_iter[:, ind] .= trials_R["iter"][R].value.log.primal_infeas_Xhat
end
obj_vals_iter[:,end] .= trials_R["iter"][10].value.log.obj_val
p_infeas_iter[:,end] .= trials_R["iter"][10].value.log.primal_infeas

labels = ["R = 10" "R = 25" "R = 50" "R = 100" "Full (implicit)"]

obj_vals_iter[obj_vals_iter .== 0] .= NaN
subopt = abs.(obj_vals_iter .- pstar) / abs(pstar)
subopt_plt = plot(
    1:T,
    subopt,
    xaxis=:log,
    yaxis=:log,
    labels=labels,
    lw=2,
    title="Suboptimality vs. Iteration",
    ylabel="Suboptimality",
    xlabel="Iteration",
    legend=:bottomleft,
    dpi=300
)

p_infeas_iter[p_infeas_iter .== 0] .= NaN
p_infeas_iter[:,end] .*= scale_X
p_infeas_plt = plot(
    1:T,
    p_infeas_iter,
    xaxis=:log,
    yaxis=:log,
    labels=["R=10" "R=25" "R=50" "R=100" "Full (implicit)"],
    lw=2,
    title="Primal Infeasibility vs. Iteration",
    ylabel="Primal Infeasibility",
    xlabel="Iteration",
    legend=:bottomleft,
    dpi=300
    # xlim=(15, T)
)

timestamp = zeros(T, length(Rs))
p_infeas_time = zeros(T, length(Rs))
for (ind, R) in enumerate(Rs)
    timestamp[:,ind] .= trials_R["time"][R].value.log.time_sec
    p_infeas_time[:, ind] .= trials_R["iter"][R].value.log.primal_infeas_Xhat
end

p_infeas_time[p_infeas_time .== 0] .= NaN
p_infeas_time_plt = plot(
    timestamp,
    p_infeas_iter,
    xaxis=:log,
    yaxis=:log,
    labels=["R=10" "R=25" "R=50" "R=100"],
    lw=2,
    title="Primal Infeasibility vs. Time",
    ylabel="Primal Infeasibility",
    xlabel="Time (s)",
    legend=:bottomleft,
    dpi=300
    # xlim=(minimum(timestamp[16:end,:]), 100)
)

subopt_time_plt = plot(
    timestamp,
    subopt,
    xaxis=:log,
    yaxis=:log,
    labels=labels,
    lw=2,
    title="Suboptimality vs. Time",
    ylabel="Suboptimality",
    xlabel="Time",
    legend=:bottomleft,
    dpi=300
    # xlim=(minimum(timestamp[16:end,:]), 100)
)

savefig(p_infeas_time_plt, joinpath(@__DIR__, "figures/infeas_R_time_g72"))
savefig(p_infeas_plt, joinpath(@__DIR__, "figures/infeas_R_iter_g72"))
savefig(subopt_time_plt, joinpath(@__DIR__, "figures/subopt_R_time_g72"))
savefig(subopt_plt, joinpath(@__DIR__, "figures/subopt_R_iter_g72"))

## Weighted gradient figures

trials_weights = BSON.load(joinpath(@__DIR__, "output/trials_weights_G72.bson"))
# long_run = BSON.load(joinpath(@__DIR__, "output/long_run.bson"))
# long_run["data"].value

δ_keys = [
    "std",
    "sa",
    "ea_0.01",
    "ea_0.05",
    "ea_0.1",
    # "ea_0.3",
    "ea_0.5",
    # "ea_0.8",
    "ea_0.9",
    # "ea_0.95",
    # "ea_0.99",
]

T = 10_000
num_trials = length(δ_keys)
obj_val = zeros(T, num_trials)
p_infeas = zeros(T, num_trials)

for (ind, k) in enumerate(δ_keys)
    obj_val[:,ind] .= trials_weights[k].value.log.obj_val
    p_infeas[:,ind] .= trials_weights[k].value.log.primal_infeas * scale_X
    ind += 1
end

# labels = [
#     "Standard",
#     "SA",
#     "EA δ = 0.01",
#     "EA δ = 0.05",
#     "EA δ = 0.1",
#     "EA δ = 0.3",
#     "EA δ = 0.5",
#     "EA δ = 0.8",
#     "EA δ = 0.9",
#     "EA δ = 0.95",
#     "EA δ = 0.99",
# ]
# labels = ["Standard" "SA" "EA δ = 0.01" "EA δ = 0.05" "EA δ = 0.1" "EA δ = 0.3" "EA δ = 0.5" "EA δ = 0.8" "EA δ = 0.9" "EA δ = 0.95" "EA δ = 0.99"]##
labels = ["Standard" "SA" "EA δ = 0.01" "EA δ = 0.05" "EA δ = 0.1" "EA δ = 0.5" "EA δ = 0.9"]##

subopt = abs.(obj_val .- pstar)/abs(pstar)
# labels = ["Standard" "EA" "EA2" "EA3" "SA" "LA" "UA"]
# linestyle = [:dash :dash :dash :dash :dash :dash]
subopt_plt = plot(
    1:T,
    subopt,
    xaxis=:log,
    yaxis=:log,
    labels=labels,
    lw=2,
    title="Suboptimality vs. Iteration",
    ylabel="Suboptimality",
    xlabel="Iteration",
    # linestyle=linestyle,
    legend=:bottomleft,
    dpi=300
)

infeas_plt = plot(
    1:T,
    p_infeas,
    xaxis=:log,
    yaxis=:log,
    labels=labels,
    lw=2,
    title="Primal Infeasibility vs. Iteration",
    ylabel="Primal Infeasibility",
    xlabel="Iteration",
    # linestyle=linestyle,
    legend=:bottomleft,
    dpi=300
)

savefig(infeas_plt, joinpath(@__DIR__, "figures/infeas_weights_g72"))
savefig(subopt_plt, joinpath(@__DIR__, "figures/subopt_weights_g72"))

##
trial_big = BSON.load(joinpath(@__DIR__, "output/big_sdp.bson"))

p_infeas_iter = trial_big["iter"].value.log.primal_infeas
p_infeas_xhat_iter = trial_big["iter"].value.log.primal_infeas_Xhat
obj_val_iter = trial_big["iter"].value.log.obj_val
cut_val_iter = trial_big["iter"].value.log.cut_value
timestamp = trial_big["time"].value.log.time_sec
obj_val_iter_xhat = trial_big["iter"].value.log.obj_val_Xhat

T = length(timestamp)
cut_val_iter[cut_val_iter.==0] .= NaN
obj_val_iter_xhat[obj_val_iter_xhat.==0] .= NaN
max_cut_plt = plot(
    1:T,
    [
        cut_val_iter,
        -obj_val_iter*1/scale_X,
        -obj_val_iter_xhat*1/scale_X
    ],
    xaxis=:log,
    # yaxis=:log,
    lw=2,
    labels=["MAXCUT" "SDP Objective Full (implicit)" "SDP Objective R = 10"],
    title="MAXCUT Value vs. Iteration",
    ylabel="MAXCUT Value (Randomized Rounding)",
    xlabel="Iteration",
    legend=:topright,
    dpi=300
)

max_cut_time_plt = plot(
    timestamp,
    [
        cut_val_iter,
        -obj_val_iter*1/scale_X,
        -obj_val_iter_xhat*1/scale_X
    ],
    xaxis=:log,
    # yaxis=:log,
    lw=2,
    labels=["MAXCUT" "SDP Objective Full (implicit)" "SDP Objective R = 10"],
    title="MAXCUT Value vs. Time",
    ylabel="MAXCUT Value (Randomized Rounding)",
    xlabel="Time (s)",
    legend=:topright,
    dpi=300
)

savefig(max_cut_plt, joinpath(@__DIR__, "figures/big_sdp_maxcut_iter"))
savefig(max_cut_time_plt, joinpath(@__DIR__, "figures/big_sdp_maxcut_time"))

##
# -obj_val_iter[end]*1/scale_X - cut_val_iter[end]


text:
"""

	X_1 ∪ X_2 ∪ ... ∪ X_n = {1, ..., p}

"""

#TODO : supprimer cette structure
struct NDisjointUnion <: Constraint

	X::Vector{Variable}
	domain::Vector{Int}

end

function NDisjointUnion(m::Model, X::Vector{Variable}, domain::UnitRange{Int})
	return NDisjointUnion(m, X, collect(domain))
end

function NDisjointUnion(m::Model, X::Vector{Variable}, domain::Vector{Int})

	#push!(m.varsInter, Variable(1:m.p, "F[$(length(m.varsInter)+1)]"))
	#push!(m.constraints, DisjointUnion(m.varsInter[end], X[1], X[2]))
	Fprec = Variable(m, 1:m.p, "F[$(length(m.varsInter)+1)]")
	DisjointUnion(m, Fprec, X[1], X[2])

	for ind in 3:length(X)
		#push!(m.varsInter, Variable(1:m.p, "F[$(length(m.varsInter)+1)]"))
		#push!(m.constraints, DisjointUnion(m.varsInter[end], m.varsInter[end-1], X[ind]))

		Fnew = Variable(m, 1:m.p, "F[$(length(m.varsInter)+1)]")
		DisjointUnion(m, Fnew, Fprec, X[ind])
		Fprec = Fnew
	end

	Fprec.lowerBound = domain
	Fprec.upperBound = domain
	Fprec.cardinalSup = length(domain)
	Fprec.cardinalInf = length(domain)
	Fprec.isFixed = true

	m.varsInter[end-(length(X)-2):end]

end


text:
function check(p, x, m)
    n = BigInt(0)
    pp = BigInt(1)
    x = reverse(x)
    for d in x
        n += parse(Int, d) * pp
        pp *= p
    end
    # println(" ", p, " ", reverse(x), " ", n)
    n ≤ m
end

function solve()
    x = readline()
    m = parse(BigInt, readline())

    if length(x) == 1
        parse(Int, x) ≤ m ? 1 : 0
    else
        max_ = 10::Int64^18 + 1
        min_min = maximum([parse(Int, a) for a in x]) + 1
        min_ = min_min

        if !check(min_, x, m)
            0
        else
            while min_ + 1 < max_
                # println(min_, " ", max_)
                p = (min_ + max_) ÷ 2
                if check(p, x, m)
                    min_ = p
                else
                    max_ = p
                end
            end
            # println("sol ", min_min, " ", min_, " ", max_)
            min_ - min_min + 1
        end
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
function generate_opinion()
    return 2 * rand() - 1
end

# this function was adapted from:
# https://www.johndcook.com/julia_rng.html
function generate_inclin_interact(lambda=log(25))
    if lambda <= 0.0
        error("mean must be positive")
    end
    -(1 / lambda) * log(rand())
end

function generate_check_regularity()
    return 1 - (rand() / 4)^2
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
"""
Bob is a lackadaisical teenager. In conversation, his responses are very limited.
Bob answers 'Sure.' if you ask him a question, such as "How are you?".
He answers 'Whoa, chill out!' if you YELL AT HIM (in all capitals).
He answers 'Calm down, I know what I'm doing!' if you yell a question at him.
He says 'Fine. Be that way!' if you address him without actually saying anything.
He answers 'Whatever.' to anything else.

Bob's conversational partner is a purist when it comes to written communication
and always follows normal rules regarding sentence punctuation in English.
"""

function bob(stimulus)
  flags = Dict("question" => 0x1b,
               "yell" => 0x10b,
               "empty" => 0x100b,
               "other" => 0x0b)

  response = Dict(0x1b => "Sure.",
                   0x10b => "Whoa, chill out!",
                   0x11b => "Calm down, I know what I'm doing!",
                   0x100b => "Fine. Be that way!",
                   0x0b => "Whatever.")

  seen = 0x0b

#  if contains(stimulus, r"\?$")
  if contains(stimulus, r"\?[-\s]*$")
    seen |= flags["question"]
  end

  if uppercase(stimulus) == stimulus && any(isletter, stimulus)
    seen |= flags["yell"]
  end

  if !contains(stimulus, r"[a-zA-Z0-9]") && !contains(stimulus, r"\?$")
    seen |= flags["empty"]
  end

  return response[seen]
end


text:
toterm(s1::Term, s2::Term) = InteractionTerm((s1, s2))
toterm(s1::InteractionTerm, s2::Term) = InteractionTerm((s1.terms..., s2))
toterm(s1::Term, s2::InteractionTerm) = InteractionTerm((s1, s2.terms...))
toterm(x::Tuple, s2::Term) = map(t -> toterm(t, s2), x)
toterm(s1::Term, x::Tuple) = map(t -> toterm(s1, t), x)
function predictor_product(predictors)
    init, rest = Iterators.peel(predictors)
    reduce((x, y) -> x + y + toterm(x, y), term.(rest), init = term(init))
end
function errorterm(ef, nested)
    terms = term(1) + term(ef)
    if length(nested) > 0
        nestedterms = predictor_product(nested)
        terms += toterm(term(ef), nestedterms)
    end
    terms
end

name(t::AbstractTerm) = "$(t.sym)"
name(t::InteractionTerm) = join(name.(t.terms), " & ")
name(t::InterceptTerm) = "Intercept"
pivoted_asgn(assign, qrp) = assign[qrp.p[1:rank(qrp.R)]]
function mqr(f::FormulaTerm, data; contrasts = Dict{Symbol, Any}())
    mf = ModelFrame(f, data; contrasts)
    mm = ModelMatrix(mf)
    qrp = qr(mm.m, Val(true))
    (qrp = qrp, mf = mf, mm = mm)
end
function aov(f::FormulaTerm, data; return_locals = false)
    qrp, mf, mm = mqr(f, data)
    aov(qrp, pivoted_asgn(mm.assign, qrp), mf.f.rhs.terms,
        getproperty(data, f.lhs.sym); return_locals)
end
aov(l, y) = aov(l..., y)
function aov(qrp, asgn, terms, y; return_locals = false)
    effects = qrp.Q' * y
    r = rank(qrp.R)
    SS = Float64[]
    names = String[]
    DF = Int[]
    for i in union(asgn)
        n = name(terms[i])
        n == "Intercept" && continue
        idxs = findall(asgn .== i)
        push!(SS, sum(abs2, @view(effects[idxs])))
        push!(names, n)
        push!(DF, length(idxs))
    end
    push!(SS, sum(abs2, @view(effects[r+1:end])))
    push!(names, "Residuals")
    push!(DF, length(effects) - r)
    MSS = SS ./ DF
    F = MSS ./ MSS[end]
    F[end] = NaN
    DFres = DF[end]
    result = DataFrame(names = names, DF = DF, SS = SS, MSS = MSS, F = F,
                       p = @.(ccdf(FDist(DF, DFres), F)))
    if return_locals
        result, (qrp = qrp, asgn = asgn, terms = terms)
    else
        result
    end
end
function aoverror(f::FormulaTerm, e, nested::T, data; return_locals = false) where T
    nested = T <: Symbol ? [nested] : nested
    ef = errorterm(e, nested)
    efit = mqr(f.lhs ~ ef, data,
        contrasts = Dict(e => HelmertCoding(),
                         [t => HelmertCoding() for t in nested]...))
    efit.mm.assign = efit.mm.assign[efit.qrp.p[1:rank(efit.qrp.R)]]
    mf = ModelFrame(f, data)
    mm = ModelMatrix(mf)
    ys = efit.qrp.Q' * getproperty(data, f.lhs.sym)
    append!(efit.mm.assign,
            fill(maximum(efit.mm.assign) + 1, length(ys) - length(efit.mm.assign)))
    Xs = efit.qrp.Q' * mm.m
    rows = [findall(efit.mm.assign .== i)
            for i in Iterators.drop(sort!(union(efit.mm.assign)), 1)]
    result = [inner(Xs, ys, row, mm.assign, mf.f.rhs.terms; return_locals)
              for row in rows]
    if return_locals
        first.(result), (Q = efit.qrp.Q,
                         locals_t = [(rest..., row)
                                     for (row, rest) in zip(rows, last.(result))])
    else
        result
    end
end
function aoverror(l, y)
    ys = l.Q' * y
    [aov(qrp, asgn, terms, ys[rows]) for (qrp, asgn, terms, rows) in l.locals_t]
end
function inner(X, y, rows, asgn, terms; return_locals = false)
    T = X[rows, :]
    cols = findall(reshape(sum(abs2, T, dims = 1), :) .> 1e-5)
    qrp = qr(T[:, cols], Val(true))
    aov(qrp, pivoted_asgn(asgn[cols], qrp), terms, y[rows]; return_locals)
end


text:
export TensorSpace,⊗,ProductSpace

#  SV is a tuple of d spaces
abstract AbstractProductSpace{SV,T,d} <: Space{T,AnyDomain,d}

spacetype{SV}(::AbstractProductSpace{SV},k) = SV.parameters[k]

# TensorIterator
# This gives the map from coefficients to the
# tensor entry of a tensor product of d spaces
# findfirst is overriden to get efficient inverse

immutable TensorIterator{DMS<:Tuple}
    dimensions::DMS
end

cache(Q::TensorIterator) = CachedIterator(Q)

Base.eltype{d,T}(::TensorIterator{NTuple{d,T}}) = NTuple{d,Int}
Base.eltype(it::TensorIterator) = NTuple{length(it.dimensions),Int}

Base.start{DMS<:NTuple{2}}(::TensorIterator{DMS}) = (1,1)
Base.start{DMS<:NTuple{3}}(::TensorIterator{DMS}) = (1,1,1)
Base.start(it::TensorIterator) = tuple(ones(Int,d)...)::eltype(it)

function Base.next(it::TensorIterator,st)
    for k=2:length(st)
        if st[k] > 1
            nst=tuple(st[1:k-2]...,st[k-1]+1,st[k]-1,st[k+1:end]...)::eltype(it)
            if all(map(≤,nst,it.dimensions)) || done(it,nst)
                return (st,nst)
            else
                return (st,next(it,nst)[2])
            end
        end
    end
    nst=tuple(ones(Int,length(it.dimensions)-1)...,st[1]+1)::eltype(it)
    if all(map(≤,nst,it.dimensions)) || done(it,nst)
        return (st,nst)
    else
        return (st,next(it,nst)[2])
    end
end

function Base.done(it::TensorIterator,st)
    for k=1:length(st)
        if st[k] ≤ it.dimensions[k]
            return false
        end
    end
    return true
end

function Base.next{d}(it::TensorIterator{NTuple{d,Infinity{Bool}}},st)
    for k=2:length(st)
        if st[k] > 1
            return (st,tuple(st[1:k-2]...,st[k-1]+1,st[k]-1,st[k+1:end]...)::NTuple{d,Int})
        end
    end
    (st,tuple(ones(Int,d-1)...,st[1]+1)::NTuple{d,Int})
end
Base.done{d}(::TensorIterator{NTuple{d,Infinity{Bool}}},st) = false

Base.next(::TensorIterator{NTuple{2,Infinity{Bool}}},st) =
    (st,st[2] == 1? (1,st[1]+1) : (st[1]+1,st[2]-1))

function Base.findfirst(::TensorIterator{NTuple{2,Infinity{Bool}}},kj::Tuple{Int,Int})
    k,j=kj
    if k > 0 && j > 0
        n=k+j-2
        (n*(n+1))÷2+k
    else
        0
    end
end

# which block of the tensor
# equivalent to sum of indices -1

block(it::TensorIterator,k) = sum(it[k])-length(it.dimensions)+1
block(ci::CachedIterator,k) = sum(ci[k])-length(ci.iterator.dimensions)+1

block(::TensorIterator{NTuple{2,Infinity{Bool}}},n) =
    floor(Integer,sqrt(2n) + 1/2)

function getindex(it::TensorIterator{NTuple{2,Infinity{Bool}}},n::Integer)
    m=block(it,n)
    p=findfirst(it,(1,m))
    j=1+n-p
    j,m-j+1
end

blockstarttuple{II}(it::TensorIterator{Tuple{II,Infinity{Bool}}},K) = (1,K)
blockstoptuple{II}(it::TensorIterator{Tuple{II,Infinity{Bool}}},K) = (K,1)

# this is for general dimension case, so we need to construct it so that
# each entry is less than the dimension
function blockstarttuple(it::TensorIterator,K)
    ret=ones(Int,length(it.dimensions))
    for k=reverse(eachindex(ret))
        ret[k]=min(K,it.dimensions[k])
        K-=ret[k]
        if K ≤ 0
            return tuple(ret...)
        end
    end
    tuple(zeros(Int,length(it.dimensions))...)
end

function blockstoptuple(it::TensorIterator,K)
    ret=ones(Int,length(it.dimensions))
    for k=eachindex(ret)
        ret[k]=min(K,it.dimensions[k])
        K-=ret[k]
        if K ≤ 0
            return tuple(ret...)
        end
    end
    tuple(zeros(Int,length(it.dimensions))...)
end

blockstart(it::TensorIterator,K) = findfirst(it,blockstarttuple(it,K))
blockstop(it::TensorIterator,K) = findfirst(it,blockstoptuple(it,K))
blockstart{II,TI<:TensorIterator}(ci::CachedIterator{II,TI},K) =
    findfirst(ci,blockstarttuple(ci.iterator,K))
blockstop{II,TI<:TensorIterator}(ci::CachedIterator{II,TI},K) =
    findfirst(ci,blockstoptuple(ci.iterator,K))

blocklength(it,K::Int) = blockstop(it,K)-blockstart(it,K)+1
blocklength(it,K::Range) = Int[blocklength(it,k) for k in K]
blocklength{II}(it::CachedIterator{II,
            TensorIterator{Tuple{Infinity{Bool},Infinity{Bool}}}},K::Int) = K
blocklength(it::TensorIterator{Tuple{Infinity{Bool},Infinity{Bool}}},K::Int) = K

# convert from block, subblock to tensor
subblock2tensor(rt::TensorIterator{Tuple{Infinity{Bool},Infinity{Bool}}},K,k) =
    (k,K-k+1)

subblock2tensor{II}(rt::CachedIterator{II,TensorIterator{Tuple{Infinity{Bool},Infinity{Bool}}}},K,k) =
    (k,K-k+1)

# TensorSpace
# represents the tensor product of several subspaces

immutable TensorSpace{SV,T,d} <:AbstractProductSpace{SV,T,d}
    spaces::SV
end

tensorizer{SV,T,d}(sp::TensorSpace{SV,T,d}) = TensorIterator(map(dimension,sp.spaces))

TensorSpace(sp::Tuple) =
    TensorSpace{typeof(sp),mapreduce(basistype,promote_type,sp),
                mapreduce(domaindimension,+,sp)}(sp)

dimension(sp::TensorSpace) = mapreduce(dimension,*,sp.spaces)

for OP in (:spacescompatible,:(==))
    @eval $OP{SV,T,d}(A::TensorSpace{SV,T,d},B::TensorSpace{SV,T,d}) =
        all(Bool[$OP(A.spaces[k],B.spaces[k]) for k=1:length(A.spaces)])
end

canonicalspace(T::TensorSpace) = TensorSpace(map(canonicalspace,T.spaces))

coefficient_type(S::TensorSpace,T) =
    mapreduce(sp->coefficient_type(sp,T),promote_type,S.spaces)

TensorSpace(A...) = TensorSpace(tuple(A...))
TensorSpace(A::ProductDomain) = TensorSpace(tuple(map(Space,A.domains)...))
⊗(A::TensorSpace,B::TensorSpace) = TensorSpace(A.spaces...,B.spaces...)
⊗(A::TensorSpace,B::Space) = TensorSpace(A.spaces...,B)
⊗(A::Space,B::TensorSpace) = TensorSpace(A,B.spaces...)
⊗(A::Space,B::Space) = TensorSpace(A,B)

domain(f::TensorSpace) = mapreduce(domain,*,f.spaces)
Space(sp::ProductDomain) = TensorSpace(sp)

*(A::Space,B::Space) = A⊗B

# every column is in the same space for a TensorSpace
#TODO: remove
columnspace(S::TensorSpace,::) = S.spaces[1]

Base.length(d::TensorSpace) = length(d.spaces)
Base.getindex(d::TensorSpace,k::Integer) = d.spaces[k]

immutable ProductSpace{S<:Space,V<:Space,T} <: AbstractProductSpace{Tuple{S,V},T,2}
    spacesx::Vector{S}
    spacey::V
end

ProductSpace(spacesx::Vector,spacey)=ProductSpace{eltype(spacesx),
                                                  typeof(spacey),
                                                  promote_type(basistype(first(spacesx)),basistype(spacey))}(spacesx,spacey)

coefficient_type(S::ProductSpace,T) =
    promote_type(coefficient_type(S.spacesx[1],T),coefficient_type(S.spacesy,T))

⊗{S<:Space}(A::Vector{S},B::Space) = ProductSpace(A,B)
domain(f::ProductSpace) = domain(f.spacesx[1])*domain(f.spacesy)

Base.getindex(d::ProductSpace,k::Integer) = k==1?d.spacesx:d.spacey

space(d::AbstractProductSpace,k) = d[k]
isambiguous(A::TensorSpace) = isambiguous(A[1])||isambiguous(A[2])

Base.transpose(d::TensorSpace) = TensorSpace(d[2],d[1])

##Transforms

plan_column_transform(S,v) = plan_transform(columnspace(S,1),v)
plan_column_itransform(S,v) = plan_itransform(columnspace(S,1),v)

function itransform!(S::TensorSpace,M::Matrix)
    n=size(M,1)

    planc=plan_itransform(space(S,1),M[:,1])
    for k=1:size(M,2)
        M[:,k]=itransform(space(S,1),M[:,k],planc)
    end

    planr=plan_itransform(space(S,2),vec(M[1,:]))
    for k=1:n
        M[k,:]=itransform(space(S,2),vec(M[k,:]),planr)
    end
    M
end

function itransform!(S::AbstractProductSpace,M::Matrix)
    n=size(M,1)

    ## The order matters
    pln=plan_column_itransform(S,n)
    for k=1:size(M,2)
        M[:,k]=itransform(columnspace(S,k),M[:,k],pln)
    end

    for k=1:n
        M[k,:]=itransform(space(S,2),vec(M[k,:]))
    end
    M
end

function transform!(S::TensorSpace,M::Matrix)
    n=size(M,1)

    planc=plan_transform(space(S,1),M[:,1])
    for k=1:size(M,2)
        M[:,k]=transform(space(S,1),M[:,k],planc)
    end

    planr=plan_transform(space(S,2),vec(M[1,:]))
    for k=1:n
        M[k,:]=transform(space(S,2),vec(M[k,:]),planr)
    end
    M
end

function transform!{T}(S::AbstractProductSpace,M::Matrix{T})
    n=size(M,1)

    ## The order matters!!
    # For Disk Space, this is due to requiring decay
    # in function
    for k=1:n
        M[k,:]=transform(space(S,2),vec(M[k,:]))
    end

    pln=plan_column_transform(S,n)
    for k=1:size(M,2)
        # col may not be full length
        col=transform(columnspace(S,k),M[:,k],pln)
        M[1:length(col),k]=col
        for j=length(col)+1:n
            M[j,k]=zero(T) # fill rest with zeros
        end
    end

    M
end

## points

points(d::Union{BivariateDomain,BivariateSpace},n,m) =
    points(d,n,m,1),points(d,n,m,2)

function points(d::BivariateSpace,n,m,k)
    ptsx=points(columnspace(d,1),n)
    ptst=points(space(d,2),m)

    promote_type(eltype(ptsx),eltype(ptst))[fromcanonical(d,x,t)[k] for x in ptsx, t in ptst]
end

##  Fun routines

function fromtensor{T}(it::TensorIterator{NTuple{2,Infinity{Bool}}},M::Matrix{T})
    ret=zeros(T,findfirst(it,(size(M,1),size(M,2))))
    for k=1:size(M,1),j=1:size(M,2)
        ret[findfirst(it,(k,j))] = M[k,j]
    end
    ret
end

# function totensor{T1,T2,T}(it::TensorIterator{Tuple{T1,T2}},M::Vector{T})
#     N=length(M)
#     inds=it[N]
#     m=inds[1]+inds[2]-1
#     ret=zeros(T,min(it.dimensions[1],m),min(it.dimensions[2],m))
#     k=1
#     for c in it
#         ret[c...] = M[k]
#         k+=1
#         if k > N
#             break
#         end
#     end
#     ret
# end

function totensor{T}(it::TensorIterator{NTuple{2,Infinity{Bool}}},M::Vector{T})
    inds=it[length(M)]
    m=inds[1]+inds[2]-1
    ret=zeros(T,m,m)
    for k=1:length(M)
        ret[it[k]...] = M[k]
    end
    ret
end

for OP in (:fromtensor,:totensor,:block,:blockstart,:blockstop)
    @eval $OP(s::Space,M) = $OP(tensorizer(s),M)
end

# TODO: remove
function totree(v::Vector)
   m=toblock(length(v))
    r=Array(Vector{eltype(v)},m)
    for k=1:m-1
        r[k]=v[fromblock(k)]
    end
    r[m]=pad!(v[fromblock(m)[1]:end],m)
    r
end

fromtree{T}(v::Vector{Vector{T}}) = vcat(v...)

function points(sp::TensorSpace,n)
    pts=Array(Tuple{Float64,Float64},0)
    for x in points(sp[1],round(Int,sqrt(n))), y in points(sp[2],round(Int,sqrt(n)))
        push!(pts,(x,y))
    end
    pts
end

function transform(sp::TensorSpace,vals)
    m=round(Int,sqrt(length(vals)))
    M=reshape(copy(vals),m,m)

    fromtensor(sp,transform!(sp,M))
end

evaluate(f::AbstractVector,S::AbstractProductSpace,x) = ProductFun(totensor(S,f),S)(x...)
evaluate(f::AbstractVector,S::AbstractProductSpace,x,y) = ProductFun(totensor(S,f),S)(x,y)

coefficientmatrix{S<:AbstractProductSpace}(f::Fun{S}) = totensor(space(f),f.coefficients)

Fun{T<:Number}(v::Vector{Vector{T}},S::TensorSpace) = Fun(fromtree(v),S)

#TODO: Implement
# function ∂(d::TensorSpace{Interval{Float64}})
#     @assert length(d.spaces) ==2
#     PiecewiseSpace([d[1].a+im*d[2],d[1].b+im*d[2],d[1]+im*d[2].a,d[1]+im*d[2].b])
# end

union_rule(a::TensorSpace,b::TensorSpace) = TensorSpace(map(union,a.spaces,b.spaces))


text:
import .QueueUnion: Range, Queues
import Printf: @printf

const EPS = 1e-10
const DEBUG = false

@deprecate(
    clip_front(elements, pqs, i, slope, offset, t),
    clip(elements, Ref(pqs, i), +slope, +offset - t, Val(true))
)

@deprecate(
    clip_back(elements, pqs, i, slope, offset, t),
    clip(elements, Ref(pqs, i), -slope, -offset + t, Val(false))
)

@deprecate(
    clip_front(qs::Queues{Event}, i::I, slope::F, offset::F, t::F) where {F,I},
    clip(qs, i, slope, offset - t, Val(true))
)

@deprecate(
    clip_back(qs::Queues{Event}, i::I, slope::F, offset::F, t::F) where {F,I},
    clip(qs, i, -slope, -offset + t, Val(false))
)

clip(qs::Queues{Event}, i::I, slope::F, offset::F, v::Val{D}) where {F,I,D} =
    clip(qs.events, Ref(qs.pq, i), slope, offset, v)

function clip(
    elements::Vector{Event},
    pq::Ref{Range},
    slope::F,
    offset::F,
    ::Val{forward},
    ::Val{check},
)::F where {F,forward,check}
    local start::Int = pq[].start
    local stop::Int = pq[].stop
    local e::Event = elements[forward ? start : stop]::Event
    @static DEBUG && local dir::String = forward ? "f" : "b"
    @static DEBUG && @printf("clip_%s: (%+g, %+.2f) [%d] %s\n",
                             dir, slope, offset, check, string(elements[start:stop]))
    while start <= stop && slope * e.x + offset < 0
        offset += intercept(e)
        slope += e.slope
        @static DEBUG && @printf(" lip_%s: (%+g, %+.2f)\n", dir, slope, offset)
        e = elements[forward ? start += 1 : stop -= 1]
    end
    local x::F = if check && abs(slope) <= EPS
        forward ? -Inf : +Inf
    else
        let x::F = -offset/slope
            elements[forward ? start -= 1 : stop += 1] = Event(x, slope)
            x
        end
    end
    @static DEBUG && @printf("  ip_%s:  --> %g\n", dir, x)
    pq[] = Range(start, stop)
    return x
end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
function analyze()

  measure, measureinfo = loadmeasure()
  @unpack prefix, analysispath, outputpath, resultid = measureinfo

  try
    PARAM[:refreshverify] && verifymeasure(; measure, measureinfo)
    PARAM[:refreshleadlag] && analyzeleadlag(; measure, measureinfo)

  catch err

    mv(outputpath, "$(analysispath)\\$prefix\\ERR_$resultid")
    #show(stdout, MIME"text/plain"(), stacktrace(catch_backtrace()))
    @error "Verificaiton failed" exception=(err, catch_backtrace())
    throw("Verification failed")
  end

end

#performs file io and initial setup for a new measure estimation
function loadmeasure(;
    burnin::DatePeriod = PARAM[:measureburnin],
    measurefilename::String = PARAM[:analysismeasurefilename],
    analysispath = PARAM[:analysispath],
    returntype::String = PARAM[:crspdatalabel],
    workingpath::String = PARAM[:workingpath],
    hfrname::String = PARAM[:hfrname],
    mfname::String = PARAM[:mfname],
    plname::String = PARAM[:plname],)

  datefrequency = Dict("w"=>:week, "m"=>:month)[PARAM[:crspfrequencysuffix]]

  prefix = (returntype === measurefilename[1:length(returntype)]
    ) ? measurefilename : "$(returntype)_$measurefilename"

  #create a directory to hold anlsyis results
  resultid = "rid$(Dates.format(now(),"yymmdd-HHMM"))"

  #load the measure file
  local measure::DataFrame
  inputpath = "$analysispath\\$prefix\\inputs"

  #each results set gets its own output path- do the below only once
  if !ispath("$analysispath\\$prefix")
    mkdir("$analysispath\\$prefix")
    mkdir("$analysispath\\$prefix\\inputs")
    mv("$analysispath\\$measurefilename.csv",  "$inputpath\\$measurefilename.csv")
    mv("$analysispath\\$measurefilename.txt",  "$inputpath\\$measurefilename.txt")
    cp("$workingpath\\$(returntype)_$(hfrname).csv", "$inputpath\\$(returntype)_$(hfrname).csv")
    cp("$workingpath\\$(returntype)_$(mfname).csv", "$inputpath\\$(returntype)_$(mfname).csv")
    #NOTE: momentum info is disabled for now- need to re-run data if going to try it
    if isfile("$workingpath\\$(returntype)_$(plname).csv")
      @warn "Momentum file detected! Note you will need to disable the error barrier in order
        to actually run the momentum tests"
      cp("$workingpath\\$(returntype)_$(plname).csv", "$inputpath\\$(returntype)_$(plname).csv")
    end
  else
    isfile("$analysispath\\$measurefilename.csv") && throw("Duplicate measure file
      $measurefilename found in analysis folder outside of input/results folder $prefix.
      This should not happen- measurefiles should correspond 1:1 with the input/results folder,
      unless the input/results folder does not exist (at which point it will be created
      immediately).")
  end

  #load the input data
  measure = CSV.File("$inputpath\\$measurefilename.csv") |> DataFrame

  #create a directory to hold anlsyis results
  resultid = "rid$(Dates.format(now(),"yymmdd-HHMM"))"
  outputpath = "$analysispath\\$prefix\\$resultid"

  @assert !ispath("$analysispath\\$prefix\\$resultid")
  mkdir(outputpath)
  cp(PARAMETER_FILE, "$outputpath\\$(resultid)_Parameters.xlsx")

  #now pre-process the columns
  #sums accross assets
  if datefrequency === :week && PARAM[:analysisaggregateweeklycols]
    assetmat = measure[!, findcols(measure, :asset)] |> Matrix{MFloat64}
    oldcols = [findcols(measure, :growth); findcols(measure, :asset);
      propertynames(measure[!, r"^V_"])]
    aggassets = sum(assetmat, dims=2) |> vec
    measure[!, Symbol(:A_, weeklyaggregatename)] = aggassets
    measure[!, Symbol(:G_, weeklyaggregatename)] = [
      missing; aggassets[2:end] ./ aggassets[1:(end-1)]; ]

    #also need to some accross volume
    volumemat = measure[!, propertynames(measure[!, r"^V_"])] |> Matrix{MFloat64}
    aggvolume = sum(volumemat, dims=2) |> vec
    measure[!, Symbol(:V_, weeklyaggregatename)] = aggvolume
    select!(measure, Not(oldcols))
  end

  #pre-process the measure data
  measure = measure[measure.date .≥ (minimum(measure.date) + burnin), :]
  derivemeasurecols!(measure; datefrequency)

  #setup neweywest SEs and moving block bootstrap
  #below was previously ceil(max(nrow(measure),nrow(funds))^0.25) |> Int
  neweylags::Int = nrow(measure)^0.25 |> ceil |> Int
  nwΣ = neweywestΣfunc(neweylags)

  bsaggfunc(v::AbstractVector) = (;se=std(skipmissing(v)), p=nonparametricp(v))
  #bsaggfunc(v) = (;se=std(skipmissing(v)), p=nonparametricp(v))

  bsmeasurecolsregex::Regex = Regex(
    PARAM[:bsmeasurecolsregex] === nothing ? "" : PARAM[:bsmeasurecolsregex])
  bsbenchmarkcolsregex::Regex = Regex(
    PARAM[:bsbenchmarkcolsregex] === nothing ? "" : PARAM[:bsbenchmarkcolsregex])
  bssamples = PARAM[:bssamples]

  #note- the below was previously #ceil(max(nrow(measure),nrow(funds))^(1/3)) |> Int
  bswidth = nrow(measure)^(1/3) |> ceil |> Int

  #we can reuse some of the bootstrap aggregation code for the null simulations
  nullar1aggfunc(v::AbstractVector,statname,stattestvals) = (;
    se=std(skipmissing(v)),
    p=nonparametricnullp(v, stattestvals[statname]))
  nullar1aggfunc(M::AbstractMatrix,statname,stattestvals) = (;
    se=std(M, dims=1) |> vec,
    p=nonparametricnullp(M, stattestvals[statname]))
  nullar1simulations = bssamples

  #lightly package the io data for use as needed
  measureinfo = (; resultid, inputpath, outputpath, returntype,
    analysispath, datefrequency, prefix,
    neweylags, nwΣ,
    bsaggfunc, bsmeasurecolsregex, bsbenchmarkcolsregex, bssamples, bswidth,
    nullar1aggfunc, nullar1simulations)

  return (measure, measureinfo)
end

#helpful for grabbing the column name sans prefix
stripcolprefix(s) = replace(s, r"^l?[A-Z0-9]+_"=>"") |> Symbol
stripcolprefix(s::Symbol) = stripcolprefix(string(s))

#compute the strategy pressure measure- AUM×turnover/aggvol
function derivemeasurecols!(measure;
    Fassets::Vector{Symbol}=findcols(measure, :asset, includelp=false),
    datefrequency::Symbol,
    lpfilterperiods::Vector{<:Union{Day, Nothing}} = PARAM[:analysislpfilterperiods],
    Ftotalvolume::Symbol = Symbol(:desc_, PARAM[:analysistotalvolume]),
    lladditionalreturncolsidx = PARAM[:lladditionalreturncolsidx],
    lladditionalnormalizers = PARAM[:lladditionalnormalizers],
    llmarketcol = PARAM[:llmarketcol])

    #other helper methods
    simpledif(v::AbstractVector) = [missing; v[2:end] .- v[1:(end-1)]]
    _logdif(v) = [missing; log.(v[2:end]) .- log.(v[1:(end-1)])]
    logdif(v) = _logdif(abs.(v)) .* sign.(v)

    #generate log growth columns
    for gcol ∈ findcols(measure, :growth)
      measure[!, gcol] = [missing; measure[2:end, gcol]]
      measure[!, Symbol(:l,gcol)] = log.(abs.(measure[!, gcol])) .* sign.(measure[!, gcol])
      measure[!, Symbol(:g,string(gcol)[2:end])] = measure[!, gcol] .- 1

      rcol = Symbol(:R_, stripcolprefix(gcol))
      acol = Symbol(:A_, stripcolprefix(gcol))
      agcol = Symbol(:ag, string(gcol)[2:end])
      ngcol = Symbol(:ng, string(gcol)[2:end])
      cggcol = Symbol(:cgg, string(gcol)[2:end])
      measure[!, agcol] = measure[!, gcol] .- measure[!, rcol] .- 1.0
      measure[!, ngcol] = measure[!, gcol] .- measure[!, llmarketcol] .- 1.0
      measure[!, cggcol] = measure[!, rcol] |> deepcopy

      #active asset growth
      measure[!, Symbol(:aA, string(gcol)[2:end])] = [missing; measure[1:(end-1), acol] .* measure[2:end, agcol]]
      measure[!, Symbol(:nA, string(gcol)[2:end])] = [missing; measure[1:(end-1), acol] .* measure[2:end, ngcol]]
    end

    #Fassetcols = assetcols(measure)
    #make growth cols
    measure = lpfilter(measure,focalcols = [
      findcols(measure, :loggrowth); findcols(measure, :asset); :mctotal])
    colstodif = findcols(measure, :loggrowth)
    transform!(measure,  colstodif .=>
      (v -> [missing; v[2:end] .- v[1:(end-1)]])
      .=> (s->Symbol(:G, s)).(colstodif))

    Fvols = propertynames(measure[!, r"^V_"])
    Fpurchases = propertynames(measure[!, r"^Vs_"])
    Fpressures = (s->Symbol(:P_, stripcolprefix(s))).(Fassets)
    Fparticipations = (s->Symbol(:p_, stripcolprefix(s))).(Fassets)
    Flgpressures = (s->Symbol(:lG, s)).(Fpressures)
    totalvolume = measure[!, Ftotalvolume]
    totalcapitalization = measure[!, :mctotal]

    Fpurchaseparticipations = (s->Symbol(:pp_, stripcolprefix(s))).(Fassets)
    Fpurchasefractions = (s->Symbol(:VsA_, stripcolprefix(s))).(Fassets)
    Fpurchasecaps = (s->Symbol(:VsMC_, stripcolprefix(s))).(Fassets)

    #compute the presure characteristic
    transform!(measure, Fvols .=>
      ((vol)-> vol ./ totalvolume) .=> Fparticipations)
    transform!(measure, Fpurchases .=>
      ((pur)-> pur ./ totalvolume) .=> Fpurchaseparticipations)
    transform!(measure, zip(Fassets, Fpurchases) .|> collect .=>
      ((asset,pur)-> pur ./ asset) .=> Fpurchasefractions)
    transform!(measure,  Fpurchases .=>
      ((pur)-> pur ./ totalcapitalization) .=> Fpurchasecaps)
    transform!(measure, Fvols .=>
      ((vol)-> log.(vol)) .=> ((s->Symbol(:l,s)).(Fvols)))
    transform!(measure, zip(Fassets, Fvols) .|> collect .=>
      ((asset, vol)->asset .* vol ./ totalvolume) .=> Fpressures)
    transform!(measure, Fpressures .=> logdif .=> Flgpressures)
    lpfilter(measure, focalcols=Flgpressures)
    Flgpressures = [Flgpressures; lpcols.(Flgpressures)...;]

    #compute a version just of the volume-pressure
    Fvpressures = (s->Symbol(:VP_, stripcolprefix(s))).(Fassets)
    Flgvpressures = (s->Symbol(:lGVP_, stripcolprefix(s))).(Fassets)
    transform!(measure, Fvols .=> (v -> logdif(v ./ totalvolume)) .=> Flgvpressures)
    lpfilter(measure, focalcols=Flgvpressures)
    Flgvpressures = [Flgvpressures; lpcols.(Flgvpressures)...]

    #compute a measure of dollar return
    Frets = (s->Symbol(:R_,s)).(stripcolprefix.(Fassets))
    Fdolrets = (s->Symbol(:X,s)).(Frets)
    transform!(measure, zip(Fassets, Frets) .|> collect .=> (
      (a,r)-> a .* r) .=> Fdolrets)

    #somewhat different calculation methedology here for the non-logged version
    #in order to avoid various operations on ratios
    lpfilter(measure, focalcols=Fvols;) #[Fvols; Ftotalvolume;])
    Fvpressures = [Fvpressures; lpcols.(Fvpressures)...;]
    Fgvpressures = (s->Symbol(:G, s)).(Fvpressures)
    Fvolsforvpressures = [Fvols; lpcols.(Fvols)...;]
    Ftotalvolumes = [Ftotalvolume for i ∈ 1:length(Fvolsforvpressures)]
    transform!(measure, zip(Fvolsforvpressures, Ftotalvolumes) .|> collect .=>
      ((v,tv) -> simpledif(v) ./ tv) .=> Fgvpressures)

    #generate the appropriate column labels
    Fmcassetcols = [Fassets; lpcols.(Fassets)...;] #all asset cols
    Fmccols = findcols(measure, :mc)
    Fmccols = [Fmccols; (Fmccol->[Fmccol for i in 1:length(lpcols(Fmccol))]).(Fmccols)...;]
    Famccols = Fassets .|>  stripcolprefix .|> (s)->Symbol(:M_, s)
    Famccols = [Famccols; lpcols.(Famccols)...;]
    Fgamccols = Famccols .|> s-> Symbol(:G, s)

    #these will be used for return predictability
        #WARNING: quick and dirty test
    #measure.mctotal = [missing; measure.mctotal[1:(end-1)]]
    lpfilter(measure, focalcols=[:mctotal]) #need to smooth the mctotal
    Fmctotcols = [:mctotal for i ∈ 1:length(Fassets)]
    Fmctotcols = [Fmctotcols; (Fmctotcol->
      [Fmctotcol for i in 1:length(lpcols(Fmctotcol))]).(Fmctotcols)...;]
    Famctotcols = Fassets .|>  stripcolprefix .|> (s)->Symbol(:MT_, s)
    Famctotcols = [Famctotcols; lpcols.(Famctotcols)...;]
    Flamctotcols = Famctotcols .|> s-> Symbol(:l, s)
    Fgamctotcols = Famctotcols .|> s-> Symbol(:G, s)

    Faamcactiveassetcols = Fassets .|> (s)->Symbol(:a, s)
    Faamcactiveassetcols = [Faamcactiveassetcols; lpcols.(Faamcactiveassetcols)...;]
    Fagamctotcols = Famctotcols .|> s-> Symbol(:aG, s)

    Fnamcactiveassetcols = Fassets .|> (s)->Symbol(:n, s)
    Fnamcactiveassetcols = [Fnamcactiveassetcols; lpcols.(Fnamcactiveassetcols)...;]
    Fngamctotcols = Famctotcols .|> s-> Symbol(:nG, s)

    #NOTE we take the ratio of the average rather than the average of the ratio,
    #while also following the Welch 2020 RoC approach

    transform!(measure,
      zip(Fmcassetcols, Fmccols) .|> collect .=>
        ((asset, mc)->(asset ./ mc)) .=> Famccols,
      zip(Fmcassetcols, Fmccols) .|> collect .=>
        ((asset, mc)->(simpledif(asset) ./ mc)) .=> Fgamccols,

      zip(Fmcassetcols, Fmctotcols) .|> collect .=>
        ((asset, mctot)->(asset ./ mctot)) .=> Famctotcols,
      zip(Fmcassetcols, Fmctotcols) .|> collect .=>
        ((asset, mctot)->log.(abs.(asset ./ mctot)) .* sign.(asset)) .=> Flamctotcols,
      zip(Fmcassetcols, Fmctotcols) .|> collect .=>
        ((asset, mctot)->(simpledif(asset) ./ mctot)) .=> Fgamctotcols,

      zip(Faamcactiveassetcols, Fmctotcols) .|> collect .=>
        ((activeasset, mctot)->(activeasset ./ mctot)) .=> Fagamctotcols,
      zip(Fnamcactiveassetcols, Fmctotcols) .|> collect .=>
        ((netasset, mctot)->(netasset ./ mctot)) .=> Fngamctotcols,
        )

    #the below algorithm normalizes for an arbitrary denominator
    for (prefix, Fdenom) ∈ zip(keys(lladditionalnormalizers), values(lladditionalnormalizers))
      @assert Fdenom ∈ propertynames(measure)
      lpfilter(measure, focalcols=[Fdenom]) #need to smooth the denom

      #form the denominator names
      Fdenoms = [Fdenom for i ∈ 1:length(Fassets)]
      Fdenoms = [Fdenoms; (Fdenomcol->
        [Fdenomcol for i in 1:length(lpcols(Fdenomcol))]).(Fdenoms)...;]

      #form the names for the new fields
      Fadenoms = Fassets .|>  stripcolprefix .|> (s)->Symbol(prefix, s)
      Fadenoms = [Fadenoms; lpcols.(Fadenoms)...;]
      Fladenoms = Fadenoms .|> s-> Symbol(:l, s)
      Fgadenoms = Fadenoms .|> s-> Symbol(:G, s)
      Fagadenoms = Fadenoms .|> s-> Symbol(:aG, s)
      Faladenoms = Fadenoms .|> s-> Symbol(:al, s)

      #normalize as necesary
      transform!(measure,
        zip(Fmcassetcols, Fdenoms) .|> collect .=>
          ((asset, denom)->(asset ./ denom)) .=> Fadenoms,
        zip(Fmcassetcols, Fdenoms) .|> collect .=>
          ((asset, denom)->log.(abs.(asset ./ denom)) .* sign.(asset)) .=> Fladenoms,
        zip(Fmcassetcols, Fdenoms) .|> collect .=>
          ((asset, denom)->(simpledif(asset) ./ denom)) .=> Fgadenoms,
        zip(Faamcactiveassetcols, Fdenoms) .|> collect .=>
          ((activeasset, denom)->(activeasset ./ denom)) .=> Fagadenoms,)

    end

    #lgvpressures doesn't include lpcols, but Fgvpressures does
    allcols = [Flgpressures; Fgvpressures; Flgvpressures]
    doubledifcols = allcols .|> s->Symbol(:G, s)
    transform!(measure, allcols .=> simpledif .=> doubledifcols)

    measure |> CSV.write(
      "$(PARAM[:testpath])\\$(PARAM[:analysismeasurefilename])_derive.csv")

    return measure
    #Fpressures =
end

#helper functions for acquiring data from the measure file
function findcols(df, coltype::Symbol; includelp=true, args...)::Vector{Symbol}

  simplepatterns = Dict(
    :loggrowth => r"^lG_",
    :growth => r"^G_",
    :growthM1 => r"^g_",
    :activegrowth => r"^ag_",
    :netgrowth => r"^ng_",
    :capitalgainsgrowth => r"^cgg_",
    :activeassetalloc => r"^aA_",
    :z => r"^Z_",
    :volume => r"^V_",
    :lvolume => r"^lV_",
    :purchases => r"^Vs_",
    :purchasefraction => r"^VsA_",
    :purchasecapitalization => r"^VsMC_",
    :return => r"^R_",
    :ff => r"^F_",
    :logpressure => r"^lP_",
    :growthpressure => r"^GP_",
    :loggrowthpressure => r"^lGP_",
    :participation => r"^p_",
    :purchaseparticipation => r"^pp_",
    :sumloggrowthpressure => r"^lGSP_",
    :growthvolumepressure => r"^GVP_",
    :loggrowthvolumepressure => r"^lGVP_",
    :assetmarketcap => r"^M_",
    :assetmarketcapgrowth => r"^GM_",

    :assetmf => r"^AN_",
    :lassetmf => r"^lAN_",
    :assetmfgrowth => r"^GAN_",
    :activelassetmf => r"^alAN_",
    :activeassetmfgrowth => r"^aGAN_",

    :assetmfequity => r"^EN_",
    :lassetmfequity => r"^lEN_",
    :assetmfequitygrowth => r"^GEN_",
    :activeassetmfequitygrowth => r"^aGEN_",

    :assetmfequityonly => r"^EON_",
    :lassetmfequityonly => r"^lEON_",
    :assetmfequityonlygrowth => r"^GEON_",
    :activeassetmfequityonlygrowth => r"^aGEON_",

    :assetmfdomesticeo => r"^DEON_",
    :lassetmfdomesticeo => r"^lDEON_",
    :assetmfdomesticeogrowth => r"^GDEON_",
    :activeassetmfdomesticeogrowth => r"^aGDEON_",

    :assetmarketcaptotal => r"^MT_",
    :lassetmarketcaptotal => r"^lMT_",
    :assetmarketcaptotalgrowth => r"^GMT_",
    :activeassetmarketcaptotalgrowth => r"^aGMT_",
    :netassetmarketcaptotalgrowth => r"^nGMT_",
    )

  local foundcols
  if coltype ∈ keys(simplepatterns)
    foundcols = propertynames(df[!, simplepatterns[coltype]])
  else
    foundcols = findcols(df, Val{coltype}(); args...)
  end

  if !includelp #sometimes we may want to exclude the filtered cols
    foundcols = setdiff(foundcols, propertynames(df[!, r"_LP[0-9]+d$"]))
  end
  return foundcols
end

findcols(measure::AbstractDataFrame, ::Val{:asset}; includelp=true) = setdiff(
  propertynames(measure[!, r"^A_"]), includelp ? [] : propertynames(measure[!, r"_LP[0-9]+d$"]))
findcols(measure::AbstractDataFrame, ::Val{:mc};
  mctype=:gross) = propertynames(measure[!, Regex("^mc$(mctype)_")])

#generic fallback
findcols(measure::AbstractDataFrame, reg::Regex; includelp=true) = setdiff(
  propertynames(measure[!, reg]), includelp ? [] : propertynames(measure[!, r"_LP[0-9]+d$"]))
findcols(measure::AbstractDataFrame, reg::String; args...) = findcols(
  measure, Regex(reg), args...)

lpcols(F::Symbol;
  lpfilterperiods::Vector{<:Union{Day, Nothing}} = PARAM[:analysislpfilterperiods]) = (
    (p->Symbol(F, lpsuffix(p))).(setdiff(lpfilterperiods, [nothing])))

differencecol(s::Symbol) = Symbol(:G, s)


text:
"""
    warp(img, tform, [indices]; kwargs...) -> imgw

Transform the coordinates of `img`, returning a new `imgw` satisfying `imgw[I] = img[tform(I)]`.

# Output

The output array `imgw` is an `OffsetArray`. Unless manually specified, `axes(imgw) == axes(img)`
does not hold in general. If you just want a plain array, you can "strip" the custom indices with
`parent(imgw)` or `OffsetArrays.no_offset_view(imgw)`.

# Arguments

- `img`: the original image that you need coordinate transformation.
- `tform`: the coordinate transformation function or function-like object, it must accept a
  [`SVector`](https://github.com/JuliaArrays/StaticArrays.jl) as input. A useful package to
  create a wide variety of such transfomrations is
  [CoordinateTransformations.jl](https://github.com/FugroRoames/CoordinateTransformations.jl).
- `indices` (Optional): specifies the output image axes.
  By default, the indices are computed in such a way that `imgw` contains all the original pixels
  in `img` using [`autorange`](@ref ImageTransformations.autorange). To do this `inv(tform)` has
  to be computed. If the given transfomration `tform` does not support `inv` then the parameter
  `indices` has to be specified manually.

# Parameters

!!! info
    To construct `method` and `fillvalue` values, you may need to load `Interpolations` package first.

- `method::Union{Degree, InterpolationType}`: the interpolation method you want to use. By default it is
  `BSpline(Linear())`. To construct the method instance, one may need to load `Interpolations`.
- `fillvalue`: the value that used to fill the new region. The default value is `NaN` if possible,
  otherwise is `0`. One can also pass the extrapolation boundary condition: `Flat()`, `Reflect()` and `Periodic()`.

# See also

There're some high-level interfaces of `warp`:

- image rotation: [`imrotate`](@ref)
- image resize: [`imresize`](@ref)

There are also lazy versions of `warp`:

- [`WarpedView`](@ref) is almost equivalent to `warp` except that it does not allocate memory.
- [`InvWarpedView(img, tform, [indices]; kwargs...)`](@ref ImageTransformations.InvWarpedView)
  is almost equivalent to `warp(img, inv(tform), [indices]; kwargs...)` except that it does not
  allocate memory.

# Extended help

## Parameters in detail

This approach is known as backward mode warping. It is called "backward" because
the internal coordinate transformation is actually an inverse map from `axes(imgr)` to `axes(img)`.

You can manually specify interpolation behavior by constructing `AbstractExtrapolation` object
and passing it to `warp` as `img`. However, this is usually cumbersome. For this reason, there
are two keywords `method` and `fillvalue` to conveniently construct an `AbstractExtrapolation`
object during `warp`.

!!! warning
    If `img` is an `AbstractExtrapolation`, then additional `method` and `fillvalue` keywords
    will be discarded.

### `method::Union{Degree, InterpolationType}`

The interpolation method you want to use to reconstruct values in the wrapped image.

Among those possible `InterpolationType` choice, there are some commonly used methods that you may
have used in other languages:

- nearest neighbor: `BSpline(Constant())`
- triangle/bilinear: `BSpline(Linear())`
- bicubic: `BSpline(Cubic(Line(OnGrid())))`
- lanczos2: `Lanczos(2)`
- lanczos3: `Lanczos(3)`
- lanczos4: `Lanczos(4)` or `Lanczos4OpenCV()`

When passing a `Degree`, it is expected to be a `BSpline`. For example, `Linear()` is equivalent to
`BSpline(Linear())`.

### `fillvalue`

In case `tform(I)` maps to indices outside the original `img`, those locations are set to a value
`fillvalue`. The default fillvalue is `NaN` if the element type of `img` supports it, and `0`
otherwise.

The parameter `fillvalue` can be either a `Number` or `Colorant`. In this case, it will be
converted to `eltype(imgr)` first. For example, `fillvalue = 1` will be converted to `Gray(1)` which
will fill the outside indices with white pixels.

Also, `fillvalue` can be extrapolation schemes: `Flat()`, `Periodic()` and `Reflect()`. The best
way to understand these schemes is perhaps try it with small example:

```jldoctest
using ImageTransformations, TestImages, Interpolations
using OffsetArrays: IdOffsetRange

img = testimage("lighthouse")

imgr = imrotate(img, π/4; fillvalue=Flat()) # zero extrapolation slope
imgr = imrotate(img, π/4; fillvalue=Periodic()) # periodic boundary
imgr = imrotate(img, π/4; fillvalue=Reflect()) # mirror boundary

axes(imgr)

# output

(IdOffsetRange(values=-196:709, indices=-196:709), IdOffsetRange(values=-68:837, indices=-68:837))
```

## The meaning of the coordinates

`imgw` keeps track of the indices that would result from applying `inv(tform)` to the indices of
`img`. This can be very handy for keeping track of how pixels in `imgw` line up with
pixels in `img`.

```jldoctest
using ImageTransformations, TestImages, Interpolations

img = testimage("lighthouse")
imgr = imrotate(img, π/4)
imgr_cropped = imrotate(img, π/4, axes(img))

# No need to manually calculate the offsets
imgr[axes(img)...] == imgr_cropped

# output
true
```

!!! tip
    For performance consideration, it's recommended to pass the `inds` positional argument to
    `warp` instead of cropping the output with `imgw[inds...]`.

# Examples: a 2d rotation

!!! note
    This example only shows how to construct `tform` and calls `warp`. For common usage, it is
    recommended to use [`imrotate`](@ref) function directly.

Rotate around the center of `img`:

```jldoctest
using ImageTransformations, CoordinateTransformations, Rotations, TestImages, OffsetArrays
using OffsetArrays: IdOffsetRange
img = testimage("lighthouse") # axes (1:512, 1:768)

tfm = recenter(RotMatrix(-pi/4), center(img))
imgw = warp(img, tfm)

axes(imgw)

# output

(IdOffsetRange(values=-196:709, indices=-196:709), IdOffsetRange(values=-68:837, indices=-68:837))
```

"""
function warp(img::AbstractExtrapolation{T}, tform, inds::Tuple = autorange(img, inv(tform))) where T
    out = similar(Array{T}, inds)
    warp!(out, img, try_static(tform, img))
end

function warp!(out, img::AbstractExtrapolation, tform)
    @inbounds for I in CartesianIndices(axes(out))
        # Backward mode:
        #   1. get the target index `I` of `out`
        #   2. maps _back_ to original index `Ĩ` of `img`
        #   3. interpolate/extrapolate the value of `Ĩ`
        #   4. this value is then assigned to `out[I]`
        # The advantage of backward mode is that all piexels
        # in the output image will be iterated once very efficiently.
        out[I] = _getindex(img, tform(SVector(I.I)))
    end
    out
end

function warp(img::AbstractArray, tform, args...; kwargs...)
    etp = box_extrapolation(img; kwargs...)
    warp(etp, try_static(tform, img), args...)
end

"""
    imrotate(img, θ, [indices]; kwargs...) -> imgr

Rotate image `img` by `θ`∈[0,2π) in a clockwise direction around its center point.

# Arguments

- `img::AbstractArray`: the original image that you need to rotate.
- `θ::Real`: the rotation angle in clockwise direction.
  To rotate the image in conter-clockwise direction, use a negative value instead.
  To rotate the image by `d` degree, use the formular `θ=d*π/180`.
- `indices` (Optional): specifies the output image axes. By default, rotated image `imgr` will not be
  cropped, and thus `axes(imgr) == axes(img)` does not hold in general.

# Parameters

!!! info
    To construct `method` and `fillvalue` values, you may need to load `Interpolations` package first.

- `method::Union{Degree, InterpolationType}`: the interpolation method you want to use. By default it is
  `Linear()`.
- `fillvalue`: the value that used to fill the new region. The default value is `NaN` if possible,
  otherwise is `0`.

This function is a simple high-level interface to `warp`, for more explaination and details,
please refer to [`warp`](@ref).

# Examples

```julia
using TestImages, ImageTransformations
img = testimage("cameraman")

# Rotate the image by π/4 in the clockwise direction
imgr = imrotate(img, π/4) # output axes (-105:618, -105:618)

# Rotate the image by π/4 in the counter-clockwise direction
imgr = imrotate(img, -π/4) # output axes (-105:618, -105:618)

# Preserve the original axes
# Note that this is more efficient than `@view imrotate(img, π/4)[axes(img)...]`
imgr = imrotate(img, π/4, axes(img)) # output axes (1:512, 1:512)
```

By default, `imrotate` uses bilinear interpolation with constant fill value (`NaN` or `0`). You can,
for example, use the nearest interpolation and fill the new region with white pixels:

```julia
using Interpolations, ImageCore
imrotate(img, π/4, method=Constant(), fillvalue=oneunit(eltype(img)))
```

And with some inspiration, maybe fill with periodic values and tile the output together to
get a mosaic:

```julia
using Interpolations, ImageCore
imgr = imrotate(img, π/4, fillvalue = Periodic())
mosaicview([imgr for _ in 1:9]; nrow=3)
```
"""
function imrotate(img::AbstractArray{T}, θ::Real, inds::Union{Tuple, Nothing} = nothing; kwargs...) where T
    # TODO: expose rotation center as a keyword
    θ = floor(mod(θ,2pi)*typemax(Int16))/typemax(Int16) # periodic discretezation
    tform = recenter(RotMatrix{2}(θ), center(img))
    # Use the `nothing` trick here because moving the `autorange` as default value is not type-stable
    inds = isnothing(inds) ? autorange(img, inv(tform)) : inds
    warp(img, tform, inds; kwargs...)
end


text:
# There is a misplaced line break on line 522, we should get rid of it

dataset_path = joinpath(dirname(@__FILE__), "cylinder-bands.data.orig")

f = open(dataset_path)
data = readlines(f)
close(f)

f = open("$dataset_path.custom", "w")
for (i, line) in enumerate(data)
  line = rstrip(line)
  if i == 522
    write(f, "$line$(rstrip(data[i+1]))\n")
  elseif i != 523
    write(f, "$line\n")
  end
end
close(f)

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
@testset "Degree Vertex Cover" begin

    g3 = star_graph(5)
    for g in testgraphs(g3)
        c = @inferred(vertex_cover(g, DegreeVertexCover()))
        @test c == [1,]
    end

    g4 = complete_graph(5)
    for g in testgraphs(g4)
        c = @inferred(vertex_cover(g, DegreeVertexCover()))
        @test length(c)== 4 #All except one vertex
    end

    #path_graph(5) with additional edge 2-5
    g5 = Graph([0 1 0 0 0; 1 0 1 0 1; 0 1 0 1 0; 0 0 1 0 1; 0 1 0 1 0])
    for g in testgraphs(g5)
        c = @inferred(vertex_cover(g, DegreeVertexCover()))
        @test sort(c) == [2, 4]
    end
end


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("lattice", "singer"), x="VoicePart", y="Height", Geom.boxplot)