text:
export PlainLogicCircuit, PlainLogicLeafNode, PlainLogicInnerNode,
    PlainLiteralNode, PlainConstantNode, Plain⋀Node, Plain⋁Node

#####################
# Plain logic nodes without additional fields
#####################

"""
Root of the plain logic circuit node hierarchy
"""
abstract type PlainLogicCircuit <: LogicCircuit end

"""
A plain logical leaf node
"""
abstract type PlainLogicLeafNode <: PlainLogicCircuit end

"""
A plain logical inner node
"""
abstract type PlainLogicInnerNode <: PlainLogicCircuit end

"""
A plain logical literal leaf node, representing the positive or negative literal of its variable
"""
struct PlainLiteralNode <: PlainLogicLeafNode
    literal::Lit
end

"""
A plain logical constant leaf node, representing true or false
"""
struct PlainConstantNode <: PlainLogicInnerNode
    constant::Bool
end

"""
A plain logical conjunction node (And node)
"""
mutable struct Plain⋀Node <: PlainLogicInnerNode
    children::Vector{PlainLogicCircuit}
end

"""
A plain logical disjunction node (Or node)
"""
mutable struct Plain⋁Node <: PlainLogicInnerNode
    children::Vector{PlainLogicCircuit}
end

#####################
# traits
#####################

@inline GateType(::Type{<:PlainLiteralNode}) = LiteralGate()
@inline GateType(::Type{<:PlainConstantNode}) = ConstantGate()
@inline GateType(::Type{<:Plain⋀Node}) = ⋀Gate()
@inline GateType(::Type{<:Plain⋁Node}) = ⋁Gate()

#####################
# methods
#####################

"Get the children of a given inner node"
@inline children(n::PlainLogicInnerNode) = n.children

"""
    conjoin(arguments::Vector{<:PlainLogicCircuit})

Conjoins Plain LogicCircuits.
"""
function conjoin(arguments::Vector{<:PlainLogicCircuit};
                 reuse=nothing)
    @assert length(arguments) > 0
    reuse isa Plain⋀Node && children(reuse) == arguments && return reuse
    return Plain⋀Node(arguments)
end

function disjoin(arguments::Vector{<:PlainLogicCircuit};
                    reuse=nothing)
    @assert length(arguments) > 0
    reuse isa Plain⋁Node && children(reuse) == arguments && return reuse
    return Plain⋁Node(arguments)
end

negate(a::PlainLiteralNode) = compile(PlainLiteralNode, -a.literal)

# claim `PlainLogicCircuit` as the default `LogicCircuit` implementation
compile(::Type{LogicCircuit}, args...) =
    compile(PlainLogicCircuit, args...)

compile(::Type{<:PlainLogicCircuit}, b::Bool) =
    PlainConstantNode(b)

compile(::Type{<:PlainLogicCircuit}, l::Lit) =
    PlainLiteralNode(l)

function compile(::Type{<:PlainLogicCircuit}, circuit::LogicCircuit)
    f_con(n) = compile(PlainLogicCircuit, constant(n))
    f_lit(n) = compile(PlainLogicCircuit, literal(n))
    f_a(_, cns) = conjoin(cns)
    f_o(_, cns) = disjoin(cns)
    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, PlainLogicCircuit)
end

fully_factorized_circuit(::Type{LogicCircuit}, n::Int) =
    fully_factorized_circuit(PlainLogicCircuit, n)

function fully_factorized_circuit(::Type{<:PlainLogicCircuit}, n::Int)
    ors = map(1:n) do v
        v = Var(v)
        pos = compile(PlainLogicCircuit, var2lit(v))
        neg = compile(PlainLogicCircuit, -var2lit(v))
        pos | neg
    end
    and = conjoin(ors)
    disjoin([and]) # see logistic circuits bias term
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
# ---
# title: 232. Implement Queue using Stacks
# id: problem232
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Stack, Design
# link: <https://leetcode.com/problems/implement-queue-using-stacks/description/>
# hidden: true
# ---
#
# Implement a first in first out (FIFO) queue using only two stacks. The
# implemented queue should support all the functions of a normal queue (`push`,
# `peek`, `pop`, and `empty`).
#
# Implement the `MyQueue` class:
#
#   * `void push(int x)` Pushes element x to the back of the queue.
#   * `int pop()` Removes the element from the front of the queue and returns it.
#   * `int peek()` Returns the element at the front of the queue.
#   * `boolean empty()` Returns `true` if the queue is empty, `false` otherwise.
#
# **Notes:**
#
#   * You must use **only** standard operations of a stack, which means only `push to top`, `peek/pop from top`, `size`, and `is empty` operations are valid.
#   * Depending on your language, the stack may not be supported natively. You may simulate a stack using a list or deque (double-ended queue) as long as you use only a stack's standard operations.
#
# **Follow-up:** Can you implement the queue such that each operation is
# **[amortized](https://en.wikipedia.org/wiki/Amortized_analysis)** `O(1)` time
# complexity? In other words, performing `n` operations will take overall `O(n)`
# time even if one of those operations may take longer.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["MyQueue", "push", "push", "peek", "pop", "empty"]
#     [[], [1], [2], [], [], []]
#     **Output**
#     [null, null, null, 1, 1, false]
#
#     **Explanation**
#     MyQueue myQueue = new MyQueue();
#     myQueue.push(1); // queue is: [1]
#     myQueue.push(2); // queue is: [1, 2] (leftmost is front of the queue)
#     myQueue.peek(); // return 1
#     myQueue.pop(); // return 1, queue is [2]
#     myQueue.empty(); // return false
#
#
#
#
# **Constraints:**
#
#   * `1 <= x <= 9`
#   * At most `100` calls will be made to `push`, `pop`, `peek`, and `empty`.
#   * All the calls to `pop` and `peek` are valid.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
using CSV
using DataFrames
using DiscreteChoiceModels

data = CSV.read("../data/apollo_modeChoiceData.csv", DataFrame)
data = data[data.RP .== 1, :]

model = multinomial_logit(
    [  # define utility functions
        1 => Coef(:tt_car) * data.time_car + Coef(:cost) * data.cost_car,
        2 => Coef(:asc_bus) + Coef(:tt_bus) * data.time_bus +
            Coef(:access) * data.access_bus + Coef(:cost) * data.cost_bus,
        3 => Coef(:asc_air) + Coef(:tt_air) * data.time_air +
            Coef(:access) * data.access_air + Coef(:cost) * data.cost_air,
        4 => Coef(:asc_rail) + Coef(:tt_rail) * data.time_rail +
            Coef(:access) * data.access_rail + Coef(:cost) * data.cost_rail
    ],
    data.choice,
    availability=[
        1 => data.av_car .== 1,
        2 => data.av_bus .== 1,
        3 => data.av_air .== 1,
        4 => data.av_rail .== 1
    ]
)

println(DiscreteChoiceModels.summary(model))


text:
"""
**Momentum constructor**

```julia
    Momentum(; η::Real=0.01, γ::Real=0.9)
```

Algorithm :

```math
\\begin{align*}
v_t =& \\gamma v_{t-1} + \\eta g_t\\\\
\\Delta x_t =& v_t
\\end{align*}
```
"""
mutable struct Momentum <: Optimizer
    opt_type::String
    t::Int64
    η::Float64
    γ::Float64
    v_t::AbstractArray
end

## Note the algorithm seems flawed, \η should be 1-\β and a supplementary global learning rate would be nice

function Momentum(; η::Real=0.01, γ::Real=0.9)
    @assert η > 0.0 "η must be greater than 0"
    @assert γ > 0.0 "γ must be greater than 0"

    Momentum("Momentum", 0, η, γ, [])
end

params(opt::Momentum) = "η=$(opt.η), γ=$(opt.γ)"

function update(opt::Momentum, g_t::AbstractArray{T}) where {T<:Real}
    # resize squares of gradients
    if opt.t == 0
        opt.v_t = zero(g_t)
    end

    # update timestep
    opt.t += 1

    opt.v_t = opt.γ * opt.v_t + opt.η * g_t

    return opt.v_t
end


text:
"""
    struct CompilerInvocation <: Any
Hold a pointer to a `clang::CompilerInvocation` object.
"""
struct CompilerInvocation
    ptr::CXCompilerInvocation
end

Base.unsafe_convert(::Type{CXCompilerInvocation}, x::CompilerInvocation) = x.ptr
Base.cconvert(::Type{CXCompilerInvocation}, x::CompilerInvocation) = x


text:
function make_row_matrix(vals,matrix)
    return hcat([v*matrix for v in vals]...)
end

function interpolation_matrix(vals,ndofs)
    return make_row_matrix(vals,diagm(ones(ndofs)))
end

function vec_to_symm_mat_converter(dim::Z) where {Z<:Integer}
    if dim == 2
        E1 = @SMatrix [1.0 0.0
                       0.0 0.0
                       0.0 1.0]
        E2 = @SMatrix [0.0 0.0
                       0.0 1.0
                       1.0 0.0]
        return [E1,E2]
    elseif dim == 3
        E1 = @SMatrix [1.0   0.0   0.0
                       0.0   0.0   0.0
                       0.0   0.0   0.0
                       0.0   1.0   0.0
                       0.0   0.0   1.0
                       0.0   0.0   0.0]
        E2 = @SMatrix [0.0   0.0   0.0
                       0.0   1.0   0.0
                       0.0   0.0   0.0
                       1.0   0.0   0.0
                       0.0   0.0   0.0
                       0.0   0.0   1.0]
        E3 = @SMatrix [0.0   0.0   0.0
                       0.0   0.0   0.0
                       0.0   0.0   1.0
                       0.0   0.0   0.0
                       1.0   0.0   0.0
                       0.0   1.0   0.0]
        return [E1,E2,E3]
    else
        throw(ArgumentError("Expected dim ∈ {1,2} got dim = $dim"))
    end
end

function symmetric_tensor_dimension(dim)
    if dim == 2
        return 3
    elseif dim == 3
        return 6
    else
        throw(ArgumentError("Expected dim ∈ {2,3}, got dim = $dim"))
    end
end

function update_mass_matrix!(matrix,basis,quad,scale,ndofs)
    for (p,w) in quad
        vals = basis(p)
        N = interpolation_matrix(vals,ndofs)
        matrix .+= scale*N'*N*w
    end
end

function update_mass_matrix!(matrix,basis,quad,linemap::LineMap,scale,ndofs)
    detjac = determinant_jacobian(linemap)
    for (p,w) in quad
        vals = basis(linemap(p...))
        N = interpolation_matrix(vals,ndofs)
        matrix .+= N'*N*detjac*scale*w
    end
end

function update_mass_matrix!(matrix,basis,quad,map::InterpolatingPolynomial,
    scale,ndofs)

    @assert length(scale) == length(quad)
    for (idx,(p,w)) in enumerate(quad)
        vals = basis(map(p))
        detjac = determinant_jacobian(map,p)
        N = interpolation_matrix(vals,ndofs)
        matrix .+= N'*N*detjac*scale[idx]*w
    end
end

function update_mass_matrix_on_faces!(matrix,basis,facequads,facemaps,
    scale,ndofs)

    @assert length(facequads) == length(facemaps)
    @assert length(facequads) == length(scale)
    for (quad,map,s) in zip(facequads,facemaps,scale)
        update_mass_matrix!(matrix,basis,quad,
            map,s,ndofs)
    end
end

function mass_matrix(basis,quad,scale,ndofs)
    nf = number_of_basis_functions(basis)
    nfndofs = nf*ndofs
    matrix = zeros(nfndofs,nfndofs)
    update_mass_matrix!(matrix,basis,quad,scale,ndofs)
    return matrix
end

function mass_matrix_on_boundary(basis,facequads,facemaps,scale,ndofs,NF)

    totaldofs = ndofs*NF
    matrix = zeros(totaldofs,totaldofs)
    update_mass_matrix_on_faces!(matrix,basis,facequads,facemaps,scale,ndofs)

    return matrix
end

function mass_matrix_on_boundary(basis,facequads,facemaps,scale,ndofs)

    NF = number_of_basis_functions(basis)
    return mass_matrix_on_boundary(basis,facequads,facemaps,scale,ndofs,NF)
end

function mass_matrix_on_boundary(basis,facequads,facemaps,facescale,
    iquad,imap,iscale,ndofs,nf)

    totaldofs = ndofs*nf
    matrix = zeros(totaldofs,totaldofs)

    update_mass_matrix_on_faces!(matrix,basis,facequads,
        facemaps,facescale,ndofs)
    update_mass_matrix!(matrix,basis,iquad,imap,iscale,ndofs)

    return matrix

end

function linear_form!(vector,rhsvals,basis,cellmap,quad)
    ndofs,nq = size(rhsvals)
    @assert length(quad) == nq

    detjac = determinant_jacobian(cellmap)
    for (idx,(p,w)) in enumerate(quad)
        vals = basis(cellmap(p))
        N = interpolation_matrix(vals,ndofs)
        vector .+= N'*rhsvals[:,idx]*detjac*w
    end
end

function linear_form(rhsvals::M,basis,quad) where {M<:AbstractMatrix}
    ndofs,nq = size(rhsvals)
    @assert length(quad) == nq

    nf = number_of_basis_functions(basis)
    rhs = zeros(nf*ndofs)

    for (idx,(p,w)) in enumerate(quad)
        vals = basis(p)
        N = interpolation_matrix(vals,ndofs)
        rhs += N'*rhsvals[:,idx]*w
    end
    return rhs
end

function linear_form(rhsvals::V,basis,quad) where {V<:AbstractVector}

    nq = length(quad)
    extrhs = repeat(rhsvals,inner=(1,nq))
    return linear_form(extrhs,basis,quad)
end


text:
module ShapeFactory

export dispatch, dispaxxx, IDispatch

include("ShapeDispatch.jl")
include("ShapeVariant.jl")
include("ShapeSimplex.jl")

# ITypeInfo machinery.

function com_getfuncdesc(type_info::Ptr{ITypeInfo}, idx::Cuint)
    desc = Ref(Ptr{FUNCDESC}(0))
    func = unsafe_load(unsafe_load(Ptr{Ptr{Ptr{Cvoid}}}(type_info)), 6)
    ccall(func, Cint, (Ptr{ITypeInfo}, Cuint, Ref{Ptr{FUNCDESC}}), type_info, idx, desc)
    return desc[]
end

function com_getfuncname(type_info::Ptr{ITypeInfo}, fd::FUNCDESC)
    name = Ref(Ptr{UInt16}(0))
    func = unsafe_load(unsafe_load(Ptr{Ptr{Ptr{Cvoid}}}(type_info)), 13)
    ccall(func, Cint, (Ptr{ITypeInfo}, Clong, Ref{Ptr{UInt16}}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}),
            type_info, fd.memid, name, C_NULL, C_NULL, C_NULL)
    return Symbol(String(name[]))
end

function com_gettypeinfo(idsp::Ptr{IDispatch})
    info = Ref(Ptr{ITypeInfo}(0))
    func = unsafe_load(unsafe_load(Ptr{Ptr{Ptr{Cvoid}}}(idsp)), 5)
    ccall(func, Cint, (Ptr{IDispatch}, Cuint, Cint, Ref{Ptr{ITypeInfo}}), idsp, 0, LCID, info)
    return info[]
end

function com_gettypeattr(type_info::Ptr{ITypeInfo})
    attr = Ref(Ptr{TYPEATTR}(0))
    func = unsafe_load(unsafe_load(Ptr{Ptr{Ptr{Cvoid}}}(type_info)), 4)
    ccall(func, Cint, (Ptr{ITypeInfo}, Ref{Ptr{TYPEATTR}}), type_info, attr)
    data = unsafe_load(attr[])
    com_releasetypeattr(type_info, attr[])
    return data
end

function com_releasetypeattr(type_info::Ptr{ITypeInfo}, attr::Ptr{TYPEATTR})
    func = unsafe_load(unsafe_load(Ptr{Ptr{Ptr{Cvoid}}}(type_info)), 20)
    ccall(func, Cint, (Ptr{ITypeInfo}, Ptr{TYPEATTR}), type_info, attr)
end

# In order to call something like cat.ActiveDocument.Part.Update(), i need to distiguish between
# properties and methods, to do so i need some information about symbol, but i cant'get it
# exclusively by symbol name. FUNCDESC retrieved by index, conversion if memid to index handled by
# ITypeInfo2 interface which seems broken in CATIA COM Implementation.
# Therefore there is some dancing with tambourine.
#
function com_gettypedata(obj::Ptr{IDispatch})
    info = com_gettypeinfo(obj)
    attr = com_gettypeattr(info)

    get!(type_cache, hash(attr.guid)) do
        tdata = TypeData(info, Dict{Symbol, FuncData}())
        for n in 0:attr.cFuncs - 1
            fd = unsafe_load(com_getfuncdesc(info, Cuint(n)))
            name = com_getfuncname(info, fd)
            meth = get!(tdata.func, name) do
                fdata = FuncData(fd.memid, 0, Array{UInt16}(undef, fd.cParams))
                if fd.cParams > 0
                    for p in 0:fd.cParams - 1
                        fdata.param_flags[p + 1] = unsafe_load(Ptr{UInt16}(fd.lprgelemdescParam + 32 * p + 24))
                    end
                end
                fdata
            end
            meth.invoke_mask |= fd.invkind
        end
        tdata
    end
end

# IDispatch machinery

function com_invoke(obj::ComObject, name::Symbol, flags::UInt16, args::Ref{DISPPARAMS}, ret::Ptr{Variant})
    fptr = unsafe_load(unsafe_load(Ptr{Ptr{Ptr{Cvoid}}}(obj.ptr)), 7)
    ccall(fptr, Cint, (Ptr{IDispatch}, Clong, Ref{GUID}, UInt32, UInt16, Ref{DISPPARAMS}, Ptr{Variant}, Ptr{Cvoid}, Ptr{Cvoid}),
    obj.ptr, obj.dat.func[name].memid, IID_NULL, LCID, flags, args, ret, C_NULL, C_NULL)
end

function com_callmethod(obj::ComObject, name::Symbol, flag::UInt16 = 0x0001)
    var = Ref(Variant(0, 0, 0, 0, 0))
    arg = Ref(DISPPARAMS(C_NULL, C_NULL, 0, 0))
    com_invoke(obj, name, flag, arg, Ptr{Variant}(pointer_from_objref(var)))
    ret = value(var[])
    variant_free(var[])
    return ret
end

# this function can be simpler, but it must handle [out] arguments.
#
function com_callfunc(obj::ComObject, name::Symbol, args::Array{Variant})
    flags = obj.dat.func[name].param_flags
    packs = deepcopy(args)
    out_n = 0
    for n in 1:length(args)
        if (flags[n] & PARAMFLAG_OUT) > 0
            com_variantset(pointer(packs, n), Int128(0))
            if (packs[n].vt & VT_ARRAY) > 0
                unsafe_store!(Ptr{Cshort}(pointer(packs, n)), VT_BYREF | VT_ARRAY | VT_VARIANT)
                com_variantset(pointer(packs, n), Ptr{Ptr{SafeArray}}(pointer(args, n) + 8))
            else
                unsafe_store!(Ptr{Cshort}(pointer(packs, n)), VT_BYREF | VT_ARRAY | VT_VARIANT)
                com_variantset(pointer(packs, n), pointer(args, n))
            end
            out_n += 1
        end
    end
    var = Ref(Variant(0, 0, 0, 0, 0))
    arg = Ref(DISPPARAMS(pointer(reverse!(packs)), C_NULL, length(args), 0))
    com_invoke(obj, name, 0x0001, arg, Ptr{Variant}(pointer_from_objref(var)))

    if out_n > 0
        ret = value.(reverse!(push!(filter(v -> (v.vt & VT_BYREF) > 0 , packs), var[])))
    else
        ret = value(var[])
    end

    map(x -> variant_free.(x), [args, packs, [var[]]])
    return ret
end

function comcall(obj::ComObject, name::Symbol, args)
    if length(args) == 0
        return com_callmethod(obj, name)
    else
        return com_callfunc(obj, name, collect(Variant.(args)))
    end
end

function Base.getproperty(idsp::Ptr{IDispatch}, name::Symbol)
    dat = com_gettypedata(idsp)
    if !(name in keys(dat.func))
        error("SF: there is no member $(name) on COM object.")
    end
    obj = ComObject(idsp, dat)
    if (dat.func[name].invoke_mask & 0x1) > 0
        return (args...) -> comcall(obj, name, args)
    else
        return com_callmethod(obj, name, 0x0002)
    end
end

# special treatment of arrays
#
function Base.getproperty(idspa::Array{Ptr{IDispatch}}, name::Symbol)
    dat = com_gettypedata(idspa[1])
    if !(name in keys(dat.func))
        error("SF: there is no member $(name) on COM object.")
    end
    objs = ComObject.(idspa, Ref(dat))
    if (dat.func[name].invoke_mask & 0x1) > 0
        return (args...) -> comcall.(objs, name, Ref(args))
    else
        return com_callmethod.(objs, name, 0x0002)
    end
end

function Base.setproperty!(idsp::Ptr{IDispatch}, name::Symbol, value)
    dat = com_gettypedata(idsp)
    if !(name in keys(dat.func))
        error("SF: there is no member $(name) on COM object.")
    end
    obj = ComObject(idsp, dat)
    var = Ref(Variant(value))
    pid = [Int32(-3)]
    arg = Ref(DISPPARAMS(Ptr{Variant}(pointer_from_objref(var)), pointer(pid), 1, 1))
    com_invoke(obj, name, 0x0004, arg, Ptr{Variant}(C_NULL))
    variant_free(var[])
end

function dispatch(name::String)

    wname = push!(transcode(UInt16, name), '\0')
    ccall((:CoInitializeEx, "ole32.dll"), Cint, (Ptr{Cvoid}, Cuint), C_NULL, 0)

    guid = Ref(GUID(0, 0, 0, (0, 0, 0, 0, 0, 0, 0, 0)))
    ccall((:CLSIDFromProgID, "ole32.dll"), Cint, (Ptr{UInt16}, Ref{GUID}), wname, guid)

    pdisp = Ref(Ptr{IDispatch}(0))
    err = ccall((:CoCreateInstance, "ole32.dll"), Cint, (Ref{GUID}, Ptr{IUnknown}, Cuint, GUID, Ref{Ptr{IDispatch}}),
            guid, C_NULL, CLSCTX_SERVER, IID_IDispatch, pdisp)

    if (err != 0)
        error("Failed to create IDispatch")
    end

    return pdisp[]
end

dispaxxx(name::String) = Ptr{IDispaxxx}(dispatch(name))

end

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory


text:
function train(X)
  l = size(X,1)
  p = size(X,2)
  W = zeros(l,l)
  for i in 1:p
    W += X[:,i]*X[:,i]'
  end
  W -= p*eye(l,l)
  return W
end

iterHopfield(W::Array{Float64,2}, x::Vector{Float64}) = sign(W*x)

hamming(a,b) = sum(abs(a-b)/2)

function findmatch(Ps::Array{Float64,2}, v::Vector{Float64})
  l = size(Ps,2)
  min = 1
  minidx = -1
  for i in 1:l
    d = hamming(Ps[:,i], v)
    if min == 0 && d == 0
      warn("find match -- double match")
      return -1
    end
    if d < min
      min = d
      minidx = i
    end
  end
  return minidx
end

function searchHopfield(Ps::Array{Float64,2}, v0::Vector{Float64}; W::Array{Float64,2}=zeros(1,0), maxiter::Int=4)
  W2 = Union{}
  if size(W,2) == 0
    W2 = train(Ps)
  else
    W2 = W
  end
  v = deepcopy(v0)
  for i in 1:maxiter
    v1 = iterHopfield(W2,v)
    hamming(v,v1) > 0 ? v = v1 : return findmatch(Ps, v1)
  end
  return -1
end

function basicTestNet()
  N = 20
  a = sign(randn(N)-1)
  b = sign(randn(N)-1)
  c = sign(randn(N))
  Ps = [a';b';c']'
  a1 = deepcopy(a)
  a1[1] = -1*a[1]
  a1[6] = -1*a[6]
  @show searchHopfield(Ps, a1)
  b1 = deepcopy(b)
  b1[6] = -1*b[6]
  @show searchHopfield(Ps, b1)
  return true
end

basicTestNet()


text:
## Forward Rules

# Note that this is type piracy as the derivative should be NaN for x == y.
function ChainRulesCore.frule(
    (_, Δx, Δy), d::Distances.Euclidean, x::AbstractVector, y::AbstractVector
)
    Δ = x - y
    D = sqrt(sum(abs2, Δ))
    if !iszero(D)
        Δ ./= D
    end
    return D, dot(Δ, Δx) - dot(Δ, Δy)
end

## Reverse Rules Delta

function ChainRulesCore.rrule(dist::Delta, x::AbstractVector, y::AbstractVector)
    d = dist(x, y)
    function evaluate_pullback(::Any)
        return NoTangent(), ZeroTangent(), ZeroTangent()
    end
    return d, evaluate_pullback
end

function ChainRulesCore.rrule(
    ::typeof(Distances.pairwise), d::Delta, X::AbstractMatrix, Y::AbstractMatrix; dims=2
)
    P = Distances.pairwise(d, X, Y; dims=dims)
    function pairwise_pullback(::Any)
        return NoTangent(), NoTangent(), ZeroTangent(), ZeroTangent()
    end
    return P, pairwise_pullback
end

function ChainRulesCore.rrule(
    ::typeof(Distances.pairwise), d::Delta, X::AbstractMatrix; dims=2
)
    P = Distances.pairwise(d, X; dims=dims)
    function pairwise_pullback(::Any)
        return NoTangent(), NoTangent(), ZeroTangent()
    end
    return P, pairwise_pullback
end

function ChainRulesCore.rrule(
    ::typeof(Distances.colwise), d::Delta, X::AbstractMatrix, Y::AbstractMatrix
)
    C = Distances.colwise(d, X, Y)
    function colwise_pullback(::Any)
        return NoTangent(), NoTangent(), ZeroTangent(), ZeroTangent()
    end
    return C, colwise_pullback
end

## Reverse Rules DotProduct

function ChainRulesCore.rrule(dist::DotProduct, x::AbstractVector, y::AbstractVector)
    d = dist(x, y)
    function evaluate_pullback(Δ::Any)
        return NoTangent(), Δ .* y, Δ .* x
    end
    return d, evaluate_pullback
end

function ChainRulesCore.rrule(
    ::typeof(Distances.pairwise),
    d::DotProduct,
    X::AbstractMatrix,
    Y::AbstractMatrix;
    dims=2,
)
    P = Distances.pairwise(d, X, Y; dims=dims)
    function pairwise_pullback_cols(Δ::Any)
        if dims == 1
            return NoTangent(), NoTangent(), Δ * Y, Δ' * X
        else
            return NoTangent(), NoTangent(), Y * Δ', X * Δ
        end
    end
    return P, pairwise_pullback_cols
end

function ChainRulesCore.rrule(
    ::typeof(Distances.pairwise), d::DotProduct, X::AbstractMatrix; dims=2
)
    P = Distances.pairwise(d, X; dims=dims)
    function pairwise_pullback_cols(Δ::Any)
        if dims == 1
            return NoTangent(), NoTangent(), 2 * Δ * X
        else
            return NoTangent(), NoTangent(), 2 * X * Δ
        end
    end
    return P, pairwise_pullback_cols
end

function ChainRulesCore.rrule(
    ::typeof(Distances.colwise), d::DotProduct, X::AbstractMatrix, Y::AbstractMatrix
)
    C = Distances.colwise(d, X, Y)
    function colwise_pullback(Δ::Any)
        return NoTangent(), NoTangent(), Δ' .* Y, Δ' .* X
    end
    return C, colwise_pullback
end

## Reverse Rules Sinus

function ChainRulesCore.rrule(s::Sinus, x::AbstractVector, y::AbstractVector)
    d = x - y
    sind = sinpi.(d)
    abs2_sind_r = abs2.(sind) ./ s.r
    val = sum(abs2_sind_r)
    gradx = twoπ .* cospi.(d) .* sind ./ (s.r .^ 2)
    function evaluate_pullback(Δ::Any)
        return (r=-2Δ .* abs2_sind_r,), Δ * gradx, -Δ * gradx
    end
    return val, evaluate_pullback
end

## Reverse Rulse SqMahalanobis

function ChainRulesCore.rrule(
    dist::Distances.SqMahalanobis, a::AbstractVector, b::AbstractVector
)
    d = dist(a, b)
    function SqMahalanobis_pullback(Δ::Real)
        a_b = a - b
        ∂qmat = InplaceableThunk(
            X̄ -> mul!(X̄, a_b, a_b', true, Δ), @thunk((a_b * a_b') * Δ)
        )
        ∂a = InplaceableThunk(
            X̄ -> mul!(X̄, dist.qmat, a_b, true, 2 * Δ), @thunk((2 * Δ) * dist.qmat * a_b)
        )
        ∂b = InplaceableThunk(
            X̄ -> mul!(X̄, dist.qmat, a_b, true, -2 * Δ), @thunk((-2 * Δ) * dist.qmat * a_b)
        )
        return Tangent{typeof(dist)}(; qmat=∂qmat), ∂a, ∂b
    end
    return d, SqMahalanobis_pullback
end

## Reverse Rules for matrix wrappers

function ChainRulesCore.rrule(::Type{<:ColVecs}, X::AbstractMatrix)
    ColVecs_pullback(Δ::Tangent) = (NoTangent(), Δ.X)
    function ColVecs_pullback(::AbstractVector{<:AbstractVector{<:Real}})
        return error(
            "Pullback on AbstractVector{<:AbstractVector}.\n" *
            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`.\n" *
            "To solve this issue overload `kernelmatrix(_diag)` for your kernel for `ColVecs`",
        )
    end
    return ColVecs(X), ColVecs_pullback
end

function ChainRulesCore.rrule(::Type{<:RowVecs}, X::AbstractMatrix)
    RowVecs_pullback(Δ::Tangent) = (NoTangent(), Δ.X)
    function RowVecs_pullback(::AbstractVector{<:AbstractVector{<:Real}})
        return error(
            "Pullback on AbstractVector{<:AbstractVector}.\n" *
            "This might happen if you try to use gradients on the generic `kernelmatrix` or `kernelmatrix_diag`.\n" *
            "To solve this issue overload `kernelmatrix(_diag)` for your kernel for `RowVecs`",
        )
    end
    return RowVecs(X), RowVecs_pullback
end


text:
using LinearAlgebra
using MixedModels
using PooledArrays
using Random
using SparseArrays
using Test

using MixedModels: rankUpdate!

@testset "mul!" begin
    for (m, p, n, q, k) in (
                            (10, 0.7, 5, 0.3, 15),
                            (100, 0.01, 100, 0.01, 20),
                            (100, 0.1, 100, 0.2, 100),
                           )
        a = sprand(m, n, p)
        b = sprand(n, k, q)
        as = sparse(a')
        bs = sparse(b')
        c = zeros(m, k)
        ab = a * b
        arbt = Array(b')
        aab = Array(a) * Array(b)
        @test aab ≈ mul!(c, a, bs', true, true)
        @test aab ≈ mul!(c, a, bs')
        @test aab ≈ mul!(c, a, arbt')
        @test aab ≈ mul!(c, a, arbt')
        @test aab ≈ mul!(fill!(c, 0.0), a, arbt', true, true)
        @test maximum(abs, mul!(c, a, arbt', -1.0, true)) ≤ sqrt(eps())
        @test maximum(abs.(ab - aab)) < 100*eps()
        @test a*bs' == ab
        @test as'*b == ab
        @test as'*bs' == ab
        f = Diagonal(rand(n))
        @test Array(a*f) == Array(a)*f
        @test Array(f*b) == f*Array(b)
    end
end

@testset "reweight!" begin
    rng = MersenneTwister(1234321)
    df = (
        Y = randn(rng, 400),
        A = repeat(PooledArray(["N","Y"]), outer=200),
        G = repeat(PooledArray(string.('A':'T')), inner = 2, outer=10),
        H = repeat(PooledArray(string.('a':'j')), inner=40),
        )
    m1 = fit(MixedModel, @formula(Y ~ 1 + A + (1+A|G) + (1+A|H)), df)
    wm1 = fit(MixedModel, @formula(Y ~ 1+A+(1+A|G)+(1+A|H)), df, wts  = ones(400))
    @test loglikelihood(wm1) ≈ loglikelihood(m1)
    MixedModels.reweight!(wm1, ones(400))
    @test loglikelihood(refit!(wm1)) ≈ loglikelihood(m1)
end

@testset "rankupdate!" begin
    x = [1 1; 1 1];
    # in Julia 1.6+, typeof(x) == Matrix{Int64}
    # in < 1.6, typeof(x) == Array{Int64, 2}
    err = ErrorException("We haven't implemented a method for $(typeof(x)), $(typeof(x)). Please file an issue on GitHub.");
    @test_throws ErrorException rankUpdate!(x, x, 1, 1);
end

#=  I don't see this testset as meaningful b/c diagonal A does not occur after amalgamation of ReMat's for the same grouping factor - D.B.
@testset "rankupdate!" begin
    @test ones(2, 2) == rankUpdate!(Hermitian(zeros(2, 2)), ones(2))
    d2 = Diagonal(fill(2., 2))
    @test Diagonal(fill(5.,2)) == rankUpdate!(Diagonal(ones(2)), d2, 1.)
    @test Diagonal(fill(-3.,2)) == rankUpdate!(Diagonal(ones(2)), d2, -1.)

    # when converting straight from diagonal to symmetric, the type is different
    @test Diagonal(fill(5.,2)) == rankUpdate!(Symmetric(Matrix(1. * I(2)), :L), d2)
    # generic method
    @test Diagonal(fill(5.,2)) == rankUpdate!(Matrix(1. * I(2)), d2)
end
=#

@testset "lmulλ!" begin
    levs(ng, tag='S') = string.(tag, lpad.(string.(1:ng), ndigits(ng), '0'))

    function gendata(rng::AbstractRNG, n::Integer, ng::Integer, nh::Integer)
        (
            Y = randn(rng, n),
            X = rand(rng, n),
            G = PooledArray(rand(rng, levs(ng, 'G'), n)),
            H = PooledArray(rand(rng, levs(nh, 'H'), n)),
        )
    end
    gendata(n::Integer, ng::Integer, nh::Integer) = gendata(MersenneTwister(42), n, ng, nh)
    gendata(n::Integer, ng::Integer) = gendata(MersenneTwister(42), n, ng, ng)

    @testset "Adjoint{T,ReMat{T,1}}, BlockedSparse{T,1,2}" begin
        # this is an indirect test of lmulΛ! for a blocking structure found in
        # an example in MixedModels.jl#123
        df = gendata(10000, 500)
        f = @formula(Y ~ (1 + X|H) + (1|G))
        m500 = fit!(LinearMixedModel(f, df))
        # the real test here isn't in the theta comparison but in that the fit
        # completes successfully
        @test length(m500.u) == 2
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IPConfigurationProfilePropertiesFormat <: SwaggerModel
    subnet::Any # spec type: Union{ Nothing, Subnet2 } # spec name: subnet
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function IPConfigurationProfilePropertiesFormat(;subnet=nothing, provisioningState=nothing)
        o = new()
        validate_property(IPConfigurationProfilePropertiesFormat, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(IPConfigurationProfilePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type IPConfigurationProfilePropertiesFormat

const _property_map_IPConfigurationProfilePropertiesFormat = Dict{Symbol,Symbol}(Symbol("subnet")=>Symbol("subnet"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_IPConfigurationProfilePropertiesFormat = Dict{Symbol,String}(Symbol("subnet")=>"Subnet2", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ IPConfigurationProfilePropertiesFormat }) = collect(keys(_property_map_IPConfigurationProfilePropertiesFormat))
Swagger.property_type(::Type{ IPConfigurationProfilePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IPConfigurationProfilePropertiesFormat[name]))}
Swagger.field_name(::Type{ IPConfigurationProfilePropertiesFormat }, property_name::Symbol) =  _property_map_IPConfigurationProfilePropertiesFormat[property_name]

function check_required(o::IPConfigurationProfilePropertiesFormat)
    true
end

function validate_property(::Type{ IPConfigurationProfilePropertiesFormat }, name::Symbol, val)
end


text:
# lattice.jl
#
# Defines the spatial lattice; supports 1D and 2D square, open and periodic boundaries
# The main data structure is the bond-spin index array, bond_spin[nBond,2]

function lattice_bond_spins(nX::Int, pbc::Bool=true)
    Nb = pbc ? nX : nX-1
    Ns = nX

    bond_spin = [(0, 0) for _ in 1:Nb]
    for i in 1:Nb
        bond_spin[i] = (i, i+1)
    end

    if pbc
        bond_spin[end] = (Nb, 1)
    end

    return bond_spin, Ns, Nb
end

# 2D square lattice
function lattice_bond_spins(nX::NTuple{2, Int}, pbc::Bool=true)
    (nX[1] == nX[2]) || throw(ArgumentError("Currently only supports square lattices!"))
    L = nX[1]

    Ns = L * L
    Nb = pbc ? 2 * Ns : 2 * (Ns - L)
    bond_spin = [(0, 0) for _ in 1:Nb]

    if pbc
        @inbounds for i in 1:Ns
            # horizontal
            bond1 = 2 * i - 1

            if mod(i, L) != 0
                bond_spin[bond1] = (i, i + 1)
            else
                bond_spin[bond1] = (i, i + 1 - L)
            end

            # vertical
            bond2 = 2 * i

            if i <= (Ns - L)
                bond_spin[bond2] = (i, i + L)
            else
                bond_spin[bond2] = (i, i + L - Ns)
            end
        end
    else
        # horizontal
        cnt = 1
        @inbounds for i in 1:div(Nb, 2)
            if mod(i, L - 1) != 0
                bond_spin[i] = (cnt, cnt + 1)
                cnt += 1
            else
                bond_spin[i] = (cnt, cnt + 1)
                cnt += 2
            end
        end

        # vertical
        cnt = 1
        @inbounds for i in (div(Nb, 2)+1):Nb
            bond_spin[i] = (cnt, cnt + L)
            cnt += 1
        end
    end

    return bond_spin, Ns, Nb
end


text:
@inline function checktail(tail::Symbol)
    if tail ∉ (:left, :right, :both)
        throw(ArgumentError("Unsupported tail specifier ($tail)"))
    end
end

"""
    logpvalue(nisect::Integer, na::Integer, nb::Integer, ntotal::Integer,
              [tail::Symbol = :right])

Log P-value for the two sets intersection.

`A` has `na` elemnts, B has `nb` elements,
they have `nisect` elements in common, and there are
`ntotal` elements in the "universe".

`tail` controls the null hypothesis:
* `:right` (default): by chance `A` and `B` would have ≥ elements in common
* `:left`: by chance `A` and `B` would have ≤ elements in common
* `:both`: by chance `A` and `B` would have either ≤ or ≥ elements in common, whichever is less probable
"""
function logpvalue(nisect::Integer, na::Integer, nb::Integer, ntotal::Integer,
                   tail::Symbol = :right)
    ((na >= 0) && (nb >= 0) && (ntotal >= 0)) ||
        throw(ArgumentError("Sets with negative number of elements"))
    ((na <= ntotal) && (nb <= ntotal)) ||
        throw(ArgumentError("Sets bigger than total number of elements"))
    # corner cases
    isect_max = min(na, nb)
    isect_min = max(0, na + nb - ntotal)
    if tail == :right
        (nisect > isect_max) && return -Inf
        (nisect <= isect_min) && return 0.0
    elseif tail == :left
        (nisect >= isect_max) && return 0.0
        (nisect < isect_min) && return -Inf
    elseif tail != :both
        # FIXME corner cases when tail == :both
        throw(ArgumentError("Unsupported tail specifier ($tail)"))
    end
    # normal cases
    distr = Distributions.Hypergeometric(na, ntotal - na, nb)
    if tail == :right
        return nisect < isect_max ? logccdf(distr, nisect-1) : logpdf(distr, nisect)
    elseif tail == :left
        return nisect > isect_min ? logcdf(distr, nisect) : logpdf(distr, nisect)
    elseif tail == :both
        return log(2.0) + min(logcdf(distr, nisect), logccdf(distr, nisect-1), log(0.5))
    else
        return NaN # noop
    end
end

# size of the sets intersection
# the sets are represented by the sorted sets of their tiles
function _isect_size(set1_tiles::AbstractVector{Int},
                     set2_tiles::AbstractVector{Int},
                     tile_sizes::Vector{Int})
    (length(set1_tiles)==0 || length(set2_tiles)==0) && return 0
    res = 0
    i = 1
    tile1 = set1_tiles[i] # current tile of the 1st set
    j = 1
    tile2 = set2_tiles[j] # current tile of the 2nd set
    @inbounds while true
        if tile1 < tile2
            i += 1
            (i <= length(set1_tiles)) || break
            tile1 = set1_tiles[i]
        elseif tile1 > tile2
            j += 1
            (j <= length(set2_tiles)) || break
            tile2 = set2_tiles[j]
        else # intersecting tile, adjust the overlap
            res += tile_sizes[tile1]
            i += 1
            j += 1
            (i <= length(set1_tiles) && j <= length(set2_tiles)) || break
            tile1 = set1_tiles[i]
            tile2 = set2_tiles[j]
        end
    end
    return res
end

# returns set×set score matrix
# a×b score is the Fisher's test P-value of their intersection
# a×a self-intersections are processed in the same way
# (might be used later in cover problem construction to penalize the use of
# multiple variables referencing the same set)
function _setXset_scores(tileXset::SparseMaskMatrix, total_size::Int,
                         set_sizes::AbstractVector{Int},
                         tile_sizes::AbstractVector{Int},
                         used_sets = Base.OneTo(size(tileXset, 2)))
    nsets = length(used_sets)
    res = zeros(Float64, nsets, nsets)
    # fill upper triangle
    @inbounds for i in 1:nsets
        set1_ix = used_sets[i]
        set1_size = set_sizes[set1_ix]
        set1_size == 0 && continue # skip empty sets
        set1_tiles = view(tileXset, :, set1_ix)
        for j in i:nsets
            set2_ix = used_sets[j]
            set2_size = set_sizes[set2_ix]
            set2_size == 0 && continue # skip empty sets
            set2_tiles = view(tileXset, :, set2_ix)
            # one-sided Fisher's P-value, right tail
            isect_size = _isect_size(set1_tiles, set2_tiles, tile_sizes)
            res[i, j] = min(logpvalue(isect_size, set1_size, set1_size + set2_size - isect_size, total_size),
                            logpvalue(isect_size, set2_size, set1_size + set2_size - isect_size, total_size))
        end
    end
    # symmetrize: copy upper triangle to lower
    @inbounds for i in 1:nsets
        for j in (i+1):nsets
            res[j, i] = res[i, j]
        end
    end
    return res
end

"""
    set_relevance(nset_observed::Integer, nset::Integer,
                  nobserved::Integer, ntotal::Integer) -> Float64

Calculates the relevance weight of the set that contains `nset` elements,
`nset_observed` of which were present (not necessarily enriched) in the data
that identified `nobserved` elements out of all known (`ntotal`).
It is used by [`SetMosaic`](@ref) to penalize the sets,
which could not be observed in the data (e.g. biological processes or pathways
that involve proteins not expressed by the cells used in the experiments).

While for [`MaskedSetMosaic`](@ref) it's recommended to use the IDs of data entities
(e.g. protein group IDs for proteomic data) to correctly count the set sizes
and estimate enrichment; `set_relevance()` should use the counts derived from
the original IDs of the annotation database (e.g. UniProt accession codes).
Otherwise it's not possible to correctly estimate the number of elements that
belong to the given annotated set, but were not observed in the data.

The returned value is the probability that no more than `nset_observed`
elements were observed at random.
"""
set_relevance(nset_observed::Integer, nset::Integer,
              nobserved::Integer, ntotal::Integer) =
    cdf(Hypergeometric(nset, ntotal - nset, nobserved), nset_observed)


text:
using JLD, DataFrames, Base.Dates

#polygon representing manhattan
city = [(40.87851017592601,-73.927903175354);(40.87743937506905,-73.92168045043945);(40.87516792197097,-73.91871929168701);(40.87393481479329,-73.91193866729736);(40.87182549927143,-73.90970706939697);(40.866016421491814,-73.91168117523193);(40.859655059077824,-73.91764640808105);(40.84855366761019,-73.9266586303711);(40.84394377141786,-73.92949104309082);(40.83550228531863,-73.9339542388916);(40.82816381252365,-73.9339542388916);(40.81972031701224,-73.93318176269531);(40.808677200957106,-73.93335342407227);(40.801920506109774,-73.92803192138672);(40.7958778790764,-73.92820358276367);(40.78210123234386,-73.94107818603516);(40.77547182731999,-73.93953323364258);(40.7701418259051,-73.94588470458984);(40.75831029512206,-73.95669937133789);(40.744395800976775,-73.96888732910156);(40.71746884168838,-73.97077560424805);(40.70927151739564,-73.9764404296875);(40.7066689811733,-73.99497985839844);(40.697299008636755,-74.01008605957031);(40.69899090674371,-74.02381896972656);(40.74907763805906,-74.01351928710938);(40.77950154452169,-73.99532318115234);(40.81822635589172,-73.96751403808594);(40.853163243121564,-73.94845962524414)]

#keep points in the dataset only if inside manhattan
function point_inside_polygon(x,y,poly)
  n = length(poly)
  inside =false

  p1x,p1y = poly[1]
  for i in 0:n
    p2x,p2y = poly[i % n + 1]
    if y > min(p1y,p2y)
      if y <= max(p1y,p2y)
        if x <= max(p1x,p2x)
          if p1y != p2y
            xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
          end
          if p1x == p2x || x <= xinters
            inside = !inside
          end
        end
      end
    end
    p1x,p1y = p2x,p2y
  end
  return inside
end

PTIME, DTIME, PLONG, PLAT, DLONG, DLAT, P1, P2, P3, TIP, P4 = 6, 7 ,11,12,13,14,16,17,18,19,20

for j in 1:12
  println("Opening data/trip_$j.csv")
  f = open("data/trip_$j.csv")
  f2 = open("data/reduced_trip_$j.csv", "w")
  write(f2, "ptime, dtime, plong, plat, dlong, dlat, fare, tip\n")
  for (i, ln) in enumerate(eachline(f))
    if i == 1
      continue
    end
    if i%1000000 == 0
      println(i)
    end
    s = split(ln,",")
    if s[PLONG] == "" || s[PLAT] == "" || s[DLONG] == "" || s[DLAT] == ""
      continue
    end
    plong = parse(Float64, s[PLONG])
    plat  = parse(Float64, s[PLAT])
    dlong = parse(Float64, s[DLONG])
    dlat  = parse(Float64, s[DLAT])
    if !(point_inside_polygon(plat,plong,city) && point_inside_polygon(dlat,dlong,city))
      continue
    end
    fare = parse(Float64,s[P1]) + parse(Float64,s[P2]) + parse(Float64,s[P3]) + parse(Float64,s[P4])
    write(f2, string(s[PTIME],",",s[DTIME],",",s[PLONG],",",s[PLAT],",",s[DLONG],",",s[DLAT],",", fare,",", s[TIP],"\n"))
  end
  close(f)
  close(f2)
end
for f in filter(x -> startswith(x,"trip") ,readdir("data"))
  rm("data/$f")
end

max_stars_repo_name: semartin/TaxiDataExtraction

max_issues_repo_name: ericytung/TaxiDataTools

max_forks_repo_name: ericytung/TaxiDataTools


text:
############################   MAIN FUNCTIONS ###############################

# Main function that takes care of translating the model into source code
# of target language
function generate_code!(ode_model::dt.OdeModel; language = "julia",
                        unit_analysis = true, name = "autogenerated_model",
                        file = "autogenerated_model", jacobian = false,
                        sensitivities = false)

  # Create dictionaries containing the default values of the inputs (at the reference scale)
  # and the conversion coefficients for the units. This prepares the inputs for the
  # code generation functions
  coef_observed = OrderedDict{ASCIIString, Float64}()
  for i in ode_model.NamesObserved
    coef_observed[i] = ode_model.Equations[i].Units.f
  end
  named_states = OrderedDict{ASCIIString, Float64}()
  coef_states = OrderedDict{ASCIIString, Float64}()
  for (key,val) in ode_model.States
    named_states[key] = val.Value * val.Units.f
    coef_states[key] = val.Units.f
  end
  named_parameters = OrderedDict{ASCIIString, Float64}()
  coef_parameters = OrderedDict{ASCIIString, Float64}()
  for (key,val) in ode_model.Parameters
    named_parameters[key] = val.Value * val.Units.f
    coef_parameters[key] = val.Units.f
  end
  forcings = OrderedDict{ASCIIString,Any}()
  coef_forcings = OrderedDict{ASCIIString,Float64}()
  c = 0
  for (key,value) in ode_model.Forcings
      c += 1
      forcings[key] = (float(value.Time), float(value.Value)*value.Units.f)
      coef_forcings[key] = value.Units.f
  end

  # Write a file with the default inputs for the model in JSON format
  write_inputs!(named_states,coef_states, named_parameters, coef_parameters,
               forcings, coef_forcings, ode_model.NamesObserved,coef_observed,
               file)

  # Sort the equations of the model
  sorted_model = OdeSorted(ode_model)

  # Check the physical dimensions of the equations in the model (if demanded)
  unit_analysis &&  check_units(sorted_model)

  # Generate the rhs function (includes observed variables)
  # If DynamicType = DAE generate semiexplicit representation of dae model
  if (sorted_model.DynamicType == "ODE")
    model_function = create_derivatives(language, sorted_model, name)
    if(language == "cpp")
      observed_function = ""
    else
      observed_function = create_observed(language, sorted_model, name)
    end
  else
    error("Only DynamicType = \"ode\" is currently supported")
  end

  # Write down the code!
  write_code!(language, ode_model.DynamicType, model_function, observed_function, name, file)

  nothing

end

############################ WRITE TIME DERIVATIVES ################################

# This function writes the source for calculating the time derivatives and observed variables
# and it dispatches the correct function depending on the language argument
function create_derivatives(language, sorted_model, name)
  if (language == "julia")
    return create_derivatives_julia(sorted_model, name)
  elseif (language == "cpp")
    return create_derivatives_rcpp(language, sorted_model, name)
  else
    error("Currently only \"julia\" and \"cpp\" are supported")
  end
end

############################ WRITE INDEPENDENT OBSERVER ################################

# This function writes the source for calculating the time derivatives and observed variables
# and it dispatches the correct function depending on the language argument
function create_observed(language, sorted_model, name)
  if (language == "julia")
    return create_observed_julia(sorted_model, name)
  else
    error("Currently only \"julia\" and \"cpp\" are supported")
  end
end

############################### WRITE CODE FILES ###################################

# Write a JSON file with all the inputs of the model.
# This should be readable in every language
function write_inputs!(named_states,coef_states, named_parameters, coef_parameters,
                       forcings, coef_forcings, observed,coef_observed,
                       file)
  f = open("$(file)_inputs.json","w")
  inputs = {"states" => named_states,
            "coef_states" => coef_states,
            "parameters" => named_parameters,
            "coef_parameters" => coef_parameters,
            "forcings" => forcings,
            "coef_forcings" => coef_forcings,
            "observed" => observed,
            "coef_observed" => coef_observed}
  JSON.print(f, inputs)
  close(f)
end

# Write a source file with all the functions associated to the model in the target language
function write_code!(language, dynamic_type, model_function, observed_function, name, file)

  if (language == "julia")
    return write_code_julia!(dynamic_type, model_function, observed_function, name, file)
  elseif (language == "cpp")
    return write_code_rcpp!(dynamic_type, model_function, name, file)
  end

end


text:
using Test
using CrypticCrosswords

include("SemanticSimilarity/runtests.jl")
include("grammar.jl")
include("derivations.jl")
include("explanations.jl")
include("clues.jl")


text:
#=
Probability distributions useful in economics.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-04

References
----------

http://en.wikipedia.org/wiki/Beta-binomial_distribution

=#

"""
The Beta-Binomial distribution

##### Fields

- `n, a, b::Float64` The three paramters to the distribution

##### Notes

See also http://en.wikipedia.org/wiki/Beta-binomial_distribution

"""
type BetaBinomial
    n::Integer
    a::Real
    b::Real
end

# moments
mean(d::BetaBinomial) = d.n * d.a / (d.a + d.b)

std(d::BetaBinomial) = sqrt(var(d))

function var(d::BetaBinomial)
    n, a, b = d.n, d.a, d.b
    top = n*a*b * (a + b + n)
    btm = (a+b)^2.0 * (a+b+1.0)
    top / btm
end

function skewness(d::BetaBinomial)
    n, a, b = d.n, d.a, d.b
    t1 = (a+b+2*n) * (b - a) / (a+b+2)
    t2 = sqrt((1+a+b) / (n*a*b * (n+a+b)))
    t1 * t2
end

"""
Evaluate the pdf of the distributions at the points 0, 1, ..., k

##### Arguments

`d::BetaBinomial`: Instance of `BetaBinomial` type

##### Returns

- `probs::vector{Float64}`: pdf of the distribution `d`, at `0:d.k`

"""
function pdf(d::BetaBinomial)
    n, a, b = d.n, d.a, d.b
    k = 0:n
    binoms = Float64[binomial(n, i) for i in k]
    probs = binoms .* beta(k .+ a, n .- k .+ b) ./ beta(a, b)
end


text:
function find_results(root_dir::String = ".";
        feedback_function::Function = function(x) end)

    fun = function(path)
        if isdir(path)
            tree = Core.ResDirTree(path);
            if isfile(Core.datafile(tree)) &&
                isfile(Core.logfile(tree)) &&
                isfile(Core.descfile(tree))

                feedback_function(path);
                return true;
            else
                return false;
            end
        else
            return false;
        end
    end

    return MyTools.filter_from!(root_dir, fun);
end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
### Irving Gomez Mendez
### January 13th 2019

using Pkg
using Plots; pyplot()
using LaTeXStrings
using Distributions
using StatsPlots
using Random
using StatsBase

# Apply a function to each value of an array and returns a new array
# containing the resulting values:
roundmap(vect, b) = map(x -> round(x, digits = b), vect)

# La densidad de la distribución Weibull de mínimos (de 2 parámetros) está dada por la expresión:

L"f(t;\eta,\beta)=\frac{\beta}{\eta}\left(\frac{t}{\eta}\right)^{\beta-1}
e^{-(t/\eta)^\beta}1\!\!1\,\,_{(0<t<\infty)}"

# Grafica f(x;θ) para varios valores de θ. ¿Existe alguna combinación de los parámetros
# que logre que f sea simétrica alrededor de la moda?

# Definimos la función de densidad
# eta parametro de escala
# beta parametrode forma
dweib_min(x, eta = 1, beta = 1) = beta/eta*(x/eta)^(beta-1)*exp(-(x/eta)^beta)

# Hacemos una secuencia de valores para x
xx = range(0.001, stop = 3, length = 100)

# La siguiente linea mandara un error
yy = dweib_min(xx)

# La manera correcta de hacerlo es
yy = map(dweib_min, xx)
roundmap(yy,3)

plot(xx, yy,
    seriestype = :scatter,
    label = "",
    m = (0.9, :RebeccaPurple, 5),
    background_color = :Lavender
    )
plot!(dweib_min,
    line = (:tomato, 1, 2),
    label = "Weibull(1,1)")

plot(dweib_min, 0, 3,
    label = "Weibull(1,1)",
    background_color = :Lavender,
    lw = 2,
    color = :RebeccaPurple)

# StatsPlots.jl extends Distributions.jl by adding a type recipe for its
# distribution types, so they can be directly interpreted as plotting data:

plot(Weibull(1,1), 0, 3,
    label = "Weibull(1,1)",
    background_color = :Lavender,
    lw = 2,
    color = :RebeccaPurple)

# Si cambiamos el parámetro de forma a 2

plot(Weibull(2,1), 0, 3,
    label = "Weibull(1,2)",
    background_color = :Lavender,
    lw = 2,
    ylabel = L"f(x,θ)",
    color = :RebeccaPurple)

# Con varios valores de beta (2, 3, 5)

plot!(Weibull(3,1),
    label = "Weibull(1,3)",
    lw = 2,
    ls = :dash,
    color = :orange)
plot!(Weibull(5,1),
    label = "Weibull(1,5)",
    lw = 2,
    ls = :dot,
    color = :green)

# Obtener información sobre funciones y sobre atributos de plot

?plot
plotattr("linestyle")

## Teorema Integral de Probabilidad

#Si U∼Unifcont[0,1], entonces la v.a. F^(-1)(U) tiene la misma distribución que X

# Funcion cuantiles
qweib_min(x, η, β) = η*(log(1/(1-x)))^(1/β)

# Función que genera una muestra (simple) de la distribucion Weibull
function rweib_min(n, η, β)
    u = rand(n) #Generamos una muestra uniforme
    out = map(x -> qweib_min(x, η, β), u)
    return(out)
end

Random.seed!(111)
muestra = rweib_min(50, 2, 6.2)

histogram(muestra,
    label = "Empirical Weibull(2,6.2)",
    normalized = true,
    background_color = :Lavender,
    fillalpha = 0.33,
    color = :RebeccaPurple)

plot!(Weibull(6.2,2),
    label = "Weibull(2,6.2)",
    lw = 2)

Empirical = ecdf(muestra)

Empirical(1.23)

plot(sort(muestra), map(Empirical, sort(muestra)),
    label = "Empirical Weibull(2,6.2)",
    ylabel = L"F_n(x,θ)",
    background_color = :Lavender,
    linetype = :steppost,
    lw = 2,
    color = :RebeccaPurple)

pweib_min(x, η, β) = cdf(Weibull(β, η), x)

plot!(x -> pweib_min(x, 2, 6.2),
    label = "Weibull(2,6.2)",
    lw = 2)

muestra2 = rand(Weibull(6.2,2), 50)

max_stars_repo_name: IrvingGomez/academic-hugo

max_issues_repo_name: IrvingGomez/academic-hugo

max_forks_repo_name: IrvingGomez/academic-hugo


text:
using Test
using jInv.Mesh
using jInv.Utils

n      = [8 15 13]

# regular mesh
domain = [0.1 1.1 -1.2 1.4 0.5 1.3]
Mr     = getRegularMesh(domain,n)

# tensor mesh
h1     = rand(n[1])
h2     = rand(n[2])
h3     = rand(n[3])
Mt     = getTensorMesh3D(h1,h2,h3)

# conductivity
sig    = rand(prod(n))
Meshes = (Mr,Mt)

massMatrices = ( (getEdgeMassMatrix,getdEdgeMassMatrix),
                 (getNodalMassMatrix,getdNodalMassMatrix),
                 (getFaceMassMatrix,getdFaceMassMatrix))

for k=1:length(Meshes)
	for j=1:length(massMatrices)
		println("\ttest $(massMatrices[j][1]) on $(typeof(Meshes[k]))")
		M  = Meshes[k]
		Mf3 = massMatrices[j][1]
		dMf3 = massMatrices[j][2]
		Av = Mf3(M,sig);
		vt = randn(size(Av,1))

		chkDir, = checkDerivative(sig->Mf3(M,sig)*vt,sig->dMf3(M,sig,vt),sig,out=false)
		@test chkDir
	end
end


text:
include("table_utils.jl")

function simulate_waveforms( mcevents::TypedTables.Table, s::Simulation{T};
                             Δt::RealQuantity = 4u"ns",
                             max_nsteps::Int = 1000,
                             verbose = false ) where {T <: SSDFloat}
    n_total_physics_events = length(mcevents)
    Δtime = T(to_internal_units(internal_time_unit, Δt))
    n_contacts = length(s.detector.contacts)
    S = Val(get_coordinate_system(s.electron_drift_field.grid))
    contacts = s.detector.contacts;
    contact_ids = Int[]
    for contact in contacts if !ismissing(s.weighting_potentials[contact.id]) push!(contact_ids, contact.id) end end
    wps_interpolated = [ interpolated_scalarfield(s.weighting_potentials[id]) for id in contact_ids ];
    e_drift_field = get_interpolated_drift_field(s.electron_drift_field);
    h_drift_field = get_interpolated_drift_field(s.hole_drift_field);

    @info "Detector has $(n_contacts) contact(s)"
    @info "Table has $(length(mcevents)) physics events ($(sum(map(edeps -> length(edeps), mcevents.edep))) single charge depositions)."

    # First simulate drift paths
    drift_paths = _simulate_charge_drifts(mcevents, s, Δt, max_nsteps, e_drift_field, h_drift_field, verbose)
    # now iterate over contacts and generate the waveform for each contact
    @info "Generating waveforms..."
    waveforms = map(
        wp ->  map(
            x -> _generate_waveform(x.dps, to_internal_units.(internal_energy_unit, x.edeps), Δt, Δtime, wp, S),
            TypedTables.Table(dps = drift_paths, edeps = mcevents.edep)
        ),
        wps_interpolated
    )
    mcevents_chns = map(
        i -> add_column(mcevents, :chnid, fill(contact_ids[i], n_total_physics_events)),
        eachindex(contact_ids)
    )
    mcevents_chns = map(
        i -> add_column(mcevents_chns[i], :waveform, ArrayOfRDWaveforms(waveforms[i])),
        eachindex(waveforms)
    )
    return vcat(mcevents_chns...)
end

function _simulate_charge_drifts( mcevents::TypedTables.Table, s::Simulation{T},
                                  Δt::RealQuantity, max_nsteps::Int,
                                  e_drift_field::Interpolations.Extrapolation, h_drift_field::Interpolations.Extrapolation,
                                  verbose::Bool ) where {T <: SSDFloat}
    return @showprogress map(mcevents) do phyevt
        _drift_charges(s.detector, s.electron_drift_field.grid, s.point_types,
                        CartesianPoint{T}.(to_internal_units.(u"m", phyevt.pos)),
                        e_drift_field, h_drift_field,
                        T(Δt.val) * unit(Δt), max_nsteps = max_nsteps, verbose = verbose)
    end
end

function _generate_waveform( drift_paths::Vector{EHDriftPath{T}}, charges::Vector{<:SSDFloat}, Δt::RealQuantity, dt::T,
                             wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}}) where {T <: SSDFloat}
    timestamps = _common_timestamps( drift_paths, dt )
    timestamps_with_units = range(zero(Δt), step = Δt, length = length(timestamps))
    signal = zeros(T, length(timestamps))
    add_signal!(signal, timestamps, drift_paths, T.(charges), wp, S)
    RDWaveform( timestamps_with_units, signal )
end


text:
module HouseholdFrames

   include( "household_frames.jl" )

end


text:
abstract type ETSModel end
import Optim: optimize
"""
SES
Description:
    - Contains parameters and data to apply Simple Exponential Smoothing (SES)
    - Subtype of ETSModel, for future proofing
Constuctor:
    - Vector{AbstractFloat} y: contains observations for forecasting
    - Integer h: Stores the number of steps after time T to forecast (T = latest time in the observed data)
        - Default = 10, must be > 0
    - Union{Nothing, Float64} alpha: smoothing parameter [1]
        - 0 < alpha < 1
        - If no value is specified, an optimal value will be chosen using LBFGS.
    - Union{Nothing, Float64} init_level: initial level value [1]
        -Inf < init_level < Inf
        - If no value is specified, an optimal value will be chosen uisng LBFGS.

EXAMPLE USAGE:
    observations = [1.0, 2.0, 3.0]
    # Inputted values will be used if they are specified. Otherwise, they will be computed.
    mdl = SES(observations, h = 5)
    mdl = SES(observations, alpha = 0.4)
    mdl = SES(observations, alpha = 0.25, init_level = 500.0)
    mdl = SES(observations, h = 15, alpha = 0.3, init_level = 750.0)

"""
mutable struct SES <: ETSModel
    y::Vector{AbstractFloat}
    h::Integer
    alpha::Union{Nothing, Float64}
    init_level::Union{Nothing, Float64}
    function SES(y = []; h = 5, alpha = nothing, init_level = nothing) # Constuctor
        # Input Cleaning
        if length(y) <= 0
            error("The input array is empty.")
        end
        if h <= 0
            error("The number of prediction steps must be positive")
        end

        cleanParams_(alpha, "alpha")

        return new(y, h, alpha, init_level)
    end
end

"""
makeForecast_()
Description:
    - Helper function
    - computes a forecast with Simple Exponential Smoothing (SES) [1]
Parameters
    - model: SES struct (defined above)
    - Note: when makeForecast_ is called, all parameters have been calibrated.
Returns:
    A tuple containing the following:
    - Vector{Float64} forecast: contains fitted values
    - SSE::Float64 SSE: the Sum of the Squares (error term)
        - used for parameter optimzation
References:
    [1] Hyndman, R.J., & Athanasopoulos, G. (2019) *Forecasting:
        principles and practice*, 3rd edition, OTexts: Melbourne,
        Australia. OTexts.com/fpp3. Accessed on April 19th 2020.
"""
function makeForecast_(model::SES) # Return vector of fitted values of length h + SSE
    h = model.h
    alpha = model.alpha
    init_level = model.init_level
    y = model.y
    data_size = length(y)

    lvls = zeros(data_size + 1) # Stores smoothed values
    forecast = zeros(data_size + h) # Stores fitted values
    lvls[1] = init_level # Set l_0
    forecast[1] = init_level
    SSE = 0.0 # Compute Sum of Squared Errors

    for i in 2:(length(lvls)) # Compute l_t's
        lvls[i] = alpha * y[i-1] + (1 - alpha) * lvls[i-1];
        forecast[i] = lvls[i]
        SSE += (forecast[i-1] - y[i-1])^2
    end

    for i in 1:h # Set forecasted values
        forecast[data_size + i] = lvls[end]
    end

    return forecast, SSE
end

"""
SSE_()
Description:
    - Returns the Sum of Squared Errors for a given alpha and init_level.
Parameters:
    - model: SES struct (defined above)
    - alpha: value of alpha used for computation
    - init_level: value of init_level used for computation
    - verbose: if verbose > 0, then the parameters will be printed as the are optimized.
Returns:
    - Float64 res: the Sum of Squared Errors for the given alpha and init_level values.
"""
function SSE_(model::SES; alpha = nothing, init_level = nothing, verbose = 0) # Helper function for optimization
    model.alpha = alpha
    model.init_level = init_level

    res = makeForecast_(model)[2]
    if verbose > 0
        println("Error: $res, Alpha: $(model.alpha), L0: $(model.init_level)")
    end

    return res
end
"""
cleanParams_()
Description:
    Throws an error if a given parameters is out of the specified bounds.
Parameters:
    - param: value of parameter to be checked
    - name: name of parameter to be checked
    - lb: lower bound of parameter, default = 0.0
    - ub: upper bound of parameter, defualt = 1.0
Note: if param is not of type Nothing, then param must be in (0.0, 1.0)
Returns:
    - nothing
"""
function cleanParams_(param::Union{Nothing, Float64}, name::String, lb::Float64 = 0.0, ub::Float64 = 1.0)
    if !isnothing(param) && (param < lb || param > ub)
        error(name * " needs to be in the range ($lb, $ub)")
    end
end
"""
fit!(model::ETSModel)
Description:
    fits an ETSModel object
    - computes initial values of parameters: 'alpha' and 'init_value'
    - modifies the 'model' directly.
    - after fit! is called, predict(model) is ready to be called
Parameters:
    - model: SES struct
    - v: verbosity, default = 0. if v > 0, verbose will be used.
Computation of model parameters:
'alpha' and 'init_value' are computed using LBFGS
Returns:
    - nothing
"""
function fit(model::SES; v=0) # Set alpha + init_level
    # Creating a copy of the user inputted params for the model as model.alpha and model.init_level will be updated later
    base_alpha = model.alpha
    base_init_level = model.init_level

    # Model to return to ensure that input model is not modified
    retModel = deepcopy(model)

    if (isnothing(retModel.alpha))
        @warn "Since no value was entered for 'alpha', it will be chosen"
    end
    if (isnothing(retModel.init_level))
        # model.init_level = compute_init_heuristic_(model.y)
        @warn "Since no value was entered for 'init_level', it will be chosen"
    end

    # Optimize predict(model) - directly updates model
    f(x) = SSE_(retModel,
                alpha=(isnothing(base_alpha) ? x[1] : base_alpha),
                init_level=(isnothing(base_init_level) ? x[2] : base_init_level),
                verbose=v)

    if isnothing(retModel.alpha) || isnothing(retModel.init_level) # if atleast one value needs to be optimized
        # lower and upper limits for alpha and init_level
        lower = [0.0, -Inf]
        upper = [1.0, Inf]
        # Initial values for alpha and init_level
        init_inputs = [0.0, model.y[1]]
        res = optimize(f, lower, upper, init_inputs)
        if v > 0
            println(res)
        end
    end

    return retModel
end
"""
predict()
Description:
    - Exposed function
    - computes a forecast with Simple Exponential Smoothing (SES) [1]
    - calls makeForecast_ a helper function, which returns a vector of fitted values and SSE error term.
Parameters
    - model: SES struct (defined above)
    - Note: when predict is called, all parameters have been calibrated.
Returns:
    A tuple containing the following:
    - Vector{Float64} forecast: contains fitted values
References:
    [1] Hyndman, R.J., & Athanasopoulos, G. (2019) *Forecasting:
        principles and practice*, 3rd edition, OTexts: Melbourne,
        Australia. OTexts.com/fpp3. Accessed on April 19th 2020.
"""
function predict(model::ETSModel)
    return makeForecast_(model)[1]
end

### Holt's Method - implementing Holt's linear trend method
mutable struct Holt <: ETSModel
    y::Vector{AbstractFloat}
    h::Integer
    alpha::Union{Nothing, Float64}
    beta::Union{Nothing, Float64}
    phi::Union{Nothing, Float64}
    init_level::Union{Nothing, Float64}
    init_trend::Union{Nothing, Float64}
    damped::Bool
    function Holt(y = []; h = 5, alpha = nothing, beta = nothing, phi = nothing, init_level = nothing, init_trend = nothing, damped = false) # Constuctor
        # Input Cleaning
        if length(y) <= 0
            error("The input array is empty.")
        end
        if h <= 0
            error("The number of prediction steps must be positive")
        end

        cleanParams_(alpha, "alpha")
        cleanParams_(beta, "beta")
        cleanParams_(phi, "phi", 0.80, 0.98)

        return new(y, h, alpha, beta, phi, init_level, init_trend)
    end
end

function makeForecast_(model::Holt) # Return vector of fitted values of length h + SSE
    h = model.h
    alpha = model.alpha
    beta = model.beta
    phi = model.phi
    init_level = model.init_level
    init_trend = model.init_trend
    y = model.y
    data_size = length(y)

    lvls = zeros(data_size + 1) # Stores level values, indices (0, T)
    trends = zeros(data_size + 1) # Stores trend values, indices (0, T)
    forecast = zeros(data_size + h) # Stores forecasted values, yhat, indices (1, T+h)

    lvls[1] = init_level # Set l_0
    trends[1] = init_trend # Set b_0

    SSE = 0.0 # Compute Sum of Squared Errors

    for i in 2:(length(lvls)) # Compute l_t's
        lvls[i] = alpha * y[i-1] + (1 - alpha) * (lvls[i-1] + trends[i-1]) # level equation
        trends[i] = beta * (lvls[i] - lvls[i-1]) + (1 - beta) * trends[i-1] # trned equation
        forecast[i-1] = lvls[i-1] + 1 * trends[i-1] # one-step ahead forecast
        SSE += (forecast[i-1] - y[i-1])^2
    end

    for i in 1:h # Set forecasted values
        forecast[data_size + i] = lvls[end] + i * trends[end]
    end

    return forecast, SSE
end

function SSE_(model::Holt; alpha = nothing, beta = nothing, phi = nothing, init_level = nothing, init_trend = nothing, verbose = 0) # Helper function for optimization
    model.alpha = alpha
    model.beta = beta
    model.phi = phi
    model.init_level = init_level
    model.init_trend = init_trend

    res = makeForecast_(model)[2]
    if verbose > 0
        println("Error: $res, Alpha: $(model.alpha), Beta*: $(model.beta), Phi: $(model.phi), L0: $(model.init_level), BO:$(model.init_trend)")
    end

    return res
end

function fit(model::Holt; v=0) # Set alpha + init_level
    # Creating a copy of the user inputted params for the model as model.alpha and model.init_level will be updated later
    base_alpha = model.alpha
    base_beta = model.beta
    base_phi = model.phi
    base_init_level = model.init_level
    base_init_trend = model.init_trend

    # copy of model to return to ensure that input model is not modified
    retModel = deepcopy(model)

    if (isnothing(retModel.alpha))
        @warn "Since no value was entered for 'alpha', it will be chosen"
    end
    if (isnothing(retModel.beta))
        @warn "Since no value was entered for 'beta', it will be chosen"
    end
    if (isnothing(retModel.init_level))
        @warn "Since no value was entered for 'init_level', it will be chosen"
    end
    if (isnothing(retModel.init_trend))
        @warn "Since no value was entered for 'init_trend', it will be chosen"
    end

    if damped
        print("The Damped Trend Method will be used.")
        if (isnothing(retModel.phi))
            @warn "Since no value was entered for `phi`, it will be chosen"
        end
    else
        if (!isnothing(retModel.phi))
            @warn "Since a value was entered for `phi`, damping will occur despite 'damped' being set to false"
        else
            base_phi = 1 # No damping is equivalent to setting the damping parameter to 1
        end
    end

    # Optimize makeForecast_(model) - directly updates model
    f(x) = SSE_(retModel,
                alpha=(isnothing(base_alpha) ? x[1] : base_alpha),
                beta=(isnothing(base_beta) ? x[2] : base_beta),
                phi=(isnothing(base_phi) ? x[3] : base_phi),
                init_level=(isnothing(base_init_level) ? x[4] : base_init_level),
                init_trend=(isnothing(base_init_trend) ? x[5] : base_init_trend),
                verbose=v)

    if isnothing(retModel.alpha) || isnothing(retModel.init_level) || isnothing(retModel.init_trend) || isnothing(retModel.beta) || (damped && isnothing(retModel.phi))# if at least one value needs to be optimized
        # lower and upper limits for alpha, beta, phi, init_level, init_trend
        lower = [0.0, 0.0, 0.8, -Inf, -Inf]
        upper = [1.0, 1.0, 0.98, Inf, Inf]

        # Initial values for alpha, beta, init_level, init_trend (optimization starting points)
        first_l0 = model.y[1]
        first_b0 = model.y[1]
        first_alpha = 0.0
        first_beta = 0.0
        first_phi = 0.9

        # Optimize the SSE of the model through the function "f"
        init_inputs = [first_alpha, first_beta, first_phi, first_l0, first_b0]
        res = optimize(f, lower, upper, init_inputs)

        if v > 0
            println(res)
        end
    end

    return retModel
end

### HoltWinters 's Method - implementing the HoltWinters Additive method
mutable struct HoltWinters <: ETSModel
    y::Vector{AbstractFloat}
    h::Integer
    alpha::Union{Nothing, Float64}
    beta::Union{Nothing, Float64}
    gamma::Union{Nothing, Float64}
    m::Integer
    init_level::Union{Nothing, Float64}
    init_trend::Union{Nothing, Float64}
    init_season::Vector{AbstractFloat}
    function HoltWinters(y = []; h = 5, alpha = nothing, beta = nothing, gamma = nothing, m, init_level = nothing, init_trend = nothing, init_season) # Constuctor
        # Input Cleaning
        if length(y) <= 0
            error("The input array is empty.")
        end
        if h <= 0
            error("The number of prediction steps must be positive")
        end
        if isnothing(m)
            error("You must provide a seasonal period 'm'")
        end
        if m <= 1
            error("'m' must be an integer greater than 1")
        end
        if length(init_season) != m
            error("The length of 'init_season' must be equal to `m`")
        end

        cleanParams_(alpha, "alpha")
        cleanParams_(beta, "beta")
        cleanParams_(gamma, "gamma")

        return new(y, h, alpha, beta, gamma, m, init_level, init_trend, init_season)
    end
end

function makeForecast_(model::HoltWinters) # Return vector of fitted values of length h + SSE
    h = model.h
    alpha = model.alpha
    beta = model.beta
    gamma = model.gamma
    m = model.m
    init_level = model.init_level
    init_trend = model.init_trend
    init_season = model.init_season
    y = model.y
    data_size = length(y)

    lvls = zeros(data_size + 1) # Stores level values, indices (0, T)
    trends = zeros(data_size + 1) # Stores trend values, indices (0, T)
    seasons = zeros(data_size + m) # Stores seasonal values
    forecast = zeros(data_size + h) # Stores forecasted values, yhat, indices (1, T+h)

    lvls[1] = init_level # Set l_0
    trends[1] = init_trend # Set b_0
    seasons[1:m] = init_season # Set s_0 ... s_m-1

    SSE = 0.0 # Compute Sum of Squared Errors
    k = trunc(Integer, (h-1)/m)

    for i in 2:(length(lvls)) # Compute l_t's
        season_idx = i + m - 1 # to start updating values after index m

        lvls[i] = alpha * (y[i-1] - seasons[season_idx - m]) + (1 - alpha) * (lvls[i-1] + trends[i-1]) # level equation
        trends[i] = beta * (lvls[i] - lvls[i-1]) + (1 - beta) * trends[i-1] # trend equation
        seasons[season_idx] = gamma * (y[i-1] - lvls[i-1] - trends[i-1]) + (1 - gamma) * seasons[season_idx - m] # season equation
        forecast[i-1] = lvls[i-1] + 1 * trends[i-1] + season[season_idx + 1 - m * (k+1)] # one-step ahead forecast
        SSE += (forecast[i-1] - y[i-1])^2
    end

    for i in 1:h # Set forecasted values
        forecast[data_size + i] = lvls[end] + i * trends[end + i - m * (k+1)]
    end

    return forecast, SSE
end

function SSE_(model::HoltWinters; alpha = nothing, beta = nothing, gamma = nothing, init_level = nothing, init_trend = nothing, verbose = 0) # Helper function for optimization
    model.alpha = alpha
    model.beta = beta
    model.gamma = gamma
    model.init_level = init_level
    model.init_trend = init_trend

    res = makeForecast_(model)[2]
    if verbose > 0
        println("Error: $res, Alpha: $(model.alpha), Beta*: $(model.beta), Gamma: $(model.gamma), L0: $(model.init_level), BO:$(model.init_trend), S0:$(model.init_season)")
    end

    return res
end

function fit(model::HoltWinters; v=0) # Set alpha + init_level
    # Creating a copy of the user inputted params for the model as model.alpha and model.init_level will be updated later
    base_alpha = model.alpha
    base_beta = model.beta
    base_gamma = model.gamma
    base_init_level = model.init_level
    base_init_trend = model.init_trend

    # copy of model to return to ensure that input model is not modified
    retModel = deepcopy(model)

    if (isnothing(retModel.alpha))
        @warn "Since no value was entered for 'alpha', it will be chosen"
    end
    if (isnothing(retModel.beta))
        @warn "Since no value was entered for 'beta', it will be chosen"
    end
    if (isnothing(retModel.gamma))
        @warn "Since no value was entered for 'gamma', it will be chosen"
    end
    if (isnothing(retModel.init_level))
        @warn "Since no value was entered for 'init_level', it will be chosen"
    end
    if (isnothing(retModel.init_trend))
        @warn "Since no value was entered for 'init_trend', it will be chosen"
    end

    # Optimize makeForecast_(model) - directly updates model
    f(x) = SSE_(retModel,
                alpha=(isnothing(base_alpha) ? x[1] : base_alpha),
                beta=(isnothing(base_beta) ? x[2] : base_beta),
                gamma=(isnothing(base_gamma) ? x[3] : base_gamma),
                init_level=(isnothing(base_init_level) ? x[4] : base_init_level),
                init_trend=(isnothing(base_init_trend) ? x[5] : base_init_trend),
                verbose=v)

    if isnothing(retModel.alpha) || isnothing(retModel.beta) || (isnothing(retModel.gamma)) || isnothing(retModel.init_level) || isnothing(retModel.init_trend) # if at least one value needs to be optimized
        # lower and upper limits for alpha, beta, gamma, init_level, init_trend
        lower = [0.0, 0.0, 0.0, -Inf, -Inf]
        upper = [1.0, 1.0, 1.0, Inf, Inf]

        # Initial values for alpha, beta, init_level, init_trend (optimization starting points)
        first_l0 = model.y[1]
        first_b0 = model.y[1]
        first_alpha = 0.0
        first_beta = 0.0
        first_gamma = 0.0

        # Optimize the SSE of the model through the function "f"
        init_inputs = [first_alpha, first_beta, first_gamma, first_l0, first_b0]
        res = optimize(f, lower, upper, init_inputs)

        if v > 0
            println(res)
        end
    end

    return retModel
end


text:
using FEASTSolver
using LinearAlgebra
using MatrixMarket
using Random
import FEASTSolver: iter_debug_print, update_R!, residuals, normalize!

function info(Λ, X, residuals, c, r)
    inside(x) = in_contour(x, c, r)
    in_eig = Λ[inside.(Λ)]
    in_res = residuals[inside.(Λ)]
    print("number inside : ")
    println(size(Λ[inside.(Λ)])[1])
    if sum(inside.(Λ)) > 0
        in_res_conv = in_res[in_res .<= 1e-5]
        in_eig_conv = in_eig[in_res .<= 1e-5]
        print("number inside converged : ")
        println(size(in_eig_conv)[1])
        print("max res inside: ")
        println(maximum(residuals[inside.(Λ)]))
        if size(in_res_conv, 1) > 0
            print("max res inside non spurious: ")
            println(maximum(in_res_conv))
        end
    end
    # display(in_eig_conv)
end

function moments_expand!(T, X::AbstractMatrix{ComplexF64}, nodes::Integer, iter::Integer;
    c=complex(0.0,0.0), r=1.0, debug=false, ϵ=10e-12, moments=2, spurious=1e-5)

    N, m₀ = size(X)
    Λ = zeros(ComplexF64, moments*m₀)
    θ = LinRange(π/nodes, 2*π-π/nodes, nodes)
    Q = zeros(ComplexF64, 2*moments, N, m₀)
    Q₀, Q₁ = zeros(ComplexF64, moments*N, moments*m₀), zeros(ComplexF64, moments*N, moments*m₀)
	Y = zeros(ComplexF64, N, m₀*moments)
    R = similar(Y, ComplexF64)

    for i=1:nodes
        z = (r*exp(θ[i]*im)+c)
		Temp = (T(z)\X) .* (r*exp(θ[i]*im)/nodes)
        Q[1,:,:] .+= Temp
        for j=1:2*moments
            Q[j,:,:] .+= Temp .* z^(j-1)
        end
    end
    for i=1:moments, j=1:moments
        Q₀[(i-1)*N+1:i*N, (j-1)*m₀+1:j*m₀] .= Q[i+j-1,:,:]
        Q₁[(i-1)*N+1:i*N, (j-1)*m₀+1:j*m₀] .= Q[i+j,:,:]
    end

    S = svd!(Q₀)
    F = eigen!(S.U' * Q₁ * S.V * Diagonal(1 ./ S.S))
    Y = S.U[1:N,:] * F.vectors
    Λ = F.values

    update_R!(Y, R, Λ, T)
    res = residuals(R, Λ, T)

	iter_debug_print(0, Λ, res, c, r, spurious)

	m₁ = moments*m₀
    Q₀, Q₁ = zeros(ComplexF64, moments*N, moments*m₁), zeros(ComplexF64, moments*N, moments*m₁)
	Q = zeros(ComplexF64, 2*moments, N, m₁)

    for i=1:nodes
        z = (r*exp(θ[i]*im)+c)
        resolvent = (1 ./(z .- Λ)) .* (r*exp(θ[i]*im)/nodes)

		Temp = (Y - T(z)\R) * Diagonal(resolvent)

        Q[1,:,:] .+= Temp
        for j=1:2*moments
            Q[j,:,:] .+= Temp .* z^(j-1)
        end
    end

    for i=1:moments, j=1:moments
        Q₀[(i-1)*N+1:i*N, (j-1)*m₁+1:j*m₁] .= Q[i+j-1,:,:]
        Q₁[(i-1)*N+1:i*N, (j-1)*m₁+1:j*m₁] .= Q[i+j,:,:]
    end

	S = svd!(Q₀)
    F = eigen!(S.U' * Q₁ * S.V * Diagonal(1 ./ S.S))
    Y = S.U[1:N,:] * F.vectors
    Λ = F.values

	R = similar(Y, ComplexF64)

	update_R!(Y, R, Λ, T)
    res = residuals(R, Λ, T)

	iter_debug_print(1, Λ, res, c, r, spurious)

	p = sortperm(res)
	res .= res[p]
	X .= Y[:, p][:,1:m₀]
	Λ .= Λ[p]
	R .= R[:, p]

    normalize!(X)
    Λ, X, res
end

#
# A0 = Matrix(mmread("data/quadraticM0.mtx"))
# A1 = Matrix(mmread("data/quadraticM1.mtx"))

Random.seed!(12345)
A0 = complex.(rand(100,100))
A1 = complex.(rand(100,100))
A0[:,1] .= 0

function T(z::ComplexF64)
	return (z + 0.4)*(z - 0.5)*(z + 0.4im)*(z-0.5im).*A1 .+ A0
end

R = 0.525
# C = complex(-0.8, 0.8)
C = complex(0, 0.0)

# e, v, res = nlfeast!(T, rand(ComplexF64,15,8), 2^6, 0, c=C, r=R, spurious=1e-3,debug=true)

e, v, res = nlfeast_moments!(T, rand(ComplexF64,100,30), 2^3, 10, c=C, r=R, moments=5, debug=true, ϵ=10e-16, store=false, spurious=1e-5)

# e, v, res = moments_expand!(T, rand(ComplexF64,15,7), 2^6, 1, c=C, r=R, moments=2, ϵ=10e-16, spurious=1e-3)
# e, v, res = beyn(T, A0, rand(ComplexF64,1000,120), 2^9; c=complex(-1.55,0.0), r=0.05)
# display(e)
# display(res)
#
inside(x) = in_contour(x, C, R)
# # print("\nmax res inside: ")
# # println(maximum(res[inside.(e)]))
# # print("number inside : ")
# # println(size(e[inside.(e)])[1])
# # # display(e)
# # display(res[inside.(e)])
#
# in_eig = e[inside.(e)]
# in_res = res[inside.(e)]
# in_res_conv = in_res[in_res .<= 1e-3]
# in_eig_conv = in_eig[in_res .<= 1e-3]
# if size(in_eig, 1) > 0
#     print("\nmax res inside: ")
#     println(maximum(res[inside.(e)]))
# end
# print("max res inside non spurious: ")
# println(maximum(in_res_conv))
# print("\nnumber inside : ")
# println(size(e[inside.(e)])[1])
# print("number inside converged : ")
# println(size(in_eig_conv)[1])

# Z = zeros(100,100)
# e, v, res = companion([(A0 .- 0.04.*A1), ((0.02im-0.02).*A1), (0.01im.*A1), ((-0.1-0.1im).*A1), A1])

info(e, v, res, C, R)

# display(res)
# println()
display(e[inside.(e)])


text:
using Test, OrdinaryDiffEq, HDF5, ModelConstructors
include(joinpath(dirname(@__FILE__), "../../../src/includeall.jl"))

regen_output = false

m = BruSan()
stategrid, _ = initialize_nojump!(m)
ode_f, ode_callback = eqcond_nojump(m)
tspan = (stategrid[:η][1], stategrid[:η][end])
θ = parameters_to_named_tuple(m.parameters)
prob = ODEProblem(ode_f, get_setting(m, :boundary_conditions)[:q][1], tspan, θ, tstops = stategrid[:η][2:end - 1])
sol = solve(prob, get_setting(m, :ode_integrator),
            reltol = get_setting(m, :ode_reltol),
            abstol = get_setting(m, :ode_abstol), callback = ode_callback)

@info "The following warning is expected."
stategrid, funcvar, derivs, endo = solve(m; nojump = true)

if regen_output
    h5open(joinpath(dirname(@__FILE__), "../../reference/models/brusan/brusan_nojump_eqm_log.h5"), "w") do file
        write(file, "q_Tsit5", sol.u)
        write(file, "q_Tsit5_solve", funcvar[:q])
        for (k, v) in endo
            write(file, "$(detexify(k))_Tsit5_solve", v)
        end
    end
end

true_soln = h5read(joinpath(dirname(@__FILE__), "../../reference/models/brusan/brusan_nojump_eqm_log.h5"), "q_Tsit5")
@testset "Solution for no jump equilibrium, log utility" begin
    @test @test_matrix_approx_eq sol.u true_soln

    file = joinpath(dirname(@__FILE__), "../../reference/models/brusan/brusan_nojump_eqm_log.h5")
    for (k, v) in endo
        @test @test_matrix_approx_eq h5read(file, "$(detexify(k))_Tsit5_solve") v
    end
end


text:
using KitBase, LinearAlgebra, Optim
using ProgressMeter: @showprogress

# one-cell simplification
begin
    quadratureorder = 5
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)
    L = 1
    ne = (L + 1)^2

    α = zeros(ne)
    u0 = [2.0, 0.0, 0.0, 0.0]
    m = KitBase.eval_spherharmonic(points, L)

    res = KitBase.optimize_closure(α, m, weights, u0, KitBase.maxwell_boltzmann_dual)
    u = KitBase.realizable_reconstruct(
        res.minimizer,
        m,
        weights,
        KitBase.maxwell_boltzmann_dual_prime,
    )
end

# multi-cell case
begin
    # space
    x0 = -1.5
    x1 = 1.5
    y0 = -1.5
    y1 = 1.5
    nx = 50#100
    ny = 50#100
    dx = (x1 - x0) / nx
    dy = (y1 - y0) / ny
    pspace = KitBase.PSpace2D(x0, x1, nx, y0, y1, ny)

    # time
    tEnd = 1.0
    cfl = 0.5
    dt = cfl / 2 * (dx * dy) / (dx + dy)

    # quadrature
    quadratureorder = 5
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)
    vspace = KitBase.VSpace1D{Float64,Int64,typeof(points),typeof(weights)}(
        -1.0,
        1.0,
        nq,
        points,
        zero(points),
        weights,
    )

    # particle
    SigmaS = 1.0 * ones(ny, nx)
    SigmaA = 0.0 * ones(ny, nx)
    SigmaT = SigmaS + SigmaA

    # moments
    L = 1
    ne = (L + 1)^2
    phi = zeros(ne, nx, ny)
    α = zeros(ne, nx, ny)
    m = KitBase.eval_spherharmonic(points, L)
end

# initial distribution
begin
    s2 = 0.03^2
    flr = 1e-4
    init_field(x, y) = max(flr, 1.0 / (4.0 * pi * s2) * exp(-(x^2 + y^2) / 4.0 / s2))
    for j = 1:nx
        for i = 1:ny
            y = y0 + (i - 0.5) * dy
            x = x0 + (j - 0.5) * dx
            # only zeroth order moment is non-zero
            phi[1, i, j] = init_field(x, y)
        end
    end
end

global t = 0.0
flux1 = zeros(ne, nx + 1, ny)
flux2 = zeros(ne, nx, ny + 1)

@showprogress for iter = 1:20
    # regularization
    Threads.@threads for j = 1:ny
        for i = 1:nx
            res = KitBase.optimize_closure(
                α[:, i, j],
                m,
                vspace.weights,
                phi[:, i, j],
                KitBase.maxwell_boltzmann_dual,
            )
            α[:, i, j] .= res.minimizer

            phi[:, i, j] .= KitBase.realizable_reconstruct(
                res.minimizer,
                m,
                vspace.weights,
                KitBase.maxwell_boltzmann_dual_prime,
            )
        end
    end

    # flux
    fη1 = zeros(nq)
    for j = 1:ny
        for i = 2:nx
            KitBase.flux_kfvs!(
                fη1,
                KitBase.maxwell_boltzmann_dual.(α[:, i-1, j]' * m)[:],
                KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:],
                vspace.u[:, 1],
                dt,
            )

            for k in axes(flux1, 1)
                flux1[k, i, j] = sum(m[k, :] .* vspace.weights .* fη1)
            end
        end
    end

    fη2 = zeros(nq)
    for i = 1:nx
        for j = 2:ny
            KitBase.flux_kfvs!(
                fη2,
                KitBase.maxwell_boltzmann_dual.(α[:, i, j-1]' * m)[:],
                KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:],
                vspace.u[:, 2],
                dt,
            )

            for k in axes(flux2, 1)
                flux2[k, i, j] = sum(m[k, :] .* (vspace.weights .* fη2))
            end
        end
    end

    # update
    for j = 2:ny-1
        for i = 2:nx-1
            for q = 1:1
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                    (SigmaS[i, j] * phi[q, i, j] - SigmaT[i, j] * phi[q, i, j]) * dt
            end

            for q = 2:ne
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                    (-SigmaT[i, j] * phi[q, i, j]) * dt
            end
        end
    end

    global t += dt
end

using Plots
contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], α[1, :, :])
contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :])


text:
"""
    Element2D(matrix::Array{Int, 2}, x::AbstractIntVar, y::AbstractIntVar, z::AbstractIntVar, trailer)

Element2D constraint, states that `matrix[x, y] == z`
"""
struct Element2D <: Constraint
    matrix::Array{Int, 2}
    x::AbstractIntVar
    y::AbstractIntVar
    z::AbstractIntVar
    n_rows_sup::Array{StateObject{Int}, 1}
    n_cols_sup::Array{StateObject{Int}, 1}
    low::StateObject{Int}
    up::StateObject{Int}
    xyz::Array{Tuple{Int, Int, Int}, 1}
    active::StateObject{Bool}

    function Element2D(matrix::Array{Int, 2}, x::AbstractIntVar, y::AbstractIntVar, z::AbstractIntVar, trailer)
        n, m = size(matrix)
        xyz = Array{Tuple{Int, Int, Int}, 1}()
        for i in 1:n
            for j in 1:m
                push!(xyz, (i, j, matrix[i, j]))
            end
        end
        sort!(xyz, by = xyz -> xyz[3])
        low = StateObject{Int}(1, trailer)
        up = StateObject{Int}(n*m, trailer)
        n_cols_sup = [StateObject{Int}(n, trailer) for i in 1:m]
        n_rows_sup = [StateObject{Int}(m, trailer) for i in 1:n]

        constraint = new(matrix, x, y, z, n_rows_sup, n_cols_sup, low, up, xyz, StateObject{Bool}(true, trailer))
        # we should prune below 1 and above (n resp. m) or add an offset
        for xi in [x, y, z]
            addOnDomainChange!(xi, constraint)
        end
        return constraint
    end
end

"""
    update_support_and_prune!(constraint::Element2D, xyz::Array{Tuple{Int, Int, Int}, 1}, lost_position::Int)

Supports are associated to x and y 's domains. Each possible value of x & y can be associated with n number of possible
values of z. If this number n equal 0, the value associated is pruned.

This function basically handle the increment, prune domains when necessary and retrieve the pruned values.
"""
function update_support_and_prune!(constraint::Element2D, xyz::Array{Tuple{Int, Int, Int}, 1}, lost_position::Int)
    setValue!(constraint.n_rows_sup[xyz[lost_position][1]], constraint.n_rows_sup[xyz[lost_position][1]].value - 1)
    setValue!(constraint.n_cols_sup[xyz[lost_position][2]], constraint.n_cols_sup[xyz[lost_position][2]].value - 1)
    prunedX, prunedY = [], []
    if constraint.n_rows_sup[xyz[lost_position][1]].value == 0
        prunedX = remove!(constraint.x.domain, xyz[lost_position][1])
    end
    if constraint.n_cols_sup[xyz[lost_position][2]].value == 0
        prunedY = remove!(constraint.y.domain, xyz[lost_position][2])
    end
    return prunedX, prunedY
end

"""
    propagate!(constraint::Element2D, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`Element2D` propagation function.

Supports (n_cols_sup & n_rows_sup) are the number of possible z values that
are in a row or column. Once a support reaches 0, the value corresponding to the index cannot be taken by
the corresponding variable.
x = [1, 2, 3]
y = [1, 2, 3]
z = [6, 7]
          | 1 | 2 | 3 | n_rows_sup
        --------------
        1 | 3   4   2    0
        2 | 6   7   1    2
        3 | 7   5   1    1

n_cols_sup  2   1   0

Here, the value 3 can be pruned from y and the value 1 can be pruned from x.

Only the lowest and the highest values of z are pruned when they are not in the matrix.
(it's be for performance reason, this implemention is inspired by miniCP)

It'd be interesting to try to benchmark with a propagation function which also pruned
non extreme values to make sure it's not better.

"""
function propagate!(constraint::Element2D, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    # get useful variables
    low = constraint.low.value
    up = constraint.up.value
    z_min = minimum(constraint.z.domain)
    z_max = maximum(constraint.z.domain)
    xyz = constraint.xyz
    n, m = size(constraint.matrix)

    # make sure x, y are valid indexes of matrix
    prunedX = vcat(removeBelow!(constraint.x.domain, 1), removeAbove!(constraint.x.domain, n))
    prunedY = vcat(removeBelow!(constraint.y.domain, 1), removeAbove!(constraint.y.domain, m))
    if !isempty(prunedX)
        addToPrunedDomains!(prunedDomains, constraint.x, prunedX)
        triggerDomainChange!(toPropagate, constraint.x)
    end
    if !isempty(prunedY)
        addToPrunedDomains!(prunedDomains, constraint.y, prunedY)
        triggerDomainChange!(toPropagate, constraint.y)
    end

    # update the supports of the cols and the rows and prune X or Y if some supports equal 0
    # lower bound
    while xyz[low][3] < z_min || !(xyz[low][1] in constraint.x.domain) || !(xyz[low][2] in constraint.y.domain)
        prunedX, prunedY = update_support_and_prune!(constraint, xyz, low)
        if !isempty(prunedX)
            addToPrunedDomains!(prunedDomains, constraint.x, prunedX)
            triggerDomainChange!(toPropagate, constraint.x)
        end
        if !isempty(prunedY)
            addToPrunedDomains!(prunedDomains, constraint.y, prunedY)
            triggerDomainChange!(toPropagate, constraint.y)
        end
        low += 1
        # @assert low <= up
        if low > up
            return false
        end
    end
    # upper bound
    while xyz[up][3] > z_max || !(xyz[up][1] in constraint.x.domain) || !(xyz[up][2] in constraint.y.domain)
        prunedX, prunedY = update_support_and_prune!(constraint, xyz, up)
        if !isempty(prunedX)
            addToPrunedDomains!(prunedDomains, constraint.x, prunedX)
            triggerDomainChange!(toPropagate, constraint.x)
        end
        if !isempty(prunedY)
            addToPrunedDomains!(prunedDomains, constraint.y, prunedY)
            triggerDomainChange!(toPropagate, constraint.y)
        end
        up -= 1
        # @assert low <= up
        if low > up
            return false
        end
    end

    # try to prune lower or upper values of z's domain
    prunedZ = vcat(removeBelow!(constraint.z.domain, xyz[low][3]), removeAbove!(constraint.z.domain, xyz[up][3]))
    if !isempty(prunedZ)
        addToPrunedDomains!(prunedDomains, constraint.z, prunedZ)
        triggerDomainChange!(toPropagate, constraint.z)
    end

    # update useful variables
    setValue!(constraint.low, low)
    setValue!(constraint.up, up)

    # deactivate the constraint if necessary
    if isbound(constraint.z)
        zv = assignedValue(constraint.z)
        if all(zv == constraint.matrix[vx, vy] for vx in constraint.x.domain for vy in constraint.y.domain)
            setValue!(constraint.active, false)
        end
    end

    # check feasibility
    return !isempty(constraint.x.domain) && !isempty(constraint.y.domain) && !isempty(constraint.z.domain)
end

variablesArray(constraint::Element2D) = [constraint.x, constraint.y, constraint.z]

function Base.show(io::IO, ::MIME"text/plain", con::Element2D)
    println(io, "Element2D constraint: maxtrix[$(con.x.id), $(con.y.id)] == $(con.z.id)", ", active = ", con.active)
    println(io, "   matrix = ", con.matrix)
    println(io, "   ", con.x)
    println(io, "   ", con.y)
    print(io, "   ", con.z)
end

function Base.show(io::IO, con::Element2D)
    print(io, typeof(con), ": maxtrix[$(con.x.id), $(con.y.id)] == $(con.z.id)")
end


text:
using AutoHashEquals

# The type parameters specify the (user-defined)
#   S: state
#   D: variable domain
#   V: objective value

"""
    Arc{S,D,V}

An arc in the decision diagram, representing a state transition.

It points to the original/previous state and also stores the decision made
(variable assignment) as well as the contribution to the objective function.
"""
@auto_hash_equals struct Arc{S,D,V}
    tail::S
    decision::D
    value::V
end

function Base.isapprox(x::Arc{S,D,V}, y::Arc{S,D,V}) where {S,D,V}
    return x.tail == y.tail && x.decision == y.decision && x.value ≈ y.value
end

"""
    Node{S,D,V}

Meta-data for a node in the decision diagram.

Stores the distance from the root node on the longest path so far, the
ingoing arc on such a path (but no other ingoing arcs) and a flag to specify
whether the state is *exact*, as opposed to *relaxed*.
"""
@auto_hash_equals struct Node{S,D,V}
    distance::V
    inarc::Union{Arc{S,D,V},Nothing}
    exact::Bool

    function Node{S,D,V}(distance, inarc=nothing, exact=true) where {S,D,V}
        new(distance, inarc, exact)
    end
end

function Base.isapprox(x::Node{S,D,V}, y::Node{S,D,V}) where {S,D,V}
    return x.distance ≈ y.distance && x.inarc ≈ y.inarc && x.exact == y.exact
end

"""
    Layer{S,D,V}

A layer of nodes in the decision diagram.

Represented by mapping from (user-defined) states to the Node meta-data. Also
has a flag `exact` to indicate whether all states are represented exactly
(neither restricted nor relaxed).
"""
@auto_hash_equals struct Layer{S,D,V}
    nodes::Dict{S,Node{S,D,V}}
    exact::Bool

    function Layer{S,D,V}(nodes=Dict(), exact=true) where {S,D,V}
        return new(nodes, exact)
    end
end

function Base.iterate(layer::Layer{S,D,V}) where {S,D,V}
    return iterate(layer.nodes)
end

function Base.iterate(layer::Layer{S,D,V}, state) where {S,D,V}
    return iterate(layer.nodes, state)
end

function Base.length(layer::Layer{S,D,V}) where {S,D,V}
    return length(layer.nodes)
end

function Base.haskey(layer::Layer{S,D,V}, state::S) where {S,D,V}
    return haskey(layer.nodes, state)
end

function Base.getindex(layer::Layer{S,D,V}, state::S) where {S,D,V}
    return getindex(layer.nodes, state)
end

function Base.setindex!(
    layer::Layer{S,D,V}, node::Node{S,D,V}, state::S
) where {S,D,V}
    return setindex!(layer.nodes, node, state)
end

"""
    Diagram{S,D,V}

A (multi-valued) decision diagram.

It's a directed acyclic graph where the nodes represent (feasible) states and
the arcs transitions triggered by decision variable assignments. Decisions are
made sequentially and arcs only connect consecutive layers. The initial layer
contains the single, given root node. All nodes in the final layer are merged to
a single terminal node.

As the variable order can be defined dynamically, the variable indices are also
stored. Note that the constructed diagram will have N+1 layers for N variables.

There is also a property `partial_sol` containing indices of variables that are
already assigned outside this diagram (in the context of branch-and-bound).
"""
@auto_hash_equals struct Diagram{S,D,V}
    partial_sol::Vector{Int}
    layers::Vector{Layer{S,D,V}}
    variables::Vector{Int}
end

function Diagram(initial::Layer{S,D,V}) where {S,D,V}
    return Diagram{S,D,V}([], [initial], [])
end

function Diagram(instance)
    state = initial_state(instance)
    S = typeof(state)
    D = domain_type(instance)
    V = value_type(instance)
    node = Node{S,D,V}(zero(V))
    root = Layer{S,D,V}(Dict(state => node))
    return Diagram(root)
end

"""
    Solution{D,V}

A feasible solution, with decisions for all variables (in order) and the
objective values.
"""
@auto_hash_equals struct Solution{D,V}
    decisions::Vector{D}  # for all variables, order 1:n
    objective::V
end

"""
    Subproblem{D,V}

A subproblem in the context of branch-and-bound, as defined by an exact node in
the diagram. It's represented by the partial solution given by a longest path
from the root that node and the current state.
"""
@auto_hash_equals struct Subproblem{S,D,V}
    # partial solution (assigned so far, in given order)
    variables::Vector{Int}
    decisions::Vector{D}
    distance::V

    # state (to complete solution)
    state::S
end


text:
function parse_line_segments(path::String)::Array{Int}
    line_segments = Array{Int}(undef, 1, 4)
    open(path) do file
        lines = readlines(file)
        for line in lines
            coordinates = split(line, " -> ")
            segment = [0 0 0 0]
            i = 1
            for coordinate in coordinates
                x = parse(Int, split(coordinate, ",")[1])
                y = parse(Int, split(coordinate, ",")[2])
                segment[i] = x
                segment[i+1] = y
                i += 2
            end
            line_segments = vcat(line_segments, segment)
        end
    end
    return line_segments[2:end, :]
end

function draw_map(line_segments)
    min, max = minimum(line_segments), maximum(line_segments)
    field = zeros(Int, max-min+1, max-min+1)
    line_segments = line_segments .+ (-min + 1)
    for i in 1:size(line_segments)[1]
        segment = line_segments[i, :]
        draw_line(field, segment)
    end
    return field
end

function draw_line(field, line_segment)
    x1, y1 = line_segment[1], line_segment[2]
    x2, y2 = line_segment[3], line_segment[4]
    x_step = x1 < x2 ? 1 : -1
    y_step = y1 < y2 ? 1 : -1
    if x1 == x2
        field[y1:y_step:y2, x1] = field[y1:y_step:y2, x1] .+ 1
    elseif y1 == y2
        field[y1, x1:x_step:x2] = field[y1, x1:x_step:x2] .+ 1
    elseif abs(x1-x2) == abs(y1-y2)
        for (x, y) in zip(x1:x_step:x2, y1:y_step:y2)
            field[y, x] += 1
        end
    end
end

function find_dangerous_areas(map)
    overlaps = map .> 1
    return sum(overlaps)
end

function main()
    path = "./inputs/puzzle.txt"
    segments = parse_line_segments(path)
    map = draw_map(segments)
    display(map)
    overlaps = find_dangerous_areas(map)
    println(overlaps)
end

main()