text:
module PhaseDiagrams

# This contains routines to read phase diagrams from disk

using ..Units
using Unitful
using DelimitedFiles, Interpolations
import Base.show
using GeoParams: AbstractMaterialParam, AbstractPhaseDiagramsStruct
import GeoParams.PerpleX_LaMEM_Diagram

export  PhaseDiagram_LookupTable, PerpleX_LaMEM_Diagram,
        ComputeDensity

"""
    Contains data of a Phase Diagram that is regularly spaced in P & T
"""
struct PhaseDiagram_LookupTable <: AbstractPhaseDiagramsStruct
    Type        ::  String
    HeaderText  ::  Any
    Name        ::  String
    rockRho     ::  Any
    meltRho     ::  Any
    meltFrac    ::  Any
    Rho         ::  Any
    rockVp      ::  Any
    rockVs      ::  Any
    rockVpVs    ::  Any
    meltVp      ::  Any
    meltVs      ::  Any
    meltVpVs    ::  Any
    Vp          ::  Any
    Vs          ::  Any
    VpVs        ::  Any
    cpxFrac     ::  Any
end

"""
    PD_Data = PerpleX_LaMEM_Diagram(fname::String; CharDim = nothing)

Reads a precomputed phase diagram in the `LaMEM/Perple_X` format (which is a phase diagram computed using `Perple_X`, but formatted in a manner that is readable using LaMEM).
The data is stored in the `PhaseDiagram_LookupTable` structure.

If the `CharDim` object is specified, the values of all diagrams will be non-dimensionalized.

# Example

```julia
julia> PD_Data = PerpleX_LaMEM_Diagram("./test_data/Peridotite.in")
Perple_X/LaMEM Phase Diagram Lookup Table:
                      File    :   ./test_data/Peridotite.in
                      T       :   293.0 - 1573.000039
                      P       :   1.0e7 - 2.9999999944e9
                      fields  :   :meltRho, :meltRho, :meltFrac, :rockRho, :Rho, :rockVp
                                  :rockVs, :rockVpVs, :meltVp, :meltVs, :meltVpVs
                                  :Vp, :Vs, :VpVs
```
Once imported, the properties on the diagram can be interpolated in a simple manner:
```
julia> PD_Data.Rho(1500,1e7)
3042.836820256982
```
This also works for vectors or arrays:
```julia
julia> T = [1500 1800; 1233 1300]
julia> P = [1e8 1e9; 1e7 1e7]
julia> rho = PD_Data.Rho.(T,P)
```
(Note the dot `.` in front of the bracket while evaluating arrays).

The fields that are available depend on what is listed in the diagram file.
The units of the fields are automatically evaluated, and employed to non-dimensionalize the parameters if `CharDim` is specified.

# Algorithm

Internally, we employ linear interpolation, as provided by the [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl) package.
Values outside the range of the diagram are set to the boundary of the diagram. The interpolation object is directly encoded in the `PhaseDiagram_LookupTable`` object.

"""
function PerpleX_LaMEM_Diagram(fname::String; CharDim = nothing)

    # Read header:
    #  the first 50 lines are comments (freely useable), followed by data
    n           = 55
    header      = open(readlines, `head -n $(n) $(fname)`)
    header_text = header[1:49]

    # Parse the names of the collumns in the data file
    fields      = split(header[49],"[")[2:end]      # this line should contain the names and units of the collumns in the file

    # Throw error message if fields are not specified
    if length(fields)==0
        error("Line 49 in the phase diagram file should contain the names of the collums")
    end

    fields_keys  = Vector{Union{Missing, Symbol}}(missing, length(fields))
    fields_units = Vector{Union{Missing, Any}   }(missing, length(fields))
    for i=1:length(fields)
        # error catch
        if length(split(fields[i],","))!= 3
            error("Error parsing collumn: $(fields[i]). Should have 3 fields (name, symbol, units)")
        end

        fields_keys[i]  = Symbol(split(fields[i],",")[1])
        unit_string     = split(split(fields[i],",")[3],"]")[1]
        fields_units[i] = uparse(unit_string)
    end

    # Determine the range
    T0      =   parse(Float64,header[50])*u"K"      # in K
    dT      =   parse(Float64,header[51])*u"K"
    numT    =   parse(Int64,  header[52])

    P0      =   parse(Float64,header[53])*u"bar"    # in bar (will be convert to Pa later)
    dP      =   parse(Float64,header[54])*u"bar"
    numP    =   parse(Int64,  header[55])

    Tvec    =   T0:dT:(T0+dT*(numT-1))              # 1D vector
    Pvec    =   P0:dP:(P0+dP*(numP-1))

    # In the LaMEM/Perple_X file format, the first 50 lines are comments
    data    =   readdlm(fname,skipstart=55,header=false);        # read numerical data

    # Shape of 2D arrays:
    siz      =   (numP, numT)

    # Initialize fields in the order they are defined in the PhaseDiagram_LookupTable structure
    Struct_Fieldnames = fieldnames(PhaseDiagram_LookupTable)[4:end] # fieldnames from structure

    # Process all fields that are present in the phase diagram (and non-dimensionalize if requested)
    Struct_Fields     =  Vector{Union{Nothing, Interpolations.Extrapolation}}(nothing, length(Struct_Fieldnames))

    # Loop through all fields and create an interpolation object for each (which is optionally non-dimensionalized)
    for (i, field) in enumerate(fields_keys)
        ind = findall(Struct_Fieldnames.==field)
        if length(ind)>0
            Struct_Fields[ind[1]]   = CreateInterpolationObject_PhaseDiagram(data[:,i], Tvec, Pvec, siz, fields_units[i], CharDim)
        end
    end

    # Some fields have melt and solid part; we can reconstruct the total part as an arithmetic average:
    Struct_Fields = ComputeTotalField_withMeltFraction(:Rho, :meltRho, :rockRho, :meltFrac, Struct_Fields, Struct_Fieldnames)
    Struct_Fields = ComputeTotalField_withMeltFraction(:Vp,  :meltVp,  :rockVp,  :meltFrac, Struct_Fields, Struct_Fieldnames)
    Struct_Fields = ComputeTotalField_withMeltFraction(:Vs,  :meltVs,  :rockVs,  :meltFrac, Struct_Fields, Struct_Fieldnames)
    Struct_Fields = ComputeTotalField_withMeltFraction(:VpVs,:meltVpVs,:rockVpVs,:meltFrac, Struct_Fields, Struct_Fieldnames)

    # Store in phase diagram structure
    PD_data = PhaseDiagram_LookupTable("Perple_X/LaMEM", header_text, fname, Struct_Fields... )

    return PD_data
end

# Print info
function show(io::IO, d::PhaseDiagram_LookupTable)
    T = d.rockRho.itp.knots[1]
    P = d.rockRho.itp.knots[2]

    println(io, "$(d.Type) Phase Diagram Lookup Table: ")
    println(io, "                      File    :   $(d.Name)")
    println(io, "                      T       :   $(minimum(T)) - $(maximum(T))")
    println(io, "                      P       :   $(minimum(P)) - $(maximum(P))")

    lst = fieldnames(typeof(d))
    str=""
    for i=4:length(lst)

        if !isnothing(getfield(d,lst[i]))
            if i==4
                str = "                      fields  :   :$(lst[i]),"
            end
            str = str*" :$(lst[i]),"
            if mod( (i-3),5)==0
                str = str[1:end-1]
                println(io, str)
                str = "                                 "
            end
        end
    end
    str = str[1:end-1]
    println(io, str)

end

# Internal routine that creates an interpolation object from a collumn of the data
function CreateInterpolationObject_PhaseDiagram(data_vec::Vector{Float64}, Tvec, Pvec, siz::Tuple{Int64, Int64}, units, CharDim)

    data_units  =   reshape(data_vec, siz)*units;      # Create 2D array

    # Convert to Pa & K
    Pvec_Pa     =   Float64.(uconvert.(u"Pa", Pvec))
    Tvec_K      =   Float64.(uconvert.(u"K",  Tvec))

    # Optional: nondimensionalize values as well as T and P
    if CharDim==nothing
        Pvec = ustrip.(Pvec_Pa);
        Tvec = ustrip.(Tvec_K);
        data = ustrip.(data_units);
    else
        Pvec = Nondimensionalize(Pvec_Pa, CharDim)
        Tvec = Nondimensionalize(Tvec_K, CharDim)
        data = Nondimensionalize(data_units, CharDim)
    end

    # Create interpolation object
    intp_data    =   LinearInterpolation((Tvec, Pvec), data,      extrapolation_bc = Flat());

    return intp_data
end

# Internal routine, which combines solid & melt properties with melt fraction
function ComputeTotalField_withMeltFraction(totalData::Symbol, meltData::Symbol, solidData::Symbol, meltFrac::Symbol, Struct_Fields, Struct_Fieldnames)

    ind_totalData = findall(Struct_Fieldnames .== totalData)
    ind_meltData  = findall(Struct_Fieldnames .== meltData )
    ind_solidData = findall(Struct_Fieldnames .== solidData)
    ind_meltFrac  = findall(Struct_Fieldnames .== meltFrac )

    # extract arrays from interpolation objects
    if  (length(ind_totalData)>0) & (length(ind_meltData)>0) &
        (length(ind_solidData)>0) & (length(ind_meltFrac)>0) &
        !(isnothing(Struct_Fields[ind_meltFrac ])) &
        !(isnothing(Struct_Fields[ind_solidData])) &
        !(isnothing(Struct_Fields[ind_meltData ]))

        if Struct_Fields[ind_meltFrac[1]]!=nothing
            ϕ = Struct_Fields[ind_meltFrac[1] ].itp.coefs      # melt fraction
            S = Struct_Fields[ind_solidData[1]].itp.coefs      # solid property
            M = Struct_Fields[ind_meltData[1] ].itp.coefs      # melt property

            # Compute average
            Result      =   (1.0 .- ϕ).*S .+ ϕ.*M

            # Create interpolation object of average
            Tvec        =   Struct_Fields[ind_meltFrac[1] ].itp.knots[1]       # temperature data
            Pvec        =   Struct_Fields[ind_meltFrac[1] ].itp.knots[2]       # pressure data
            intp_Result =   LinearInterpolation((Tvec, Pvec), Result,  extrapolation_bc = Flat());

            # assign
            Struct_Fields[ind_totalData[1]] = intp_Result
        else
            # use solid fraction
            S = Struct_Fields[ind_solidData[1]]      # solid property

            # assign
            Struct_Fields[ind_totalData[1]] = S

        end

    end

    return Struct_Fields
end

end


text:
## TODO Update Refactor

using Plots
using PortfolioOpt
using PortfolioOpt.TestUtils

include("./examples/literature_examples/example_utils.jl")

############ Read Prices #############
Prices = get_test_data()
numD, numA = size(Prices) # A: Assets    D: Days

############# Calculating returns #####################
returns_series = percentchange(Prices)
returns = values(returns_series)
# risk free asset
rf = fill(3.2e-6, numD) # readcsv(".\\rf.csv")

############ Efficient frontier ###################
t = 100
k_back = 60
Σ, r̄ = mean_variance(returns[(t - k_back):(t - 1), :])

range = -0.01:0.0005:0.025;
len = size(range, 1);
E = zeros(len);
σ = zeros(len);
E_quad = zeros(len);
σ_quad = zeros(len);
E_quad_Rf = zeros(len);
σ_quad_Rf = zeros(len);

x_sharpe = max_sharpe(Σ, r̄, rf[t])
E_sharpe = r̄'x_sharpe
σ_sharpe = sqrt(x_sharpe'Σ * x_sharpe)
iter = 0
for R in range
    global iter += 1
    # analytical
    x = mean_variance_noRf_analytical(Σ, r̄, R)
    # quadratic optimization no risk free asset
    model, w = base_model(numA)
    po_minvar_limitmean_noRf!(model, w, Σ, r̄, R)
    x_quad, obj, r = compute_solution(model, w)
    # quadratic optimization with risk free asset
    model, w = base_model(numA)
    po_minvar_limitmean_Rf!(model, w, Σ, r̄, R, rf[t], 1)
    x_quad_Rf, obj, r = compute_solution(model, w)

    E[iter] = sum(x'r̄)
    σ[iter] = sqrt((x'Σ * x)[1])
    E_quad[iter] = r[1]
    σ_quad[iter] = sqrt(x_quad'Σ * x_quad)
    E_quad_Rf[iter] = R
    σ_quad_Rf[iter] = sqrt(x_quad_Rf'Σ * x_quad_Rf)
end

# plot frontier
plt = plot(
    σ_quad,
    E_quad;
    title="Efficient Frontier Mean-Variance",
    xlabel="σ",
    ylabel="r",
    label="Quadratic",
    legend=:outertopright,
);
plot!(plt, [0; σ_sharpe], [rf[t]; E_sharpe]; label="Sharpe");
plot!(plt, σ, E; label="Analytical");
plot!(plt, σ_quad_Rf, E_quad_Rf; label="Quadratic with RF");
plt

ext: jl


text:
function doctest(lib::PolyhedraLibrary)
    A = [1 1;1 -1;-1 0]
    b = [1,0,0]
    V = [1//2 1//2; 0 1; 0 0]
    hr = hrep(A, b)
    p = polyhedron(hr, lib)
    @test p ⊆ HalfSpace([1, -1], 0)
    @test !(p ⊆ HalfSpace([1, 1], 0))
    @test !(p ⊆ HyperPlane([-1, 0], 0))
    inequality_fulltest(p, A, b, BitSet())
    generator_fulltest(p, V)
    @test hrepiscomputed(p)
    @test vrepiscomputed(p) # computation triggered in generator_fulltest
    @test [1//4, 1//2] in p
    @test ininterior([1//4, 1//2], p)
    @test inrelativeinterior([1//4, 1//2], p)
    @test !([1, 1//2] in p)
    @test !ininterior([1, 1//2], p)
    @test !inrelativeinterior([1, 1//2], p)
    @test [0, 1] in p
    @test !ininterior([0, 1], p)
    @test !inrelativeinterior([0, 1], p)
    @test !(Ray([0, 1]) in p)
    @test !ininterior(Ray([0, 1]), p)
    @test !inrelativeinterior(Ray([0, 1]), p)
    Ap = reshape([1, -1], 2, 1)
    bp = [1, 0]
    inequality_fulltest(eliminate(p, [1]), Ap, bp, BitSet())
    inequality_fulltest(project(p, [2]), Ap, bp, BitSet())
    inequality_fulltest(eliminate(p, [1], ProjectGenerators()), Ap, bp, BitSet())
    inequality_fulltest(project(p, [2], ProjectGenerators()), Ap, bp, BitSet())
    if supportselimination(p, BlockElimination())
        inequality_fulltest(eliminate(p, [1], BlockElimination()), Ap, bp, BitSet())
        inequality_fulltest(project(p, [2], BlockElimination()), Ap, bp, BitSet())
    end
    @test_throws DimensionMismatch project(p, ones(1, 1))
    @test_throws DimensionMismatch project(p, ones(2, 3))
    inequality_fulltest(project(p, reshape([0, 1], 2, 1)), Ap, bp, BitSet())
    inequality_fulltest(fixandeliminate(p, 2, eltype(p) <: AbstractFloat ? 1/4 : 1//4), reshape([1, -1], 2, 1), [1/4, 0], BitSet())
    # tries to create CDDPolyhedron of BigFloat, needs to let CDD decide of the promotion
    #inequality_fulltest(fixandeliminate(p, 2, 1/4, reshape([1, -1], 2, 1), [1/4, 0], BitSet())
    inequality_fulltest(translate(p, -[0, 1]), [1 1; 1 -1;-1 0], [0, 1, 0], BitSet())
    # Testing 2D plotting
    RecipesBase.apply_recipe(Dict{Symbol, Any}(), p)[1].args == ([0, 0.5, 0.5, 0, 0], [1, 0.5, 0.5, 0, 1])
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 TopologyParameters <: SwaggerModel
    targetResourceGroupName::Any # spec type: Union{ Nothing, String } # spec name: targetResourceGroupName
    targetVirtualNetwork::Any # spec type: Union{ Nothing, SubResource } # spec name: targetVirtualNetwork
    targetSubnet::Any # spec type: Union{ Nothing, SubResource } # spec name: targetSubnet

    function TopologyParameters(;targetResourceGroupName=nothing, targetVirtualNetwork=nothing, targetSubnet=nothing)
        o = new()
        validate_property(TopologyParameters, Symbol("targetResourceGroupName"), targetResourceGroupName)
        setfield!(o, Symbol("targetResourceGroupName"), targetResourceGroupName)
        validate_property(TopologyParameters, Symbol("targetVirtualNetwork"), targetVirtualNetwork)
        setfield!(o, Symbol("targetVirtualNetwork"), targetVirtualNetwork)
        validate_property(TopologyParameters, Symbol("targetSubnet"), targetSubnet)
        setfield!(o, Symbol("targetSubnet"), targetSubnet)
        o
    end
end # type TopologyParameters

const _property_map_TopologyParameters = Dict{Symbol,Symbol}(Symbol("targetResourceGroupName")=>Symbol("targetResourceGroupName"), Symbol("targetVirtualNetwork")=>Symbol("targetVirtualNetwork"), Symbol("targetSubnet")=>Symbol("targetSubnet"))
const _property_types_TopologyParameters = Dict{Symbol,String}(Symbol("targetResourceGroupName")=>"String", Symbol("targetVirtualNetwork")=>"SubResource", Symbol("targetSubnet")=>"SubResource")
Base.propertynames(::Type{ TopologyParameters }) = collect(keys(_property_map_TopologyParameters))
Swagger.property_type(::Type{ TopologyParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TopologyParameters[name]))}
Swagger.field_name(::Type{ TopologyParameters }, property_name::Symbol) =  _property_map_TopologyParameters[property_name]

function check_required(o::TopologyParameters)
    true
end

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


text:
# FASTA Reader
# ============

"""
    FASTAReader(input::IO; index=nothing)
    FASTAReader{S}(input::IO; index=nothing)

Create a data reader of the FASTA file format.

When type parameter `S` is specified, the reader reads sequences in that type;
otherwise the reader tries to infer the sequence type based on the frequencies
of characters from the input.

# Arguments
* `input`: data source
* `index=nothing`: filepath to a random access index (currently *fai* is supported)
"""
type FASTAReader{S<:Sequence} <: Bio.IO.AbstractReader
    state::Ragel.State
    seqbuf::BufferedOutputStream{BufferedStreams.EmptyStream}
    index::Nullable{FASTAIndex}

    function FASTAReader(input::BufferedInputStream, index)
        return new(Ragel.State(fasta_machine.start_state, input),
                   BufferedOutputStream(), index)
    end
end

function Bio.IO.stream(reader::FASTAReader)
    return reader.state.stream
end

function Base.eltype{S}(::Type{FASTAReader{S}})
    return FASTASeqRecord{S}
end

function FASTAReader(input::IO; index=nothing)
    # alphabet type will be inferred from data
    return FASTAReader{BioSequence}(input, index=index)
end

function (::Type{FASTAReader{S}}){S<:Sequence}(input::IO; index=nothing)
    if isa(index, AbstractString)
        index = FASTAIndex(index)
    else
        if index != nothing
            error("unrecognizable index argument")
        end
    end
    return FASTAReader{S}(BufferedInputStream(input), index)
end

function Base.getindex(reader::FASTAReader, name::AbstractString)
    if isnull(reader.index)
        error("no index")
    end
    seekrecord(reader.state.stream, get(reader.index), name)
    reader.state.cs = fasta_machine.start_state
    return read(reader)
end

# Predict sequence type based on character frequencies in `seq[start:stop]`.
function predict(seq::Vector{UInt8}, start, stop)
    # count characters
    a = c = g = t = u = n = alpha = 0
    for i in start:stop
        @inbounds x = seq[i]
        if x == 0x41 || x == 0x61
            a += 1
        elseif x == 0x43 || x == 0x63
            c += 1
        elseif x == 0x47 || x == 0x67
            g += 1
        elseif x == 0x54 || x == 0x74
            t += 1
        elseif x == 0x55 || x == 0x75
            u += 1
        elseif x == 0x4e || x == 0x6e
            n += 1
        end
        if 0x41 ≤ x ≤ 0x5a || 0x61 ≤ x ≤ 0x7a
            alpha += 1
            if alpha ≥ 300 && t + u > 0 && a + c + g + t + u + n == alpha
                # pretty sure that the sequence is either DNA or RNA
                break
            end
        end
    end

    # the threshold (= 0.95) is somewhat arbitrary
    if (a + c + g + t + u + n) / alpha > 0.95
        if t ≥ u
            return DNASequence
        else
            return RNASequence
        end
    else
        return AminoAcidSequence
    end
end

info("compiling old FASTA")
const fasta_machine = (function ()
    re = Automa.RegExp

    lf          = re"\n"
    newline     = re"\r?" * lf
    hspace      = re"[ \t\v]"
    whitespace  = re.space() | newline
    identifier  = re.rep1(re.any() \ re.space())
    description = (re.any() \ hspace) * re"[^\r\n]*"
    letters     = re.rep1(re.any() \ re.space() \ re">")
    sequence    = re.cat(re.rep(whitespace), re.opt(letters), re.rep(re.rep1(whitespace) * letters))
    record      = re.cat('>', identifier, re.opt(re.rep1(hspace) * description), newline, sequence, re.rep(whitespace))
    fasta       = re.rep(whitespace) * re.rep(record)

    lf.actions[:enter]          = [:count_line]
    identifier.actions[:enter]  = [:mark]
    identifier.actions[:exit]   = [:identifier]
    description.actions[:enter] = [:mark]
    description.actions[:exit]  = [:description]
    letters.actions[:enter]     = [:mark]
    letters.actions[:exit]      = [:letters]
    record.actions[:exit]       = [:record]

    return Automa.compile(fasta)
end)()

@inline function anchor!(stream, p)
    stream.anchor = p
end

@inline function upanchor!(stream)
    @assert stream.anchor != 0 "upanchor! called with no anchor set"
    anchor = stream.anchor
    stream.anchor = 0
    return anchor
end

const fasta_actions = Dict(
    :count_line  => :(linenum += 1),
    :mark        => :(anchor!(stream, p)),
    :identifier  => :(copy!(output.name, stream.buffer, upanchor!(stream), p - 1)),
    :description => :(copy!(output.metadata.description, stream.buffer, upanchor!(stream), p - 1)),
    :letters     => :(append!(reader.seqbuf, stream.buffer, upanchor!(stream), p - 1)),
    :record      => :(found_record = true; @escape))

function Base.read!(reader::FASTAReader, output::FASTASeqRecord)
    return _read!(reader, reader.state, output)
end

@eval function _read!(reader::FASTAReader, state::Ragel.State, output::FASTASeqRecord)
    cs = state.cs
    linenum = state.linenum
    stream = state.stream
    data = stream.buffer
    p = stream.position
    p_end = stream.available
    p_eof = -1
    found_record = false

    while true
        $(Automa.generate_exec_code(fasta_machine, actions=fasta_actions, code=:goto, check=false))

        state.cs = cs
        state.finished = cs == 0
        state.linenum = linenum
        stream.position = p

        if cs < 0
            error("FASTA file format error on line ", linenum)
        elseif found_record
            if seqtype(typeof(output)) == ReferenceSequence
                output.seq = ReferenceSequence(reader.seqbuf.buffer, 1, length(reader.seqbuf))
            elseif seqtype(typeof(output)) == BioSequence
                ET = predict(reader.seqbuf.buffer, 1, length(reader.seqbuf))
                if ET == typeof(output.seq)
                    resize!(output.seq, length(reader.seqbuf))
                    encode_copy!(output.seq, 1, reader.seqbuf.buffer, 1, length(reader.seqbuf))
                else
                    output.seq = ET(reader.seqbuf.buffer, 1, length(reader.seqbuf))
                end
            else
                resize!(output.seq, length(reader.seqbuf))
                encode_copy!(output.seq, 1, reader.seqbuf.buffer, 1, length(reader.seqbuf))
            end
            empty!(reader.seqbuf)
            break
        elseif cs == 0
            throw(EOFError())
        elseif p > p_eof ≥ 0
            error("incomplete FASTA input on line ", linenum)
        else
            hits_eof = BufferedStreams.fillbuffer!(stream) == 0
            p = stream.position
            p_end = stream.available
            if hits_eof
                p_eof = p_end
            end
        end
    end

    return output
end


text:
################################################################################

module HiddenMarkovModelReaders

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

# dependencies
using StructArrays
using Parameters

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

# hidden Markov model functions
export setup, process!

# hidden Markov model object & parameters
export HMM, HMMParams

# distance functions
export euclDist, bhattDist

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

# hidden markov model object
include("struct.jl")

# methods
include("reset.jl")
include("feed.jl")
include("backTrace.jl")
include("sort.jl")

# distances
include("distance.jl")

# setup
include("setup.jl")

# process
include("process.jl")

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

end

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

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_forks_repo_name: DanielRivasMD/HiddenMarkovModelReaders


text:
# """
#     An Uniform playing strategy for Bandits: Can be used as benchmark
# """

struct  UniformStrategy <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    function UniformStrategy( noOfArms::Int64 )
        new( noOfArms,
             0,
             0
        )
    end
end

function get_arm_index( agent::UniformStrategy )
    agent.lastPlayedArm = rand( 1:agent.noOfArms )
    return agent.lastPlayedArm
end

function update_reward!( agent::UniformStrategy, r::Real )
    agent.noOfSteps += 1
    # Do nothing

    nothing
end

function reset!( agent::UniformStrategy )
    agent.lastPlayedArm = 0
    agent.noOfSteps = 0
    # Do nothing

    nothing
end

function info_str( agent::UniformStrategy, latex::Bool )
    return @sprintf( "Uniform Strategy" )
end


text:
struct LagrangianEntry{E} <: Entry
    index::Int
    e::E
end
struct LagrangianHessian{F}
    f::F
end
@inline (laghess::LagrangianHessian)(z,x,p,l,s) = _eval_laghess(laghess.f,z,x,p,l,s)
@inline function _eval_laghess(f::Field1{E,I},z,x,p,l,s) where {E,I}
    _eval_laghess(inner(f),z,x,p,l,s)
    @simd for i in eachindex(f.es)
        f.es[i](z,x,p,l)
    end
end
@inline _eval_laghess(h::H,z,x,p,l,s) where H <: Hessian = h(z,x,p,s)
@inline (lagentry::LagrangianEntry{E})(z,x,p,l) where E = lagentry.e(z,x,p,l[index(lagentry)])

LagrangianHessian(obj,grad,con,jac,indexer=nothing) = LagrangianHessian(_LagrangianHessian(obj,grad,con,jac,indexer))

function _LagrangianHessian(obj,grad,con::Field1{E1,I1},jac::Field1{E2,I2},indexer = nothing) where {E1,E2,I1,I2}
    Field1(_LagrangianHessian(obj,grad,inner(con),inner(jac),indexer),
           [LagrangianEntry(index(con.es[i]),Hessian(con.es[i].e,jac.es[i].e,indexer)) for i in eachindex(con.es)])
end

function _LagrangianHessian(obj,grad,con::Field1{E1,Nothing},jac::Field1{E2,Nothing},indexer = nothing) where {E1,E2}
    Field1(Hessian(obj,grad,indexer),[LagrangianEntry(index(con.es[i]),Hessian(con.es[i].e,jac.es[i].e,indexer)) for i in eachindex(con.es)])
end

function _LagrangianHessian(obj,grad,con::FieldNull,jac::FieldNull,indexer = nothing) where {E1,E2}
    Hessian(obj,grad,indexer)
end


text:
"""
    IndexedGraph

A representation of an adjaacency graph optimized for redistricting applications.
Node-level population is stored in a special array; node neighbors and mappings
between node pairs and edge indices are also broken out into special arrays for
ease of access. Auxiliary metadata is stored as a dictionary in the `attributes`
column.

Importantly, this structure does _not_ store plan-level information, such as
a mapping between nodes and districts. This higher-level information should be
stored in a separate `Plan` structure.
"""
struct IndexedGraph
    n_nodes::Int
    n_edges::Int
    neighbors_per_node::Array{Int}
    node_neighbors::Array{Int, 2}
    population::Array{Int}
    src_dst_to_edge::SparseMatrixCSC{Int, Int}
    edges::Array{Int, 2}
    attributes::Array{Dict{String, Any}}
end

"""
    IndexedGraph (raw, pop_col)

Generate an IndexedGraph from a NetworkX-format adjacency graph dictionary.
This is typically loaded from a JSON file.
"""
function IndexedGraph(raw::Dict, pop_col::AbstractString)::IndexedGraph
    # Generate the base SimpleGraph.
    n_nodes = length(raw["nodes"])
    population = zeros(Int, n_nodes)
    graph = SimpleGraph(n_nodes)
    for (index, node) in enumerate(raw["nodes"])
        population[index] = node[pop_col]
    end
    for (index, edges) in enumerate(raw["adjacency"])
        for edge in edges
            if edge["id"] + 1 > index
                add_edge!(graph, index, edge["id"] + 1)
            end
        end
    end

    # Build graph indices.
    n_edges = ne(graph)
    all_edges = zeros(Int, 2, n_edges)
    cut_edges = BitSet(Int[])

    src_dst_to_edge = spzeros(Int, n_nodes, n_nodes)
    for (index, edge) in enumerate(edges(graph))
        all_edges[1, index] = src(edge)
        all_edges[2, index] = dst(edge)
        src_dst_to_edge[src(edge), dst(edge)] = index
        src_dst_to_edge[dst(edge), src(edge)] = index
    end

    node_neighbors = zeros(Int, 20, n_nodes) # TODO: fix this constant
    neighbors_per_node = zeros(Int, n_nodes)
    for index in 1:n_nodes
        for (neighbor_idx, neighbor) in enumerate(neighbors(graph, index))
            node_neighbors[neighbor_idx, index] = neighbor
            neighbors_per_node[index] += 1
        end
    end

    # Attach node-level attribute data.
    attributes = Array{Dict{String, Any}}(undef, n_nodes)
    for (index, node) in enumerate(raw["nodes"])
        attributes[index] = node
    end

    return IndexedGraph(n_nodes, n_edges, neighbors_per_node, node_neighbors,
                        population, src_dst_to_edge, all_edges, attributes)
end

"""
    edge_index(graph, src, dst)

Return the index of the edge connecting node `src` to node `dst`.
"""
function edge_index(graph::IndexedGraph, src::Int, dst::Int)::Int
    return graph.src_dst_to_edge[src, dst]
end


text:
# v0.12 deprecations
@deprecate Dropout(p, dims) Dropout(p; dims=dims)
@deprecate InstanceNorm(λ, β, γ, μ, σ², ϵ, momentum) InstanceNorm(λ, β, γ, μ, σ², ϵ, momentum, true, true, nothing)
@deprecate BatchNorm(λ, β, γ, μ, σ², ϵ, momentum) BatchNorm(λ, β, γ, μ, σ², ϵ, momentum, true, true, nothing)
@deprecate GroupNorm(G, λ, β, γ, μ, σ², ϵ, momentum) GroupNorm(G, λ, β, γ, μ, σ², ϵ, momentum, nothing)
@deprecate outdims(f, inputsize) outputsize(f, inputsize)
@deprecate Conv(; weight,  bias, activation=identity, kws...) Conv(weight, bias, activation; kws...)
@deprecate ConvTranspose(; weight, bias, activation=identity, kws...) ConvTranspose(weight, bias, activation; kws...)
@deprecate DepthwiseConv(; weight, bias, activation=identity, kws...) DepthwiseConv(weight, bias, activation; kws...)

function Base.getproperty(a::Dense, s::Symbol)
  if s === :W
    Base.depwarn("field name dense.W is deprecated in favour of dense.weight", :Dense)
    return getfield(a, :weight)
  elseif s === :b
    Base.depwarn("field name dense.b is deprecated in favour of dense.bias", :Dense)
    return getfield(a, :bias)
  end
  return getfield(a, s)
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    Nonlinear

!!! warning
    The Nonlinear submodule is experimental. Until this message is removed,
    breaking changes may be introduced in any minor or patch release of
    MathOptInterface.
"""
module Nonlinear

import Base.Meta: isexpr
import ForwardDiff
import ..MathOptInterface
import OrderedCollections: OrderedDict
import SparseArrays

using SpecialFunctions

const MOI = MathOptInterface

# Override basic math functions to return NaN instead of throwing errors.
# This is what NLP solvers expect, and sometimes the results aren't needed
# anyway, because the code may compute derivatives wrt constants.
import NaNMath:
    sin,
    cos,
    tan,
    asin,
    acos,
    acosh,
    atanh,
    log,
    log2,
    log10,
    lgamma,
    log1p,
    pow,
    sqrt

include("univariate_expressions.jl")
include("operators.jl")
include("types.jl")
include("parse.jl")
include("model.jl")
include("evaluator.jl")

include("ReverseAD/ReverseAD.jl")

end  # module

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
const MISSING_KEY = NaN # use Inf for old version

const DEFAULT_DISCRETIZERS = Dict{Symbol,AbstractDiscretizer}()
    DEFAULT_DISCRETIZERS[:f_turnrate_250ms      ] = LinearDiscretizer([-0.025,-0.02,-0.015,-0.01,-0.005,0.005,0.01,0.015,0.02,0.025], Int)
    DEFAULT_DISCRETIZERS[:f_turnrate_500ms      ] = LinearDiscretizer([-0.025,-0.02,-0.015,-0.01,-0.005,0.005,0.01,0.015,0.02,0.025], Int)
    DEFAULT_DISCRETIZERS[:f_accel_250ms         ] = LinearDiscretizer([-5.00,-3.045338252335337,-2.0,-1.50593693678868,-0.24991770599882523,0.06206203400761478,0.2489269478410686,0.5,3.0], Int)
    DEFAULT_DISCRETIZERS[:f_accel               ] = LinearDiscretizer([-5.00,-3.045338252335337,-2.0,-1.50593693678868,-0.24991770599882523,0.06206203400761478,0.2489269478410686,0.5,3.0], Int)
    DEFAULT_DISCRETIZERS[:f_accel_500ms         ] = LinearDiscretizer([-1.0,-0.25,-0.08,0.0,0.08,0.25,1.0], Int)
    DEFAULT_DISCRETIZERS[:f_des_angle_250ms     ] = LinearDiscretizer([-0.25,-0.05,-0.03555530775035057,-0.02059935606012066,-0.005175260331415599,0.00642120613623413,0.02143002425291505,0.05,0.25], Int)
    DEFAULT_DISCRETIZERS[:f_des_angle           ] = LinearDiscretizer([-0.6,-0.05,-0.03555530775035057,-0.02059935606012066,-0.005175260331415599,0.00642120613623413,0.02143002425291505,0.05,0.75], Int)
    DEFAULT_DISCRETIZERS[:f_des_angle_500ms     ] = LinearDiscretizer([-0.025,-0.01,-0.005,0.005,0.01,0.025], Int)
    DEFAULT_DISCRETIZERS[:f_des_speed_250ms     ] = LinearDiscretizer([-2.0,-1.0,-0.5,0.0,0.5,1.0,2.0]+29.06, Int)
    DEFAULT_DISCRETIZERS[:f_des_speed_500ms     ] = LinearDiscretizer([-2.0,-1.0,-0.5,0.0,0.5,1.0,2.0]+29.06, Int)
    DEFAULT_DISCRETIZERS[:f_acc_control_250ms   ] = LinearDiscretizer([-0.2,0.0,0.05,0.2,0.4,0.6], Int)
    DEFAULT_DISCRETIZERS[:f_acc_control_500ms   ] = LinearDiscretizer([-0.2,0.0,0.05,0.2,0.4,0.6], Int)
    DEFAULT_DISCRETIZERS[:f_deltaY_250ms        ] = LinearDiscretizer([-0.4,-0.2,-0.1,0.1,0.2,0.4], Int)
    DEFAULT_DISCRETIZERS[:lanechange2s          ] = LinearDiscretizer([-0.5,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:yaw                   ] = LinearDiscretizer([-0.05,-0.02,-0.0065,0.0065,0.02,0.05], Int)
    DEFAULT_DISCRETIZERS[:posFyaw               ] = LinearDiscretizer([-0.05,-0.02,-0.0065,0.0065,0.02,0.05], Int)
    DEFAULT_DISCRETIZERS[:posFt                 ] = LinearDiscretizer([-2.5,-1.5,-1,1,1.5,2.5], Int)
    DEFAULT_DISCRETIZERS[:speed                 ] = LinearDiscretizer([0,10,20,22.5,25,30], Int)
    DEFAULT_DISCRETIZERS[:delta_speed_limit     ] = LinearDiscretizer([-1.0,-0.5,-0.0,0.5,1.0], Int)
    DEFAULT_DISCRETIZERS[:posFy                 ] = LinearDiscretizer([-2.5,-1.5,-0.5,0.5,1.5,2.5], Int)
    DEFAULT_DISCRETIZERS[:velFx                 ] = LinearDiscretizer([0,10,20,22.5,25,30], Int)
    DEFAULT_DISCRETIZERS[:velFy                 ] = LinearDiscretizer([-0.6,-0.4,-0.2,-0.1,0.1,0.2,0.4,0.6], Int)
    DEFAULT_DISCRETIZERS[:velBx                 ] = LinearDiscretizer([0,10,20,22.5,25,30], Int)
    DEFAULT_DISCRETIZERS[:velBy                 ] = LinearDiscretizer([-0.6,-0.4,-0.2,-0.1,0.1,0.2,0.4,0.6], Int)
    DEFAULT_DISCRETIZERS[:velFs                 ] = LinearDiscretizer([0,10,20,22.5,25,30], Int)
    DEFAULT_DISCRETIZERS[:velFt                 ] = LinearDiscretizer([-0.6,-0.4,-0.2,-0.1,0.1,0.2,0.4,0.6], Int)
    DEFAULT_DISCRETIZERS[:turnrate              ] = LinearDiscretizer([-0.04,-0.015,-0.005,0.005,0.015,0.04], Int)
    DEFAULT_DISCRETIZERS[:turnrate_global       ] = LinearDiscretizer([-0.04,-0.015,-0.005,0.005,0.015,0.04], Int)
    DEFAULT_DISCRETIZERS[:acc                   ] = LinearDiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int)
    DEFAULT_DISCRETIZERS[:accFx                 ] = LinearDiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int)
    DEFAULT_DISCRETIZERS[:accFy                 ] = LinearDiscretizer([-0.5,-0.15,-0.05,0.05,0.15,0.5], Int)
    DEFAULT_DISCRETIZERS[:accBx                 ] = LinearDiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int)
    DEFAULT_DISCRETIZERS[:accBy                 ] = LinearDiscretizer([-0.5,-0.15,-0.05,0.05,0.15,0.5], Int)
    DEFAULT_DISCRETIZERS[:accFs                 ] = LinearDiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int)
    DEFAULT_DISCRETIZERS[:accFt                 ] = LinearDiscretizer([-0.5,-0.15,-0.05,0.05,0.15,0.5], Int)
    DEFAULT_DISCRETIZERS[:cl                    ] = CategoricalDiscretizer([1.0,2.0,3.0,4.0], Int)
    DEFAULT_DISCRETIZERS[:d_cl                  ] = LinearDiscretizer([-2.5,-1.5,-1,1,1.5,2.5], Int)
    DEFAULT_DISCRETIZERS[:d_ml                  ] = LinearDiscretizer([-4,-3,-2,-1.5,-1,0], Int)
    DEFAULT_DISCRETIZERS[:d_mr                  ] = LinearDiscretizer([0,1,1.5,2,3,4], Int)
    DEFAULT_DISCRETIZERS[:d_merge               ] = LinearDiscretizer([0.0,50,75,100], Int)
    DEFAULT_DISCRETIZERS[:d_split               ] = LinearDiscretizer([0.0,50,75,100], Int)
    DEFAULT_DISCRETIZERS[:ttcr_ml               ] = datalineardiscretizer([0.0,6,9.6,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:ttcr_mr               ] = datalineardiscretizer([0.0,6,9.6,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:est_ttcr              ] = datalineardiscretizer([-10.0, -9.6, -6,0,6,9.6,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:scene_velFx           ] = LinearDiscretizer([0.0,10,20,22.5,25,30], Int)
    DEFAULT_DISCRETIZERS[:a_req_stayinlane      ] = datalineardiscretizer([0,0.05,0.15,2], Int, missing_key=NaN) # NOTE: used to not be this
    DEFAULT_DISCRETIZERS[:sumo                  ] = LinearDiscretizer([-5.00,-1.0,-0.5,0.0,0.5,3.0], Int)
    DEFAULT_DISCRETIZERS[:idm                   ] = LinearDiscretizer([-5.00,-1.0,-0.5,0.0,0.5,3.0], Int)
    DEFAULT_DISCRETIZERS[:n_lane_right          ] = CategoricalDiscretizer([0.0,1.0,2.0,3.0,4.0], Int)
    DEFAULT_DISCRETIZERS[:n_lane_left           ] = CategoricalDiscretizer([0.0,1.0,2.0,3.0,4.0], Int)
    DEFAULT_DISCRETIZERS[:has_lane_left         ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:has_lane_right        ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:curvature             ] = LinearDiscretizer([-0.002,-0.0005,0.0005,0.002], Int)
    DEFAULT_DISCRETIZERS[:has_front             ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:d_y_front             ] = datalineardiscretizer([-2,-0.5,0.5,2], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:v_x_front             ] = datalineardiscretizer([-1.5,-0.25,-0.1,0.1,0.25,1.5], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:d_x_front             ] = datalineardiscretizer([0.0,5,10,30,100], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:v_y_front             ] = datalineardiscretizer([-3,-1.5,-0.5,0.5,1.5,3], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:yaw_front             ] = datalineardiscretizer([-0.1,-0.02,0.02,0.1], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:a_req_front           ] = datalineardiscretizer([-2,-1,-0.25,0], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:turnrate_front        ] = datalineardiscretizer([-0.04,-0.015,-0.005,0.005,0.015,0.04], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:ttc_x_front           ] = datalineardiscretizer([0.0,2,4,9,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:timegap_x_front       ] = datalineardiscretizer([0.0,2,4,9,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:has_rear              ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:d_y_rear              ] = datalineardiscretizer([-2,-0.5,0.5,2], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:v_x_rear              ] = datalineardiscretizer([-1.5,-0.25,-0.1,0.1,0.25,1.5], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:d_x_rear              ] = datalineardiscretizer([0.0,5,10,30,100], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:v_y_rear              ] = datalineardiscretizer([-3,-1.5,-0.5,0.5,1.5,3], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:yaw_rear              ] = datalineardiscretizer([-0.1,-0.02,0.02,0.1], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:a_req_rear            ] = datalineardiscretizer([0.0,1,2,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:turnrate_rear         ] = datalineardiscretizer([-0.04,-0.015,-0.005,0.005,0.015,0.04], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:ttc_x_rear            ] = datalineardiscretizer([0.0,2,4,8,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:timegap_x_rear        ] = datalineardiscretizer([0.0,2,4,9,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:d_y_left              ] = datalineardiscretizer([0.0,2,5,8], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:v_x_left              ] = datalineardiscretizer([-1.5,-0.25,-0.1,0.1,0.25,1.5], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:d_x_left              ] = datalineardiscretizer([-50.0,-20,-10,10,20,50], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:v_y_left              ] = datalineardiscretizer([-3,-1.5,-0.5,0.5,1.5,3], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:yaw_left              ] = datalineardiscretizer([-0.1,-0.02,0.02,0.1], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:a_req_left            ] = datalineardiscretizer([-2,-1,-0.05,0.05,1,2], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:turnrate_left         ] = datalineardiscretizer([-0.04,-0.015,-0.005,0.005,0.015,0.04], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:ttc_x_left            ] = datalineardiscretizer([0.0,2,4,8,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:timegap_x_left        ] = datalineardiscretizer([0.0,2,4,9,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:d_y_right             ] = datalineardiscretizer([-8.0,-5,-2,0], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:v_x_right             ] = datalineardiscretizer([-1.5,-0.25,-0.1,0.1,0.25,1.5], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:d_x_right             ] = datalineardiscretizer([-50.0,-20,-10,10,20,50], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:v_y_right             ] = datalineardiscretizer([-3,-1.5,-0.5,0.5,1.5,3], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:yaw_right             ] = datalineardiscretizer([-0.1,-0.02,0.02,0.1], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:a_req_right           ] = datalineardiscretizer([-2,-1,-0.05,0.05,1,2], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:turnrate_right        ] = datalineardiscretizer([-0.04,-0.015,-0.005,0.005,0.015,0.04], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:ttc_x_right           ] = datalineardiscretizer([0.0,2,4,8,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:timegap_x_right       ] = datalineardiscretizer([0.0,2,4,9,10], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:time_consecutive_brake] = datalineardiscretizer([0.0,0.251,0.501,1.001,2.0], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:time_consecutive_accel] = datalineardiscretizer([0.0,0.251,0.501,1.001,2.0], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:time_consecutive_throttle] = datalineardiscretizer([-2.0,-1.001,-0.501,-0.251,0.0,0.251,0.501,1.001,2.0], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:timesincelanecrossing ] = datalineardiscretizer([0.0,0.501,10.0], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:osg_isoccupied_f      ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:osg_isoccupied_fr     ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:osg_isoccupied_r      ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:osg_isoccupied_br     ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:osg_isoccupied_b      ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:osg_isoccupied_bl     ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:osg_isoccupied_l      ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:osg_isoccupied_fl     ] = CategoricalDiscretizer([0.0,1.0], Int)
    DEFAULT_DISCRETIZERS[:osg_time_f            ] = LinearDiscretizer([0.0,4.0,8.0,12.0], Int)
    DEFAULT_DISCRETIZERS[:osg_time_fr           ] = LinearDiscretizer([0.0,4.0,8.0,12.0], Int)
    DEFAULT_DISCRETIZERS[:osg_time_r            ] = LinearDiscretizer([0.0,4.0,8.0,12.0], Int)
    DEFAULT_DISCRETIZERS[:osg_time_br           ] = LinearDiscretizer([0.0,4.0,8.0,12.0], Int)
    DEFAULT_DISCRETIZERS[:osg_time_b            ] = LinearDiscretizer([0.0,4.0,8.0,12.0], Int)
    DEFAULT_DISCRETIZERS[:osg_time_bl           ] = LinearDiscretizer([0.0,4.0,8.0,12.0], Int)
    DEFAULT_DISCRETIZERS[:osg_time_l            ] = LinearDiscretizer([0.0,4.0,8.0,12.0], Int)
    DEFAULT_DISCRETIZERS[:osg_time_fl           ] = LinearDiscretizer([0.0,4.0,8.0,12.0], Int)
    DEFAULT_DISCRETIZERS[:pastacc250ms          ] = LinearDiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int)
    DEFAULT_DISCRETIZERS[:pastacc500ms          ] = LinearDiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int)
    DEFAULT_DISCRETIZERS[:pastacc750ms          ] = LinearDiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int)
    DEFAULT_DISCRETIZERS[:pastacc1s             ] = LinearDiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int)

    DEFAULT_DISCRETIZERS[:past_1_f_accel        ] = datalineardiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:past_2_f_accel        ] = datalineardiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:past_3_f_accel        ] = datalineardiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:past_4_f_accel        ] = datalineardiscretizer([-1,-0.25,-0.08,0.08,0.25,1], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:past_1_f_des_angle    ] = datalineardiscretizer([-0.6,-0.05,-0.035,-0.021,-0.005,0.006,0.02,0.05,0.75], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:past_2_f_des_angle    ] = datalineardiscretizer([-0.6,-0.05,-0.035,-0.021,-0.005,0.006,0.02,0.05,0.75], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:past_3_f_des_angle    ] = datalineardiscretizer([-0.6,-0.05,-0.035,-0.021,-0.005,0.006,0.02,0.05,0.75], Int, missing_key=NaN)
    DEFAULT_DISCRETIZERS[:past_4_f_des_angle    ] = datalineardiscretizer([-0.6,-0.05,-0.035,-0.021,-0.005,0.006,0.02,0.05,0.75], Int, missing_key=NaN)

    DEFAULT_DISCRETIZERS[:pastturnrate250ms     ] = LinearDiscretizer([-0.02,-0.01,-0.005,0.005,0.01,0.02], Int)
    DEFAULT_DISCRETIZERS[:pastturnrate500ms     ] = LinearDiscretizer([-0.02,-0.01,-0.005,0.005,0.01,0.02], Int)
    DEFAULT_DISCRETIZERS[:pastturnrate750ms     ] = LinearDiscretizer([-0.02,-0.01,-0.005,0.005,0.01,0.02], Int)
    DEFAULT_DISCRETIZERS[:pastturnrate1s        ] = LinearDiscretizer([-0.02,-0.01,-0.005,0.005,0.01,0.02], Int)
    DEFAULT_DISCRETIZERS[:pastvelFy250ms        ] = LinearDiscretizer([-1.0,-0.5,-0.1,0.1,0.5,1.0], Int)
    DEFAULT_DISCRETIZERS[:pastvelFy500ms        ] = LinearDiscretizer([-1.0,-0.5,-0.1,0.1,0.5,1.0], Int)
    DEFAULT_DISCRETIZERS[:pastvelFy750ms        ] = LinearDiscretizer([-1.0,-0.5,-0.1,0.1,0.5,1.0], Int)
    DEFAULT_DISCRETIZERS[:pastvelFy1s           ] = LinearDiscretizer([-1.0,-0.5,-0.1,0.1,0.5,1.0], Int)
    DEFAULT_DISCRETIZERS[:pastd_cl250ms         ] = LinearDiscretizer([-0.3,-0.2,-0.1,0.1,0.2,0.3], Int)
    DEFAULT_DISCRETIZERS[:pastd_cl500ms         ] = LinearDiscretizer([-0.5,-0.3,-0.15,0.15,0.3,0.5], Int)
    DEFAULT_DISCRETIZERS[:pastd_cl750ms         ] = LinearDiscretizer([-0.75,-0.5,-0.25,0.25,0.5,0.75], Int)
    DEFAULT_DISCRETIZERS[:pastd_cl1s            ] = LinearDiscretizer([-0.75,-0.5,-0.25,0.25,0.5,0.75], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx250ms         ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx500ms         ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx750ms         ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx1s            ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx1500ms        ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx2s            ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx2500ms        ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx3s            ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFx4s            ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy250ms         ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy500ms         ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy750ms         ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy1s            ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy1500ms        ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy2s            ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy2500ms        ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy3s            ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxaccFy4s            ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate250ms      ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate500ms      ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate750ms      ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate1s         ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate1500ms     ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate2s         ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate2500ms     ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate3s         ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:maxturnrate4s         ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx250ms        ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx500ms        ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx750ms        ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx1s           ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx1500ms       ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx2s           ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx2500ms       ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx3s           ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFx4s           ] = LinearDiscretizer([-1.5,-0.5,-0.15,0.15,0.5,1.5], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy100ms        ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy150ms        ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy200ms        ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy250ms        ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy500ms        ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy750ms        ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy1s           ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy1500ms       ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy2s           ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy2500ms       ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy3s           ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanaccFy4s           ] = LinearDiscretizer([-2,-0.5,-0.15,0.15,0.5,2], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate250ms     ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate500ms     ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate750ms     ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate1s        ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate1500ms    ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate2s        ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate2500ms    ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate3s        ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:meanturnrate4s        ] = LinearDiscretizer([-0.1,-0.03,-0.004,0.004,0.03,0.1], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx250ms         ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx500ms         ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx750ms         ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx1s            ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx1500ms        ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx2s            ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx2500ms        ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx3s            ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFx4s            ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy250ms         ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy500ms         ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy750ms         ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy1s            ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy1500ms        ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy2s            ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy2500ms        ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy3s            ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdaccFy4s            ] = LinearDiscretizer([0,0.1,0.2,0.5], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate250ms      ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate500ms      ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate750ms      ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate1s         ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate1500ms     ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate2s         ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate2500ms     ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate3s         ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)
    DEFAULT_DISCRETIZERS[:stdturnrate4s         ] = LinearDiscretizer([0,0.01,0.02,0.04], Int)


text:
__precompile__()
module MathOptInterfaceCPLEX

export CPLEXOptimizer

using CPLEX
const CPX = CPLEX
using MathOptInterface
const MOI = MathOptInterface
using LinQuadOptInterface
const LQOI = LinQuadOptInterface

const SUPPORTED_OBJECTIVES = [
    LQOI.Linear,
    LQOI.Quad
]

const SUPPORTED_CONSTRAINTS = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    (LQOI.Linear, LQOI.IV),
    (LQOI.Quad, LQOI.EQ),
    (LQOI.Quad, LQOI.LE),
    (LQOI.Quad, LQOI.GE),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.SinVar, MOI.ZeroOne),
    (LQOI.SinVar, MOI.Integer),
    (LQOI.VecVar, LQOI.SOS1),
    (LQOI.VecVar, LQOI.SOS2),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros)
]

mutable struct CPLEXOptimizer <: LQOI.LinQuadOptimizer
    LQOI.@LinQuadOptimizerBase
    env
    params::Dict{String,Any}
    CPLEXOptimizer(::Void) = new()
end

LQOI.LinQuadModel(::Type{CPLEXOptimizer},env) = CPX.Model(env)

function CPLEXOptimizer(;kwargs...)

    env = CPX.Env()
    m = CPLEXOptimizer(nothing)
    m.env = env
    m.params = Dict{String,Any}()
    MOI.empty!(m)
    for (name,value) in kwargs
        m.params[string(name)] = value
        CPX.cpx_setparam!(m.inner.env, string(name), value)
    end
    return m
end

function MOI.empty!(m::CPLEXOptimizer)
    MOI.empty!(m,m.env)
    for (name,value) in m.params
        CPX.cpx_setparam!(m.inner.env, string(name), value)
    end
end

LQOI.lqs_supported_constraints(s::CPLEXOptimizer) = SUPPORTED_CONSTRAINTS
LQOI.lqs_supported_objectives(s::CPLEXOptimizer) = SUPPORTED_OBJECTIVES

#=
    inner wrapper
=#

#=
    Main
=#

# LinQuadSolver # Abstract type
# done above

# LQOI.lqs_setparam!(env, name, val)
# TODO fix this one
LQOI.lqs_setparam!(m::CPLEXOptimizer, name, val) = CPX.cpx_setparam!(m.inner, string(name), val)

# LQOI.lqs_setlogfile!(env, path)
# TODO fix this one
LQOI.lqs_setlogfile!(m::CPLEXOptimizer, path) = CPX.setlogfile(m.inner, path::String)

# LQOI.lqs_getprobtype(m)
# TODO - consider removing, apparently useless

#=
    Constraints
=#

# LQOI.lqs_chgbds!(m, colvec, valvec, sensevec)
LQOI.lqs_chgbds!(instance::CPLEXOptimizer, colvec, valvec, sensevec) = CPX.cpx_chgbds!(instance.inner, colvec, valvec, sensevec)

# LQOI.lqs_getlb(m, col)
LQOI.lqs_getlb(instance::CPLEXOptimizer, col) = CPX.cpx_getlb(instance.inner, col)
# LQOI.lqs_getub(m, col)
LQOI.lqs_getub(instance::CPLEXOptimizer, col) = CPX.cpx_getub(instance.inner, col)

# LQOI.lqs_getnumrows(m)
LQOI.lqs_getnumrows(instance::CPLEXOptimizer) = CPX.cpx_getnumrows(instance.inner)

# LQOI.lqs_addrows!(m, rowvec, colvec, coefvec, sensevec, rhsvec)
LQOI.lqs_addrows!(instance::CPLEXOptimizer, rowvec, colvec, coefvec, sensevec, rhsvec) = CPX.cpx_addrows!(instance.inner, rowvec, colvec, coefvec, sensevec, rhsvec)

# LQOI.lqs_getrhs(m, rowvec)
LQOI.lqs_getrhs(instance::CPLEXOptimizer, row) = CPX.cpx_getrhs(instance.inner, row)

# colvec, coef = LQOI.lqs_getrows(m, rowvec)
# TODO improve
function LQOI.lqs_getrows(instance::CPLEXOptimizer, idx)
    return CPX.cpx_getrows(instance.inner, idx)
end

# LQOI.lqs_getcoef(m, row, col) #??
# TODO improve
LQOI.lqs_getcoef(instance::CPLEXOptimizer, row, col) = CPX.cpx_getcoef(instance.inner, row, col)

# LQOI.lqs_chgcoef!(m, row, col, coef)
# TODO SPLIT THIS ONE
LQOI.lqs_chgcoef!(instance::CPLEXOptimizer, row, col, coef)  = CPX.cpx_chgcoef!(instance.inner, row, col, coef)

# LQOI.lqs_delrows!(m, row, row)
LQOI.lqs_delrows!(instance::CPLEXOptimizer, rowbeg, rowend) = CPX.cpx_delrows!(instance.inner, rowbeg, rowend)

# LQOI.lqs_chgctype!(m, colvec, typevec)
# TODO fix types
LQOI.lqs_chgctype!(instance::CPLEXOptimizer, colvec, typevec) = CPX.cpx_chgctype!(instance.inner, colvec, typevec)

# LQOI.lqs_chgsense!(m, rowvec, sensevec)
# TODO fix types
LQOI.lqs_chgsense!(instance::CPLEXOptimizer, rowvec, sensevec) = CPX.cpx_chgsense!(instance.inner, rowvec, sensevec)

const VAR_TYPE_MAP = Dict{Symbol,Cchar}(
    :CONTINUOUS => Cchar('C'),
    :INTEGER => Cchar('I'),
    :BINARY => Cchar('B')
)
LQOI.lqs_vartype_map(m::CPLEXOptimizer) = VAR_TYPE_MAP

# LQOI.lqs_addsos(m, colvec, valvec, typ)
LQOI.lqs_addsos!(instance::CPLEXOptimizer, colvec, valvec, typ) = CPX.add_sos!(instance.inner, typ, colvec, valvec)
# LQOI.lqs_delsos(m, idx, idx)
LQOI.lqs_delsos!(instance::CPLEXOptimizer, idx1, idx2) = CPX.cpx_delsos!(instance.inner, idx1, idx2)

const SOS_TYPE_MAP = Dict{Symbol,Symbol}(
    :SOS1 => :SOS1,#Cchar('1'),
    :SOS2 => :SOS2#Cchar('2')
)
LQOI.lqs_sertype_map(m::CPLEXOptimizer) = SOS_TYPE_MAP

# LQOI.lqs_getsos(m, idx)
# TODO improve getting processes
function LQOI.lqs_getsos(instance::CPLEXOptimizer, idx)
    indices, weights, types = CPX.cpx_getsos(instance.inner, idx)

    # types2 = Array{Symbol}(length(types))
    # for i in eachindex(types)
    #     if types[i] == Cchar('1')
    #         types2[i] = :SOS1
    #     elseif types[i] == Cchar('2')
    #         types2[i] = :SOS2
    #     end
    # end

    return indices, weights, types == Cchar('1') ? :SOS1 : :SOS2
end
# LQOI.lqs_getnumqconstrs(m)
LQOI.lqs_getnumqconstrs(instance::CPLEXOptimizer) = CPX.cpx_getnumqconstrs(instance.inner)

# LQOI.lqs_addqconstr(m, cols,coefs,rhs,sense, I,J,V)
LQOI.lqs_addqconstr!(instance::CPLEXOptimizer, cols,coefs,rhs,sense, I,J,V) = CPX.cpx_addqconstr!(instance.inner, cols,coefs,rhs,sense, I,J,V)

# LQOI.lqs_chgrngval
LQOI.lqs_chgrngval!(instance::CPLEXOptimizer, rows, vals) = CPX.cpx_chgrngval!(instance.inner, rows, vals)

const CTR_TYPE_MAP = Dict{Symbol,Cchar}(
    :RANGE => Cchar('R'),
    :LOWER => Cchar('L'),
    :UPPER => Cchar('U'),
    :EQUALITY => Cchar('E')
)
LQOI.lqs_ctrtype_map(m::CPLEXOptimizer) = CTR_TYPE_MAP

#=
    Objective
=#

# LQOI.lqs_copyquad(m, intvec,intvec, floatvec) #?
LQOI.lqs_copyquad!(instance::CPLEXOptimizer, I, J, V) = CPX.cpx_copyquad!(instance.inner, I, J, V)

# LQOI.lqs_chgobj(m, colvec,coefvec)
function LQOI.lqs_chgobj!(instance::CPLEXOptimizer, colvec, coefvec)
    ncols = CPX.cpx_getnumcols(instance.inner)
    new_colvec = collect(1:ncols)
    new_coefvec = zeros(ncols)
    for (ind,val) in enumerate(colvec)
        new_coefvec[val] = coefvec[ind]
    end
    # this only sums to obj
    CPX.cpx_chgobj!(instance.inner, new_colvec, new_coefvec)
end
# LQOI.lqs_chgobjsen(m, symbol)
# TODO improve min max names
LQOI.lqs_chgobjsen!(instance::CPLEXOptimizer, symbol) = CPX.cpx_chgobjsen!(instance.inner, symbol)

# LQOI.lqs_getobj(m)
LQOI.lqs_getobj(instance::CPLEXOptimizer) = CPX.cpx_getobj(instance.inner)

# lqs_getobjsen(m)
LQOI.lqs_getobjsen(instance::CPLEXOptimizer) = CPX.cpx_getobjsen(instance.inner)

#=
    Variables
=#

# LQOI.lqs_getnumcols(m)
LQOI.lqs_getnumcols(instance::CPLEXOptimizer) = CPX.cpx_getnumcols(instance.inner)

# LQOI.lqs_newcols!(m, int)
LQOI.lqs_newcols!(instance::CPLEXOptimizer, int) = CPX.cpx_newcols!(instance.inner, int)

# LQOI.lqs_delcols!(m, col, col)
LQOI.lqs_delcols!(instance::CPLEXOptimizer, col, col2) = CPX.cpx_delcols!(instance.inner, col, col2)

# LQOI.lqs_addmipstarts(m, colvec, valvec)
LQOI.lqs_addmipstarts!(instance::CPLEXOptimizer, colvec, valvec)  = CPX.cpx_addmipstarts!(instance.inner, colvec, valvec)

#=
    Solve
=#

# LQOI.lqs_mipopt!(m)
LQOI.lqs_mipopt!(instance::CPLEXOptimizer) = CPX.cpx_mipopt!(instance.inner)

# LQOI.lqs_qpopt!(m)
LQOI.lqs_qpopt!(instance::CPLEXOptimizer) = CPX.cpx_qpopt!(instance.inner)

# LQOI.lqs_lpopt!(m)
LQOI.lqs_lpopt!(instance::CPLEXOptimizer) = CPX.cpx_lpopt!(instance.inner)

const TERMINATION_STATUS_MAP = Dict(
    CPX.CPX_STAT_OPTIMAL                => MOI.Success,
    CPX.CPX_STAT_UNBOUNDED              => MOI.UnboundedNoResult,
    CPX.CPX_STAT_INFEASIBLE             => MOI.InfeasibleNoResult,
    CPX.CPX_STAT_INForUNBD              => MOI.InfeasibleOrUnbounded,
    CPX.CPX_STAT_OPTIMAL_INFEAS         => MOI.Success,
    CPX.CPX_STAT_NUM_BEST               => MOI.NumericalError,
    CPX.CPX_STAT_ABORT_IT_LIM           => MOI.IterationLimit,
    CPX.CPX_STAT_ABORT_TIME_LIM         => MOI.TimeLimit,
    CPX.CPX_STAT_ABORT_OBJ_LIM          => MOI.ObjectiveLimit,
    CPX.CPX_STAT_ABORT_USER             => MOI.Interrupted,
    CPX.CPX_STAT_OPTIMAL_FACE_UNBOUNDED => MOI.UnboundedNoResult,
    CPX.CPX_STAT_ABORT_PRIM_OBJ_LIM     => MOI.ObjectiveLimit,
    CPX.CPX_STAT_ABORT_DUAL_OBJ_LIM     => MOI.ObjectiveLimit,
    CPX.CPXMIP_OPTIMAL                  => MOI.Success,
    CPX.CPXMIP_OPTIMAL_TOL              => MOI.Success,
    CPX.CPXMIP_INFEASIBLE               => MOI.InfeasibleNoResult,
    CPX.CPXMIP_SOL_LIM                  => MOI.SolutionLimit,
    CPX.CPXMIP_NODE_LIM_FEAS            => MOI.NodeLimit,
    CPX.CPXMIP_NODE_LIM_INFEAS          => MOI.NodeLimit,
    CPX.CPXMIP_TIME_LIM_FEAS            => MOI.TimeLimit,
    CPX.CPXMIP_TIME_LIM_INFEAS          => MOI.TimeLimit,
    CPX.CPXMIP_FAIL_FEAS                => MOI.OtherError,
    CPX.CPXMIP_FAIL_INFEAS              => MOI.OtherError,
    CPX.CPXMIP_MEM_LIM_FEAS             => MOI.MemoryLimit,
    CPX.CPXMIP_MEM_LIM_INFEAS           => MOI.MemoryLimit,
    CPX.CPXMIP_ABORT_FEAS               => MOI.Interrupted,
    CPX.CPXMIP_ABORT_INFEAS             => MOI.Interrupted,
    CPX.CPXMIP_OPTIMAL_INFEAS           => MOI.Success,
    CPX.CPXMIP_FAIL_FEAS_NO_TREE        => MOI.MemoryLimit,
    CPX.CPXMIP_FAIL_INFEAS_NO_TREE      => MOI.MemoryLimit,
    CPX.CPXMIP_UNBOUNDED                => MOI.UnboundedNoResult,
    CPX.CPXMIP_INForUNBD                => MOI.InfeasibleOrUnbounded
)

# LQOI.lqs_terminationstatus(m)
function LQOI.lqs_terminationstatus(model::CPLEXOptimizer)
    m = model.inner

    code = CPX.cpx_getstat(m)
    mthd, soltype, prifeas, dualfeas = CPX.cpx_solninfo(m)

    if haskey(TERMINATION_STATUS_MAP, code)
        out = TERMINATION_STATUS_MAP[code]

        if code == CPX.CPX_STAT_UNBOUNDED && prifeas > 0
            out = MOI.Success
        elseif code == CPX.CPX_STAT_INFEASIBLE && dualfeas > 0
            out = MOI.Success
        end
        return out
    else
        error("Status $(code) has not been mapped to a MOI termination status.")
    end
end

function LQOI.lqs_primalstatus(model::CPLEXOptimizer)
    m = model.inner

    code = CPX.cpx_getstat(m)
    mthd, soltype, prifeas, dualfeas = CPX.cpx_solninfo(m)

    out = MOI.UnknownResultStatus

    if soltype in [CPX.CPX_NONBASIC_SOLN, CPX.CPX_BASIC_SOLN, CPX.CPX_PRIMAL_SOLN]
        if prifeas > 0
            out = MOI.FeasiblePoint
        else
            out = MOI.InfeasiblePoint
        end
    end
    if code == CPX.CPX_STAT_UNBOUNDED #&& prifeas > 0
        out = MOI.InfeasibilityCertificate
    end
    return out
end
function LQOI.lqs_dualstatus(model::CPLEXOptimizer)
    m = model.inner

    code = CPX.cpx_getstat(m)
    mthd, soltype, prifeas, dualfeas = CPX.cpx_solninfo(m)
    if !LQOI.hasinteger(model)
        if soltype in [CPX.CPX_NONBASIC_SOLN, CPX.CPX_BASIC_SOLN]
            if dualfeas > 0
                out = MOI.FeasiblePoint
            else
                out = MOI.InfeasiblePoint
            end
        else
            out = MOI.UnknownResultStatus
        end
        if code == CPX.CPX_STAT_INFEASIBLE && dualfeas > 0
            out = MOI.InfeasibilityCertificate
        end
        return out
    end
    return MOI.UnknownResultStatus
end

# LQOI.lqs_getx!(m, place)
LQOI.lqs_getx!(instance::CPLEXOptimizer, place) = CPX.cpx_getx!(instance.inner, place)

# LQOI.lqs_getax!(m, place)
LQOI.lqs_getax!(instance::CPLEXOptimizer, place) = CPX.cpx_getax!(instance.inner, place)

# LQOI.lqs_getdj!(m, place)
LQOI.lqs_getdj!(instance::CPLEXOptimizer, place) = CPX.cpx_getdj!(instance.inner, place)

# LQOI.lqs_getpi!(m, place)
LQOI.lqs_getpi!(instance::CPLEXOptimizer, place) = CPX.cpx_getpi!(instance.inner, place)

# LQOI.lqs_getobjval(m)
LQOI.lqs_getobjval(instance::CPLEXOptimizer) = CPX.cpx_getobjval(instance.inner)

# LQOI.lqs_getbestobjval(m)
LQOI.lqs_getbestobjval(instance::CPLEXOptimizer) = CPX.cpx_getbestobjval(instance.inner)

# LQOI.lqs_getmiprelgap(m)
LQOI.lqs_getmiprelgap(instance::CPLEXOptimizer) = CPX.cpx_getmiprelgap(instance.inner)

# LQOI.lqs_getitcnt(m)
LQOI.lqs_getitcnt(instance::CPLEXOptimizer)  = CPX.cpx_getitcnt(instance.inner)

# LQOI.lqs_getbaritcnt(m)
LQOI.lqs_getbaritcnt(instance::CPLEXOptimizer) = CPX.cpx_getbaritcnt(instance.inner)

# LQOI.lqs_getnodecnt(m)
LQOI.lqs_getnodecnt(instance::CPLEXOptimizer) = CPX.cpx_getnodecnt(instance.inner)

# LQOI.lqs_dualfarkas(m, place)
LQOI.lqs_dualfarkas!(instance::CPLEXOptimizer, place) = CPX.cpx_dualfarkas!(instance.inner, place)

# LQOI.lqs_getray(m, place)
LQOI.lqs_getray!(instance::CPLEXOptimizer, place) = CPX.cpx_getray!(instance.inner, place)

MOI.free!(instance::CPLEXOptimizer) = CPX.free_model(instance.inner)

"""
    writeproblem(m: :MOI.AbstractOptimizer, filename::String)
Writes the current problem data to the given file.
Supported file types are solver-dependent.
"""
writeproblem(instance::CPLEXOptimizer, filename::String, flags::String="") = CPX.write_model(instance.inner, filename)

LQOI.lqs_make_problem_type_continuous(instance::CPLEXOptimizer) = CPX._make_problem_type_continuous(instance.inner)
end # module


text:
import GPUArrays: allowscalar, @allowscalar

function _getindex(xs::CuArray{T}, i::Integer) where T
  buf = Mem.view(buffer(xs), (i-1)*sizeof(T))
  return Mem.download(T, buf)[1]
end

function _setindex!(xs::CuArray{T}, v::T, i::Integer) where T
  buf = Mem.view(buffer(xs), (i-1)*sizeof(T))
  Mem.upload!(buf, T[v])
end


text:
mutable struct FriendSets
    var::Vector{Symbol}
    sets::Vector{String}
    counts::Vector{Int64}
    items::Vector{Vector{Int64}}
    FriendSets(var, sets, counts, items) = new(var, sets, counts, items)
    FriendSets() = new(Symbol[], String[], Int64[], Vector{Vector{Int64}}(undef, 0))
end


text:
includet(pwd()*"/src/functions/hull2array.jl")
includet(pwd()*"/src/functions/calcSim.jl")
includet(pwd()*"/src/functions/memoryPart.jl")

function give_me_my_name(f::Function)
           return String(Symbol(f))
       end

function getdists()
    # Distance functions from the Distances package
    dists = unique(sort(String.(names(Distances))))
    # Fiter and remove those starting with uppercase and _
    dists = dists[collect(.!(isuppercase(only(dists[n][1])) || dists[n][1] == "_" for n in 1:length(dists)))]

    # functions to remove
    removable = ["bregman" "colwise" "colwise!" "evaluate" "haversine" "mahalanobis" "minkowski" "pairwise" "pairwise!" "peuclidean" "renyi_divergence" "result_type" "rogerstanimoto" "spherical_angle" "sqmahalanobis" "wminkowski" "wcityblock" "weuclidean" "whamming" "wsqeuclidean"]
    for r in removable
        dists = dists[collect(.!(dists[n] == r for n in 1:length(dists)))]
    end

    # weigthed metrics
    weighted = ["wcityblock" "wcityblockvar" "wcityblockinv" "weuclidean" "weuclideanvar" "weuclideaninv" "whamming" "whammingvar" "whamminginv" "wsqeuclidean" "wsqeuclideanvar" "wsqeuclideaninv"]
    for w in weighted
        push!(dists, w)
    end

        return dists
    end

files = ["UCData" "AffProp" "GAPM"  "GAPM_Vec" "mUCEval" "plotAffProp" "UC" "UCEvaluation"];
    for f in files
        includet("functions/"*f*".jl")
    end

function evaluateScenario(ω)
        # println(ω)
        for t in T
            set_normalized_rhs(PowerBalance_Second[t], -(SystemLoad[t] - sum(RES[t,ω])) )
            for g in G
                set_normalized_rhs(PGen_2nd[g,t], Pexp[g,t])
                set_normalized_rhs(PMax2nd[g,t], V[g,t]*GenData.Pmax[g])
                set_normalized_rhs(PMin2nd[g,t], -V[g,t]*GenData.Pmin[g])
                set_normalized_rhs(RplusLim[g,t], Rup[g,t])
                set_normalized_rhs(RminusLim[g,t], Rdw[g,t])
                end
            end

            # Evaluating 2nd stage
            optimize!(mUCEval)
            termination_status(mUCEval)
            OpCost[ω] = JuMP.objective_value.(mUCEval)

            if ALL_DUALS == false # Use only elec balance dual
                λ[ω,:] =  -collect(shadow_price.(PowerBalance_Second[t]) for t in T)
                else
                    λ[ω,:] =
                        [reshape(-collect(shadow_price.(PowerBalance_Second[t]) for t in T),1,:)                                                    reshape(-collect(shadow_price.(PMax2nd[g,t]) for g in G, t in T),1,:)                                                    reshape(-collect(shadow_price.(PMin2nd[g,t]) for g in G, t in T),1,:)                                                    reshape(-collect(shadow_price.(RplusLim[g,t]) for g in G, t in T),1,:)                                                    reshape(-collect(shadow_price.(RminusLim[g,t]) for g in G, t in T),1,:)                                                    reshape(-collect(shadow_price.(Ramps_2ndUp[g,t]) for g in G, t in T),1,:)                                                    reshape(-collect(shadow_price.(Ramps_2ndDw[g,t]) for g in G, t in T),1,:)                                                    reshape(-collect(shadow_price.(PGen_2nd[g,t]) for g in G, t in T),1,:)
                                ]

                end
end

function latex_sci_not( x , ndec; font="sf" )
  xchar = strip(Formatting.sprintf1("%17.$(ndec)e",x))
  data = split(xchar,"e")
  inonzero = findfirst( i -> i != '0', data[2][2:length(data[2])])
  if font == "sf"
    f = "\\textrm{\\sffamily "
    fe = "\\textrm{\\sffamily\\scriptsize "
  else
    f = "{"
    fe = "{"
  end
  if inonzero == nothing
    string = latexstring("$f$(data[1])}")
  else
    if data[2][1] == '-'
      string = latexstring("$f$(data[1])}\\times $f 10}^{$fe$(data[2][1])$(data[2][inonzero+1:length(data[2])])}}")
    else
      string = latexstring("$f$(data[1])}\\times $f 10}^{$fe$(data[2][inonzero+1:length(data[2])])}}")
    end
  end
  return string
    end

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
module ZEON
using ..AGFFileReader: Glass, GlassID, AGF
export E48R, _480R, _330R

const E48R = Glass(GlassID(AGF, 1776), 1, 2.30319693, -0.00117825616, 0.0140597151, 2.97687093e-5, 3.28913151e-5, -1.61443145e-6, NaN, NaN, NaN, NaN, 0.4, 1.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.0, -1.0, -1.0, 60.0, -1.0, 0.0, -1.0, [(0.4, 0.989, 3.0), (0.42, 0.993, 3.0), (0.44, 0.996, 3.0), (0.46, 0.998, 3.0), (0.48, 0.999, 3.0), (0.5, 1.0, 3.0), (0.55, 1.0, 3.0), (0.6, 1.0, 3.0), (0.65, 1.0, 3.0), (0.7, 1.0, 3.0), (0.75, 0.999, 3.0), (0.8, 1.0, 3.0), (0.9, 0.987, 3.0), (1.0, 0.997, 3.0), (1.1, 0.996, 3.0), (1.2, 0.595, 3.0), (1.3, 0.97, 3.0), (1.4, 0.78, 3.0), (1.5, 0.909, 3.0), (1.6, 0.948, 3.0)], 1.53117, -1.0, -1.0, 0.0, 56.0, 0.0, 1.01, 0.0)
const _480R = Glass(GlassID(AGF, 1777), 1, 2.28919759, -0.00190608244, 0.0107336897, 0.00116337387, -0.000142408784, 7.73782254e-6, NaN, NaN, NaN, NaN, 0.4, 1.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.0, -11.0, -1.0, 60.0, -1.0, 0.0, -1.0, [(0.4, 0.99, 3.0), (0.42, 0.991, 3.0), (0.44, 0.993, 3.0), (0.46, 0.994, 3.0), (0.48, 0.994, 3.0), (0.5, 0.995, 3.0), (0.55, 0.996, 3.0), (0.6, 0.997, 3.0), (0.65, 0.997, 3.0), (0.7, 0.998, 3.0), (0.75, 0.995, 3.0), (0.8, 0.998, 3.0), (0.9, 0.982, 3.0), (1.0, 0.994, 3.0), (1.1, 0.995, 3.0), (1.2, 0.66, 3.0), (1.3, 0.97, 3.0), (1.4, 0.797, 3.0), (1.5, 0.904, 3.0), (1.6, 0.943, 3.0)], 1.5253, -1.0, -1.0, 0.0, 55.9, 0.0, 1.01, 0.0)
const _330R = Glass(GlassID(AGF, 1778), 1, 2.24171666, -0.00174741324, 0.011152331, 0.000782763713, -7.9413741e-5, 3.93670794e-6, NaN, NaN, NaN, NaN, 0.4, 1.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.0, -1.0, -1.0, 90.0, -1.0, 0.0, -1.0, [(0.4, 0.997, 3.0), (0.42, 0.998, 3.0), (0.44, 0.999, 3.0), (0.46, 0.999, 3.0), (0.48, 0.999, 3.0), (0.5, 1.0, 3.0), (0.55, 1.0, 3.0), (0.6, 1.0, 3.0), (0.65, 1.0, 3.0), (0.7, 1.0, 3.0), (0.75, 0.999, 3.0), (0.8, 1.0, 3.0), (0.9, 0.987, 3.0), (1.0, 0.996, 3.0), (1.1, 0.999, 3.0), (1.2, 0.55, 3.0), (1.3, 0.974, 3.0), (1.4, 0.764, 3.0), (1.5, 0.908, 3.0), (1.6, 0.941, 3.0)], 1.50942, -1.0, -1.0, 0.0, 56.0, 0.0, 0.95, 0.0)
end # module

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat


text:
using Plots
using MLFlowClient
using Random

# Parameters
αs = [0.0, 0.9, 0.98]
n = 100

"""Method that generates price paths of length `n` based on `α`"""
function getpricepath(α, n)
    x = zeros(n + 1)
    x[1] = 0.0
    for t in 1:n
        x[t+1] = α * x[t] + rand()
    end
    x
end
p = plot()

# Create MLFlow instance
mlf = MLFlow("http://localhost:5000")

# Initiate new experiment
experiment_id = createexperiment(mlf; name="price-paths")

# Create a run in the new experiment
exprun = createrun(mlf, experiment_id)

# Log parameters and their values
for (idx, α) in enumerate(αs)
    logparam(mlf, exprun, "alpha$(idx)", string(α)) # MLFlow only supports string parameter values
end

# Obtain pricepaths
pricepaths = [getpricepath(α, n) for α in αs]

# Log pricepaths in MLFlow
for (idx, pricepath) in enumerate(pricepaths)
    plot!(p,
        pricepath,
        title="Random price paths",
        label="alpha = $(αs[idx])",
        xlabel="Timestep",
        ylabel="Price"
    )

    logmetric(mlf, exprun, "pricepath$(idx)", pricepath)
end

# Save the price path plot as an image
plotfilename = "pricepaths-plot.png"
png(plotfilename)

# Upload the plot as an artifact associated with the MLFlow experiment's run
logartifact(mlf, exprun, plotfilename)

# remote temporary plot which was already uploaded in MLFlow
rm(plotfilename)

# complete the experiment
updaterun(mlf, exprun, "FINISHED")

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
function test_EuclideanKernel_1()

    k = EuclideanKernel()
    p = rand()
    q = randn(10)
    x = k(p, q)
    x_correct = p * q
    return x ≈ x_correct

end

function test_EuclideanKernel_2()

    k = EuclideanKernel()
    p = rand()
    q = randn()
    x = k(p, q)
    x_correct = p * q
    return x ≈ x_correct

end

function test_GaussianKernel_1()

    k = GaussianKernel(2)
    p = rand()
    q = randn(10)
    return k(p, q) ≈ k(q, p)

end

function test_GaussianKernel_2()

    Λ = 2
    k = GaussianKernel(Λ)
    p = rand()
    q = randn()
    x = k(p, q)
    x_correct = exp(-0.5 * abs2((p - q) / Λ))
    return x ≈ x_correct

end

@testset "Kernels" begin

    @test test_EuclideanKernel_1()
    @test test_EuclideanKernel_2()
    @test test_GaussianKernel_1()
    @test test_GaussianKernel_2()

end


text:
# Base on the TypedPolynomials/abstract/substition.jl written by Robin Deits
export subs

# TODO Vararg{<:...} -> Vararg{...}
const Substitution = Pair{<:AbstractVariable}
const MultiSubstitution{N} = Pair{<:Tuple{Vararg{AbstractVariable, N}}, <:Tuple{Vararg{<:Any, N}}}
const MultiVectorSubstitution = Pair{<:Tuple{Vararg{AbstractVariable}}, <:AbstractVector}

# When the variables are promoted to be in the same vector they could be promoted into a monomial
const VectorMultiSubstitution = Pair{<:AbstractVector{<:AbstractMonomialLike}, <:Tuple}
const VectorMultiVectorSubstitution = Pair{<:AbstractVector{<:AbstractMonomialLike}, <:AbstractVector}
_monov2vart(s::Pair{<:AbstractVector{<:AbstractMonomial}}) = variable.(Tuple(s.first)) => s.second
_monov2vart(s) = s

const AbstractMultiSubstitution = Union{MultiSubstitution, MultiVectorSubstitution, VectorMultiVectorSubstitution, VectorMultiSubstitution}
const AbstractSubstitution = Union{Substitution, AbstractMultiSubstitution}
const Substitutions = Tuple{Vararg{AbstractSubstitution}}

abstract type AbstractSubstitutionType end
struct Subs <: AbstractSubstitutionType end
struct Eval <: AbstractSubstitutionType end
const AST = AbstractSubstitutionType

"""
    subs(polynomial, (x, y)=>(1, 2))

is equivalent to:

    subs(polynomial, (x=>1, y=>2))
"""
substitute(st::AST, p::APL, s::AbstractMultiSubstitution) = substitute(st, p, pairzip(_monov2vart(s)))

# Evaluate the stream
# If I do s2..., then
# subs(x, x=>x+y, y=>2) would call substitute(Subs(), x+y, y=>2)
# subs(x, x=>1, y=>2) would call substitute(Subs(), 1, y=>2)
# so it would force use to also define
# sustitute(::AST, ..., ::AbstractSubstitution...) for Any, APL and RationalPoly.
#substitute(st::AST, p::APL, s1::AbstractSubstitution, s2::AbstractSubstitution...) = substitute(st, substitute(st, p, s1), s2...)
substitute(st::AST, p::APL, s1::AbstractSubstitution, s2::AbstractSubstitution...) = substitute(st, substitute(st, p, s1), s2)

## Variables
substitute(st::AST, v::AbstractVariable, s::Substitutions) = substitute(st, v, s...)

## Monomials
powersubstitute(st::AST, s::Substitutions, p::Tuple{AbstractVariable, Integer}) = substitute(st, p[1], s...)^p[2]
powersubstitute(st::AST, s::Substitutions, p::Tuple{AbstractVariable, Integer}, p2...) = powersubstitute(st, s, p) * powersubstitute(st, s, p2...)
substitute(st::AST, m::AbstractMonomial, s::Substitutions) = powersubstitute(st, s, powers(m)...)

## Terms
substitute(st::AST, t::AbstractTerm, s::Substitutions) = coefficient(t) * substitute(st, monomial(t), s)

## Polynomials
_polynomial(α) = α
_polynomial(p::APL) = polynomial(p)
function substitute(st::AST, p::AbstractPolynomial, s::Substitutions)
    if iszero(p)
        _polynomial(substitute(st, zeroterm(p), s))
    else
        ts = terms(p)
        r1 = substitute(st, ts[1], s)
        R = Base.promote_op(+, typeof(r1), typeof(r1))
        result::R = convert(R, r1)
        for i in 2:length(ts)
            result += substitute(st, ts[i], s)
        end
        result
    end
end

## Fallbacks
substitute(st::AST, p::APL, s::Substitutions) = substitute(st, polynomial(p), s)
substitute(st::AST, q::RationalPoly, s::Substitutions) = substitute(st, q.num, s) / substitute(st, q.den, s)

# subs(x, x=>x+y, y=>2) would call substitute(Subs(), x+y, y=>2)
#substitute(st::AST, p::Union{APL, RationalPoly}, s::AbstractSubstitution...) = substitute(st, p, s)

# Substitute Arrays
function substitute(st::AST, A::AbstractArray{<:APL}, s::Substitutions)
    map(p -> substitute(st, p, s), A)
end
## Everything else
substitute(::AST, x, s::Substitutions) = x
# subs(x, x=>1, y=>2) would call substitute(Subs(), 1, y=>2)
#substitute(::AST, x, s::AbstractSubstitution...) = x

"""
    subs(p, s::AbstractSubstitution...)

Apply the substitutions `s` to `p`.
Use `p(s...)` if we are sure that all the variables are substited in `s`.

The allowed substutions are:

* `v => p` where `v` is a variable and `p` a polynomial, e.g. `x => 1` or `x => x^2*y + x + y`.
* `V => P` where `V` is a tuple or vector of variables and `P` a tuple or vector of polynomials, e.g. `(x, y) => (y, x)` or `(y, x) => (2, 1)`.

The order of the variables is lexicographic with the name with TypedPolynomials and by order of creation with DynamicPolynomials.
Since there is no guarantee on the order of the variables, substitution directly with a tuple or a vetor is not allowed.
You can use `p(variables(p) => (1, 2))` instead if you are sure of the order of the variables (e.g. the name order matches the creation order).

### Examples

```julia
p = 3x^2*y + x + 2y + 1
p(x => 2, y => 1) # Return type is Int
subs(p, x => 2, y => 1) # Return type is Int in TypedPolynomials but is a polynomial of Int coefficients in DynamicPolynomials
subs(p, y => x*y^2 + 1)
p(y => 2) # Do not do that, this works fine with TypedPolynomials but it will not return a correct result with DynamicPolynomials since it thinks that the return type is `Int`.
```

"""
subs(p, s::AbstractSubstitution...) = substitute(Subs(), p, s)

(p::RationalPoly)(s::AbstractSubstitution...) = p.num(s...) / p.den(s...)


text:
module TestOperators

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running operator tests")

## init test values
vals = [1 1 1 1;
        2 2 2 2;
        3 3 3 3;
        4 4 4 4]
dats = Date[Date(2013, 7, ii) for ii=1:4]
nams = [:A, :B, :C, :D]
valsDf = TimeData.composeDataFrame(vals, nams)

###############
## operators ##
###############

## almost all tests for Timematr should also work for Timenum without
## missing observations

## separate testing
## - matrix-mult
## - comparisons for Timedata objects
macro test_basic_operators(myType)
    esc(quote

        td = $(myType)(vals, nams, dats)

        ## mathematical operators
        td[2:end, :] .- td[1:(end-1), :]

        @test isequal(td .+ td, 2.*td)
        @test isequal(td, 3.*(td./3))
        @test isequal(get((3./td).*td, 1, 1), 3)

        ## arithmetics with scalar values
        @test isequal(td .+ 1,  1 .+ td)

        @test isequal((td .- 1.5).*(-1), 1.5 .- td)
        @test isequal(td .* 3, 3 .* td)

        ## unary operators / mathematical functions
        pres_msUnitary_functions = [:abs, :sign,
                                    :exp, :log, :sqrt, :gamma,
                                    :(+), :(-)
                                    ]
        for f in pres_msUnitary_functions
            eval(:($(f)(td)))
        end

        ## comparisons
        td .> td
        td .< td
        td .<= td
        td .>= td
        td .== td
        td .!= td

        td .< 0.5
        td .== 0.3
        td .!= 0.3

        kk = (td .> 0.5) | (td .> 0.7)
        kk2 = td .> 0.5
        @test isequal(kk, kk2)

        td1 = exp(td)
        td2 = log(td1)
        @test_approx_eq convert(Array{Float64, 1}, get(td)[:]) convert(Array{Float64, 1}, get(td2)[:])

        ## rounding functions
        roundFunc = [:round, :ceil, :floor, :trunc]
        for f in roundFunc
            eval(:($(f)(td)))
            eval(:($(f)(td, 2)))
        end
    end)
end

for t = (:Timenum, :Timematr)
    eval(macroexpand(:(@test_basic_operators($t))))
end

############
## test ! ##
############

df = DataFrame(a = @data([NA, true, false]),
               b = @data([true, NA, false]))
dfNegated = DataFrame(a = @data([NA, false, true]),
               b = @data([false, NA, true]))

td = Timedata(df)
@test isequal(!td, Timedata(dfNegated))

end


text:
primes = [2]

function updateprimes(n)
  global primes
  isprime = true
  for p in primes
    if n % p == 0
      isprime = false
      break
    end
  end
  if isprime
    push!(primes, n)
  end
end

function primefactorize(n)
  global primes
  pf = repeat([0], n)
  for p in primes
    m = n
    while m % p == 0
      pf[p] += 1
      m /= p
    end
  end
  return pf
end

n, m = 2, 3
pfn = primefactorize(n)
dv = 2
while dv <= 500
 global n, m, pfn, pfm, dv
 updateprimes(m)
 pfm = primefactorize(m)
 if m % 2 == 0
  pfm[2] -= 1
 end
 push!(pfn, 0)
 dv = prod(map(x -> x + 1, filter(x -> x > 0, (pfn + pfm))))
 n = m
 pfn = pfm
 m += 1
end
t = Int((n-1) * n / 2)
println(t)

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
include("Elasticity_110.jl")
using JuLIPMaterials.CLE.Elasticity_110:
      elastic_moduli,
      voigt_moduli,
      fourth_order_basis,
      sextic_roots,
      A_coefficients,
      D_coefficients,
      little_a

using JuLIP.Potentials: fcut, fcut_d

u_edge_fcc_110(x, y, b, C::Array{T,4}; TOL=1e-4) where {T} =
         u_edge(x, y, b, voigt_moduli(C), TOL=TOL)

"""
* `u_edge{T}(x, y, b, C::Array{T,N}) -> ux, uy`

* `C` can be either 3 x 3 x 3 x 3 (`elastic_moduli`) or 6 x 6 (`voigt_moduli`)

This function computes the anisotropic CLE solution for an in-plane
edge dislocation. The elastic moduli are taken to within `TOL` accuracy (
   this is a keyword argument)
"""
function u_edge_fcc_110(x, y, b, Cv::Array{T,2}, a; TOL = 1e-4) where {T}
   Cv = copy(Cv)
   test1 = Cv[1,1]
   test2 = Cv[1,2]
   test3 = Cv[6,6]
   # >>>>>>>>> START DEBUG >>>>>>>>
   #Zero out elasticity tensor from any numerical artifacts
   Cv = zeros(6,6)
   Cv[1,1] = 1.0*test1
   Cv[1,2] = 1.0*test2
   Cv[6,6] = 1.0*test3
   Cv[2,2] = Cv[3,3] = Cv[1,1]
   Cv[1,3] = Cv[2,3] = Cv[2,1] = Cv[3,1] = Cv[3,2] = Cv[1,2]
   Cv[4,4] = Cv[5,5] = Cv[6,6]
   # <<<<<<<<< END DEBUG <<<<<<<<<

   #Compute Anisotropic solution from Chou and Sha, J. App. Phys 42 (7) 2625
   #This uses the elasticity tensor in the usual coordinates.
   #Note that the rotated tensor agrees with the values shown on 2625
   #Same problem as with Hirth and Lothe--either formulas are incorrect or problem
   #with numeric instability.  Complex log formulas seem to work
   bar = sqrt( .5*Cv[1,1]*(Cv[1,1] +Cv[1,2] + 2.0*Cv[4,4]) )
   C = ( (bar + Cv[1,2])*(bar - Cv[1,2] - 2.0*Cv[4,4] )  )/(bar*Cv[4,4])
   lem = sqrt(bar/Cv[1,1]+0im)

   delta2 = sqrt(-C)
   delta1 = sqrt(C+4)

   #Now compute using Hirth and Lothe which should be valid for 110 dislocation
   #This should use K instead of Cv
   #Something seems to be either wrong with these formulas or numeric issues are arising
   #Skip to general formula using complex logs

   K = fourth_order_basis(Cv,a)
   Cvoigt = round.(K, digits=3)
   c̄11 = sqrt(K[1,1]*K[2,2])    # (13-106)
   lam = (K[1,1]/K[2,2])^(1/4)
     ϕ = 0.5 * acos( (K[1,2]^2 + 2*K[1,2]*K[6,6] - c̄11^2) / (2.0*c̄11*K[6,6]) )
   apple = - lam*(K[6,6]*exp(im*ϕ) + c̄11*exp(-im*ϕ))/(K[1,2]+K[6,6])
   dodo = im*b[1,1]/(2.0*c̄11*sin(2*ϕ))*(K[1,2] - c̄11*exp(2*im*ϕ))

   #Should test this against solving the full linear system
   A = Complex{Float64}[0 0 0; 0 0 0; 0 0 0]
   A = [0 1 1; 0 apple -apple; 0 0 0]
   #Set up for burgers vector in x1 direction only
   D = zeros(6,1)
   D = [ 0; 0; real(dodo); imag(dodo); -real(dodo); -imag(dodo)]
   p = sextic_roots(K)
   #Do some funny business with root switching here.  This could be what is screwing things up.
   #May have to change the roots elsewhere instead
   p[3,1] = -p[2,1]

   ux = real.( im/(2*π)*(A[1,2]*(D[3] + D[4]*im)*log.(x+p[2]*y) + A[1,3]*(D[5] + D[6]*im)*log.(x+p[3]*y)   ))

   uy = real.( im/(2*π)*(A[2,2]*(D[3] + D[4]*im)*log.(x+p[2]*y) + A[2,3]*(D[5] + D[6]*im)*log.(x+p[3]*y)   ))

   return ux, uy
end

mutable struct EdgeCubic{T1,T2,T3,T4,T5}
    A::Matrix{T1}
    D::Vector{T2}
    p::Vector{T3}
    x0::T4
    b::T5
end

function EdgeCubic(b::Real, Cv::Array{T,2}, a::Real; x0 = zeros(3)) where {T}
    # clean up the tensor
   test1 = Cv[1,1]
   test2 = Cv[1,2]
   test3 = Cv[6,6]
   Cv = zeros(6,6)
   Cv[1,1] = test1
   Cv[1,2] = test2
   Cv[6,6] = test3
   #Zero out elasticity tensor from any numerical artifacts
   Cv[2,2] = Cv[3,3] = Cv[1,1]
   Cv[1,3] = Cv[2,3] = Cv[2,1] = Cv[3,1] = Cv[3,2] = Cv[1,2]
   Cv[4,4] = Cv[5,5] = Cv[6,6]

   # parameters needed for the Hirth/Lothe solution
   bar = sqrt( .5*Cv[1,1]*(Cv[1,1] +Cv[1,2] + 2.0*Cv[4,4]) )
   C = ( (bar + Cv[1,2])*(bar - Cv[1,2] - 2.0*Cv[4,4] )  )/(bar*Cv[4,4])
   lem = (bar/Cv[1,1])^(1/2)
   delta2 = sqrt(-C)
   delta1 = sqrt(Complex(C+4))

   # Now compute using Hirth and Lothe which should be valid for 110 dislocation
   # This should use K instead of Cv
   # Something seems to be either wrong with these formulas or numeric issues are arising
   # Skip to general formula using complex logs
   K = fourth_order_basis(Cv,a; Tr = [0 1/√2  -1/√2; 1 0 0; 0 1/√2 1/√2])
   Cvoigt = round.(K, 8)

   # more parameters
   c̄11 = sqrt(K[1,1]*K[2,2])    # (13-106)
   lam = (K[1,1]/K[2,2])^(1/4)
   ϕ = 0.5 * acos( (K[1,2]^2 + 2*K[1,2]*K[6,6] - c̄11^2) / (2.0*c̄11*K[6,6]) )
   apple = - lam*(K[6,6]*exp(im*ϕ) + c̄11*exp(-im*ϕ))/(K[1,2]+K[6,6])
   dodo = im*b[1,1]/(2.0*c̄11*sin(2*ϕ))*(K[1,2] - c̄11*exp(2*im*ϕ))

   #Should test this against solving the full linear system
   A = Complex{Float64}[0 1 1; 0 apple -apple; 0 0 0]
   #Set up for burgers vector in x1 direction only
   D = zeros(6,1)
   D = [ 0; 0; real(dodo); imag(dodo); -real(dodo); -imag(dodo)]
   p = sextic_roots(K)
   #Do some funny business with root switching here.  This could be what is screwing things up.
   #May have to change the roots elsewhere instead
   p[3,1] = -p[2,1]

   return EdgeCubic(A, D, p, x0[1:2], b)
end

#Create function to smoothly blend from 0 to 1
eta(Y::AbstractVector, r0=0.1, r1=1.8) = 1.0 - fcut(norm(Y), r0, r1)

grad_eta(Y::AbstractVector, r0=0.1, r1=1.8) = - fcut_d(norm(Y), r0, r1) * (Y/norm(Y))

#xi correction solver from Ehrlacher, Ortner, Shapeev
#
# TODO: we could make this robust
#
function xi_solver(Y::Vector, b::Float64; TOL = 1e-9, maxnit = 500)
    ξ1(x::Real, y::Real, b) = x - b * eta(Y) * angle(x + im * y) / (2*π)
    dξ1(x::Real, y::Real, b) = 1 + b * eta(Y) * y / (x^2 + y^2) / (2*π) + b * grad_eta(Y)[1] * angle(x + im * y)/ (2*π)
    y = Y[2]
    x = y
    for n = 1:maxnit
        f = ξ1(x, y, b) - Y[1]
        if abs(f) <= TOL; break; end
        x = x - f / dξ1(x, y, b)
    end
    if abs(ξ1(x, y, b) - Y[1]) > TOL
        warn("newton solver did not converge; returning input")
        return Y
    end
    return [x, y]
end

#return the inverse of the xi derivative
function xi_deriv_inv(Y::Vector, b::Float64)
    dξ1(x::Real, y::Real, b) = 1 + b * y / (x^2 + y^2) / (2*π) + b * grad_eta(Y)[1] * angle(x + im * y)/ (2*π)
    dξ2(x::Real, y::Real, b) = -b * x / (x^2 + y^2) / (2*π) + b * grad_eta(Y)[2] * angle(x + im * y)/ (2*π)
    D_xi_inv = (1/dξ1(Y[1],Y[2], b))*[1 -dξ2(Y[1],Y[2], b) 0;0 dξ1(Y[1],Y[2], b) 0; 0 0 0]
    return D_xi_inv
end

function evaluate(U::EdgeCubic, X::AbstractVector)
    x = X[1] - U.x0[1]
    y = X[2] - U.x0[2]
    x, y = xi_solver([x'; y'], -U.b)
    A, D, p = U.A, U.D, U.p
    ux = real( im/(2*π)*(A[1,2]*(D[3] + D[4]*im)*log(x+p[2]*y) + A[1,3]*(D[5] + D[6]*im)*log(x+p[3]*y) ))
    uy = real( im/(2*π)*(A[2,2]*(D[3] + D[4]*im)*log(x+p[2]*y) + A[2,3]*(D[5] + D[6]*im)*log(x+p[3]*y) ))
    return [ux, uy, 0.0]
end

function jacobian(U::EdgeCubic, X::AbstractVector)
    x = X[1] - U.x0[1]
    y = X[2] - U.x0[2]
    x, y = xi_solver([x'; y'], -U.b)
    D_xi_inv = xi_deriv_inv([x'; y'], -U.b)
    A, D, p = U.A, U.D, U.p
    uxx = real( im/(2*π)*(A[1,2]*(D[3] + D[4]*im)/(x+p[2]*y) + A[1,3]*(D[5] + D[6]*im)/(x+p[3]*y) ) )
    uxy = real( im/(2*π)*(A[1,2]*(D[3] + D[4]*im)*p[2]/(x+p[2]*y) + A[1,3]*(D[5] + D[6]*im)*p[3]/(x+p[3]*y) ) )
    uyx = real( im/(2*π)*(A[2,2]*(D[3] + D[4]*im)/(x+p[2]*y) + A[2,3]*(D[5] + D[6]*im)/(x+p[3]*y)   ))
    uyy = real( im/(2*π)*(A[2,2]*(D[3] + D[4]*im)*p[2]/(x+p[2]*y) + A[2,3]*(D[5] + D[6]*im)*p[3]/(x+p[3]*y)   ))
    return [uxx uxy 0.0; uyx uyy 0.0; 0.0 0.0 0.0]*D_xi_inv
end

function (U::EdgeCubic)(X::AbstractVector)
    return evaluate(U, X), jacobian(U, X)
end


text:
using ITensors
using ITensorsGateEvolution

function main()

N = 8

osX = ProductOps()
for n in 1:N
  osX *= "X", n
end

osZ = ProductOps()
for n in 1:N
  osZ *= "Z", n
end

osSw = ProductOps()
for n in 1:N-2
  osSw *= "Sw", n, n+2
end

osCx = ProductOps()
for n in 1:N-3
  osCx *= "Cx", n, n+3
end

osT = ProductOps()
for n in 1:N-3
  osT *= "T", n, n+1, n+3
end

osRx = ProductOps()
for n in 1:N
  osRx *= "Rx", n, (θ = π,)
end

os_noise = ProductOps()
for n in 1:N-4
  os_noise *= "noise", n, n+2, n+4
end

os = osX * osSw * osRx * osZ * osCx * osT * os_noise

#os = ProductOps()
#os *= "T", 1, 3, 4

@show os

s = siteinds("qubit", N; conserve_qns = false)
gates = ops(s, os)

M0 = MPO(s, "Id")

# Apply the gates

s0 = siteinds(M0)

M = apply(gates, M0; apply_dag = true, cutoff = 1e-15, maxdim = 500)
@show dim(s[1])^(N ÷ 2)
@show maxlinkdim(M)

s = siteinds(M)
for n in 1:N
  @assert hassameinds(s[n], s0[n])
end

println("Evolution complete, test the result")

set_warn_itensor_order!(18)

prodM = apply(gates, prod(M0); apply_dag = true)

@show prod(M) ≈ prodM
@show norm(prod(M) - prodM) / norm(prodM)

reset_warn_itensor_order!()

return M

end

M = main();


text:
## -- Changepoint.jl tests

    A = [randn(100).-2; randn(100).+2]
    nsteps = 10000
    burnin = 1000

    dist = changepoint(A, 10000; np=1)
    @test isapprox(nanmean(dist[burnin:end]), 101, atol=4)
    dist = changepoint(A, 10000; npmin=1, npmax=5)[burnin:end,:]
    @test isapprox(nanmean(dist[dist.>0]), 101, atol=4)

    dist = changepoint(A, ones(200), 10000; np=1)
    @test isapprox(nanmean(dist[burnin:end]), 101, atol=4)
    dist = changepoint(A, ones(200), 10000; npmin=1, npmax=5)[burnin:end,:]
    @test isapprox(nanmean(dist[dist.>0]), 101, atol=4)

## --


text:
import Base: mapreduce, sum
import Statistics: mean

@eval begin
    @explicit_intercepts(
        mapreduce,
        Tuple{Function, $plustype, AbstractArray{<:∇Scalar}},
        [false, false, true],
        (dims=:,),
    )
    function ∇(
        ::typeof(mapreduce),
        ::Type{Arg{3}},
        p, y, ȳ, f,
        ::$plustype,
        A::AbstractArray{<:∇Scalar};
        dims=:,
    )
        hasmethod(∇, Tuple{typeof(f), Type{Arg{1}}, ∇Scalar}) ?
            broadcast((An, ȳn)->ȳn * ∇(f, Arg{1}, An), A, ȳ) :
            broadcast((An, ȳn)->ȳn * fmad(f, (An,), Val{1}), A, ȳ)
    end
end

@explicit_intercepts(
    sum,
    Tuple{Function, AbstractArray{<:∇Scalar}},
    [false, true],
    (dims=:,),
)
@explicit_intercepts(
    sum,
    Tuple{AbstractArray{<:∇Scalar}},
    [true],
    (dims=:,),
)
function ∇(
    ::typeof(sum),
    ::Type{Arg{2}},
    p, y, ȳ,
    f::Function,
    A::AbstractArray{<:∇Scalar};
    dims=:,
)
    # Just pass through to mapreduce
    return ∇(mapreduce, Arg{3}, p, y, ȳ, f, Base.add_sum, A; dims=dims)
end
function ∇(
    ::typeof(sum),
    ::Type{Arg{1}},
    p, y, ȳ,
    A::AbstractArray{<:∇Scalar};
    dims=:,
)
    # Again pass through to mapreduce, using identity as the mapped function
    return ∇(mapreduce, Arg{3}, p, y, ȳ, identity, Base.add_sum, A; dims=dims)
end
# Specialize on sum(abs2, x) as it's a common pattern with a simple derivative
function ∇(
    ::typeof(sum),
    ::Type{Arg{2}},
    p, y, ȳ,
    ::typeof(abs2),
    A::AbstractArray{<:∇Scalar};
    dims=:,
)
    return 2ȳ .* A
end

@explicit_intercepts(
    mean,
    Tuple{Function, AbstractArray{<:∇Scalar}},
    [false, true],
    #(dims=:,)  # https://github.com/JuliaLang/julia/issues/31412
)
@explicit_intercepts(
    mean,
    Tuple{AbstractArray{<:∇Scalar}},
    [true],
    (dims=:,)
)

_denom(x, dims::Colon) = length(x)
_denom(x, dims::Integer) = size(x, dims)
_denom(x, dims) = mapreduce(i->size(x, i), Base.mul_prod, unique(dims), init=1)

function ∇(
    ::typeof(mean),
    ::Type{Arg{2}},
    p, y, ȳ,
    f::Function,
    x::AbstractArray{<:∇Scalar},
)
    return ∇(sum, Arg{2}, p, y, ȳ, f, x; dims=:) / _denom(x, :)
end
function ∇(
    ::typeof(mean),
    ::Type{Arg{1}},
    p, y, ȳ,
    x::AbstractArray{<:∇Scalar};
    dims=:,
)
    return ∇(sum, Arg{1}, p, y, ȳ, x; dims=dims) ./ _denom(x, dims)
end


text:
############################################################################################
# Check types for Markov Kernel
markov_pf = ParticleFilter(_rng, markov_objective)
_vals, _diag = propose(_rng, markov_pf, markov_objective)
@testset "PF - Types - Markov" begin
    @test eltype(markov_latent) == eltype(markov_pf.particles.val)
    @test typeof(markov_latent) == typeof(markov_pf.particles.buffer.parameter.val)
    @test typeof(markov_latent) == typeof(markov_pf.particles.buffer.proposal)
    ## Check predictions
    #@test typeof(_diag.prediction) == Tuple{eltype(markov_latent),eltype(markov_data)}
end

############################################################################################
# Check types for SemiMarkov Kernel
semimarkov_pf = ParticleFilter(_rng, semimarkov_objective)
_vals, _diag = propose(_rng, semimarkov_pf, semimarkov_objective)
@testset "PF - Types - SemiMarkov" begin
    @test eltype(semimarkov_latent) == eltype(semimarkov_pf.particles.val)
    @test typeof(semimarkov_latent) == typeof(semimarkov_pf.particles.buffer.parameter.val)
    @test typeof(semimarkov_latent) == typeof(semimarkov_pf.particles.buffer.proposal)
    ## Check predictions
    #@test typeof(_diag.prediction) == Tuple{eltype(semimarkov_latent),eltype(semimarkov_data)}
end