text:
using TimetableSolver
using JuMP
using ConstraintSolver
using DataStructures: OrderedDict

# TODO: add return type for functions

const CS = ConstraintSolver

"""
    VarRefOrExpr

Represents either a model variable reference or expression reference.

Alias for union of VariableRef and AffExpr.
"""
VarRefOrExpr = Union{VariableRef,AffExpr}

#=
Each variable is a string of the form:
    `{div}_{row}_{period}_{type}`,
where div is the str representation of division,
row is the number of row, period is the number of period,
and type is either 'subject' or 'teacher'.

Each domain value represents a subject or teacher, so its value
should be a str of the form: `{subject}` or `{teacher}`.

However, all the strings are hidden and ints are used for solving directly.
Mappings are generated and stored in VariableData struct.
=#

"""
    VariableData

Immutable struct to store data related to divisions, variables and values,
and their mappings from ints to strings, because the solver only works on ints.

Fields:
- `divvarval_maps::OrderedDict{Symbol,OrderedDict}`: Keys are a symbol `:subject` or `:teacher`.
  Values are map of ints (div) to 2-tuple of vecs of ints (var, val).
- `divvarval_mapstrs::OrderedDict{Symbol,OrderedDict}`: Same as `div_var_val_map` but with strings instead of ints.

- `alldivs::Vector{Int}`: Vector of all divs as ints.
- `allvars::Vector{Int}`: Vector of all vars as ints. (no `:subject` `:teacher` bifurcation)
- `allvals::Vector{Int}`: Vector of all vals as ints. (no `:subject` `:teacher` bifurcation)

- `divmap::OrderedDict{Int,String}`: Map from div ints to strings.
- `varmap::OrderedDict{Int,String}`: Map from variable ints to strings.
- `valmap::OrderedDict{Int,String}`: Map from values ints to strings.
- `inversedivmap::OrderedDict{String,Int}`: Inverted `divmap`.
- `inversevarmap::OrderedDict{String,Int}`: Inverted `varmap`.
- `inversevalmap::OrderedDict{String,Int}`: Inverted `valmap`.

- `schedule::Schedule`: Schedule type instance.

"""
struct VariableData
    divvarval_maps::OrderedDict{Symbol,OrderedDict}
    divvarval_mapstrs::OrderedDict{Symbol,OrderedDict}

    alldivs::Vector{Int}
    allvars::Vector{Int}
    allvals::Vector{Int}

    divmap::OrderedDict{Int,String}
    varmap::OrderedDict{Int,String}
    valmap::OrderedDict{Int,String}
    inversedivmap::OrderedDict{String,Int}
    inversevarmap::OrderedDict{String,Int}
    inversevalmap::OrderedDict{String,Int}

    schedule::Schedule
end

"""
    VariableData(schedule::Schedule)::VariableData

Return `VariableData` instance with all the fields filled based on the schedule.
"""
function VariableData(schedule::Schedule)::VariableData
    # get divvarval_mapstrs
    divvarval_mapstrs = get_divvarval_mapstrs(schedule)

    # extract str representations of all data
    alldiv_strs = Vector{String}([])
    allvar_strs = OrderedDict{Symbol,Vector{String}}(
        :subject => [],
        :teacher => []
    )
    allval_strs = OrderedDict{Symbol,Vector{String}}(
        :subject => [],
        :teacher => []
    )
    # symbol is `:subject` or `:teacher`
    # symbol_map is div_var_val_map_strs but for either subject or teacher only
    for (symbol, symbolmap) in divvarval_mapstrs
        for (div_str, varval_strs) in symbolmap
            push!(alldiv_strs, div_str)
            for var_str in varval_strs[1]
                # push the var to the proper key in all_var_strs
                push!(allvar_strs[symbol], var_str)
            end
            for val_str in varval_strs[2]
                # push the val to the proper key in all_var_strs
                push!(allval_strs[symbol], val_str)
            end
        end
    end

    # remove duplicate values from allval_strs
    allval_strs[:subject] = unique(allval_strs[:subject])
    allval_strs[:teacher] = unique(allval_strs[:teacher])

    # remove duplicate values from alldiv_strs
    alldiv_strs = unique(alldiv_strs)

    # define int representations of all data
    # and int to str mappings
    alldivs = Vector{Int}(1:length(alldiv_strs))
    divmap::OrderedDict{Int,String} = OrderedDict(
        enumerate(alldiv_strs)
    )
    inversedivmap::OrderedDict{String,Int} = OrderedDict(v => k for (k, v) in divmap)

    allvars = Vector{Int}(1:(length(allvar_strs[:subject]))+length(allvar_strs[:teacher]))
    # nested list comp
    # return var_str for var_str in symbol_var_strs which are in all_var_strs
    varmap::OrderedDict{Int,String} = OrderedDict(
        enumerate([v for (_, var_strs) in allvar_strs for v in var_strs])
    )
    inversevarmap::OrderedDict{String,Int} = OrderedDict(v => k for (k, v) in varmap)

    allvals = Vector{Int}(1:(length(allval_strs[:subject]))+length(allval_strs[:teacher]))
    # nested list comp (same as for `var_map`)
    valmap::OrderedDict{Int,String} = OrderedDict(
        enumerate([v for (_, val_strs) in allval_strs for v in val_strs])
    )
    inversevalmap::OrderedDict{String,Int} = OrderedDict(v => k for (k, v) in valmap)

    # dict to store ordereddict of all vars and vals mapped to divs
    # dict keys are :subject and :teacher
    divvarval_maps::OrderedDict{Symbol,OrderedDict} = OrderedDict(
        :subject => OrderedDict(),
        :teacher => OrderedDict()
    )

    # generate values for div_var_val_map
    # symbol is `:subject` or `:teacher`
    # symbol_map is divvarval_mapstrs but for either subject or teacher only
    for (symbol, symbolmap) in divvarval_mapstrs
        for (i, (div_str, (var_strs, val_strs))) in enumerate(symbolmap)
            # int representing current div
            div::Int = inversedivmap[div_str]

            # ints representing current vars
            vars = Vector{Int}([inversevarmap[v] for v in var_strs])
            vals = Vector{Int}([inversevalmap[v] for v in val_strs])

            # add key value pair to div_var_val_map
            divvarval_maps[symbol][div] = (vars, vals)
        end
    end

    return VariableData(
        divvarval_maps, divvarval_mapstrs,
        alldivs, allvars, allvals,
        divmap, varmap, valmap, inversedivmap, inversevarmap, inversevalmap,
        schedule
    )
end

"""
    get_divvarval_mapstrs(schedule::Schedule)::OrderedDict{Symbol,OrderedDict}

Return `OrderedDict` mapping symbol (either `:subject` or `:teacher`), to an `OrderedDict`.

The first-order values (`OrderedDict`), are mappings from division string, to a 2-tuple of variables and values.
"variables" and "values" are vectors of variable strings and value strings.
"""
function get_divvarval_mapstrs(schedule::Schedule)::OrderedDict{Symbol,OrderedDict}
    # dict to store ordereddict of all vars and vals mapped to divs
    # dict keys are :subject and :teacher
    divvarval_maps::OrderedDict{
        Symbol,OrderedDict{String,Tuple{Vector{String},Vector{String}}}
    } = OrderedDict(
        :subject => OrderedDict(),
        :teacher => OrderedDict()
    )

    # fill divvarval_maps with data
    for (div, tt) in schedule.data
        subjectvals::Vector{String} = []
        teachervals::Vector{String} = []

        # add dom vals for both teachers and subjects
        for (subject, teachers) in tt.subjectteachers
            push!(subjectvals, subject)
            for t in teachers
                push!(teachervals, t)
            end
        end

        subjectvars::Vector{String} = []
        teachervars::Vector{String} = []

        for (i, row) in enumerate(tt.data)
            for j in (1:length(row))
                # add subject variable
                subject = "$(div)_$(i)_$(j)_subject"
                push!(subjectvars, subject)
                # add teacher variable
                teacher = "$(div)_$(i)_$(j)_teacher"
                push!(teachervars, teacher)
            end
        end

        # only use unique vals
        subjectvals = unique(subjectvals)
        teachervals = unique(teachervals)

        # create a key for every div, in teacher and subject
        divvarval_maps[:subject][div] = (subjectvars, subjectvals)
        divvarval_maps[:teacher][div] = (teachervars, teachervals)
    end

    return divvarval_maps
end

# TODO: add kwargs arg to this and other function to pass kwargs to Model()
"""
    get_model(vardata::VariableData, timeout::Float64=Inf, all_solutions::Bool=false)::Tuple{Model,OrderedDict}

Return a JuMP model created using `vardata`, and an OrderedDict of the variables defined in the model.

# Arguments
- `vardata::VariableData`: VariableData object representing the variables and constraints to be added to the model.
- `timeout::Float64=Inf`: Float representing the maximum time to run the model.
- `all_solutions::Bool=false`: Bool representing whether to return all solutions or just the first.

# Notes
- Create a model with the correct variables, and adds constraints to the model.
- The variables defined are anonymous, so they are stored in an OrderedDict mapped to their int representation.
"""
function get_model(vardata::VariableData, timeout::Float64 = Inf, all_solutions::Bool = false
)::Tuple{Model,OrderedDict}
    # create a constraint solver model and set ConstraintSolver as the optimizer
    m = Model(optimizer_with_attributes(CS.Optimizer,
        # logging vector containing :Info and/or :Table or empty
        "logging" => [],

        # whether to return all possible solutions
        "all_solutions" => all_solutions,

        # seed for random operations given for reproducibility
        "seed" => 0,
        # default to :Auto
        # "traverse_strategy" => :BFS,
        # "traverse_strategy"=> :DFS,
        # "traverse_strategy"=> :DBFS,

        # default :Auto
        # "branch_split"=>:Smallest,
        # "branch_split"=>:Biggest,
        # "branch_split" => :InHalf,

        "branch_strategy" => :IMPS, # default

        # simplify alldifferent, etc. constraints
        "simplify" => true, # default

        # max seconds before exiting with status `MOI.TIME_LIMIT`
        "time_limit" => timeout,

        # backtracking options
        "backtrack" => true, # default
        "backtrack_sorting" => true, # default

        # "lp_optimizer" => cbc_optimizer,
        # "lp_optimizer" => glpk_optimizer,
        # "lp_optimizer" => ipopt_optimizer,
    ))

    # define variables
    modelvars = define_variables!(m, vardata)

    # define constraints
    define_constraints!(m, modelvars, vardata)

    return m, modelvars
end

"""
    define_variables!(m::Model, vardata::VariableData)::OrderedDict{Int,VarRefOrExpr}

Define variables in the model `m` using the `vardata` object.
Return an OrderedDict mapping the int representation of the variable to the variable reference/expression reference.

# Notes
- Modify the model `m` in place.
- See [`VarRefOrExpr`](@ref).
"""
function define_variables!(
    m::Model, vardata::VariableData)::OrderedDict{Int,VarRefOrExpr}
    divvarval_maps = vardata.divvarval_maps

    # mapping variable int repr to its model variableref
    modelvars = OrderedDict{Int,VarRefOrExpr}()

    for (symbol, symbolmap) in divvarval_maps
        for (div, (vars, vals)) in symbolmap
            # define a set of vars at a time
            currentmodelvars = @variable(m, [vars], CS.Integers(vals))
            # allocate them to the value of var in modelvars one by one
            for v in vars
                modelvars[v] = currentmodelvars[v]
            end
        end
    end

    return modelvars
end

"""
    define_subjectconstraints!(m::Model, modelvars::OrderedDict, vardata::VariableData)::Nothing

Define subject constraints in the model `m` using the `vardata` object and the model variables `modelvars`.

# Notes
- Modify the model `m` in place.
"""
function define_subjectconstraints!(m::Model, modelvars::OrderedDict, vardata::VariableData)::Nothing
    schedule = vardata.schedule
    divvarval_maps = vardata.divvarval_maps
    inversedivmap = vardata.inversedivmap
    inversevalmap = vardata.inversevalmap
    for (div_str, tt) in schedule.data
        # define terms
        div = inversedivmap[div_str]
        counts = OrderedDict{Int,Int}(
            inversevalmap[subject_str] => count
            for (subject_str, count) in tt.subjectcounts
        )
        vars = divvarval_maps[:subject][div][1]
        vals = collect(Int, keys(counts))
        # tl;dr (var::Int = 2) -> (var::Binary[] = [0, 1, 0]) index shows the value
        # create a matrix of binary variables:
        # `vars` as indices and `vals` as number of values
        # so instead of setting an int as value of each var
        # each var has `vals` corresponding binary variables
        binvars = @variable(m, [vars, vals], Bin)
        # add constraint to only choose one val
        @constraint(m, [i = vars], sum(binvars[i, :]) == 1)
        # add constraint to only choose that val `counts[val]` times
        @constraint(m, [j = vals], sum(binvars[:, j]) == counts[j])
        # array of expressions that actually store the value of var as int
        exprs = @expression(m, [i = vars], sum(j * binvars[i, j] for j in vals))
        # replace the corresponding vars in modelvars with exprs
        for var in vars
            # delete the var from the model
            delete(m, modelvars[var])
            # set the var to an expr defined above
            modelvars[var] = exprs[var]
        end
    end
end

"""
    define_subjectteacherconstraints!(m::Model, modelvars::OrderedDict, vardata::VariableData)::Nothing

Define subject-teacher constraints in the model `m` using the `vardata` object and the model variables `modelvars`.

# Notes
- Modify the model `m` in place.
"""
function define_subjectteacherconstraints!(m::Model, modelvars::OrderedDict, vardata::VariableData)::Nothing
    divvarval_maps = vardata.divvarval_maps
    inversevalmap = vardata.inversevalmap
    schedule = vardata.schedule

    # mapping div to corresponding subjectteachers
    subjectteacher_maps = OrderedDict{Int,OrderedDict}()
    for (div_str, tt) in schedule.data
        div = vardata.inversedivmap[div_str]
        # subject teacher mappings as ints
        subjectteachers = OrderedDict{Int,Vector{Int}}()
        # get the subjectteacher value for this timetable
        subjectteacher_strs = tt.subjectteachers
        # convert strs to ints, and add to `subjectteachers`
        for (subject, teachers) in subjectteacher_strs
            subject_int = inversevalmap[subject]
            teacher_ints = [inversevalmap[t] for t in teachers]
            subjectteachers[subject_int] = teacher_ints
        end
        subjectteacher_maps[div] = subjectteachers

    end

    # iterate over every set of subject variables (by div)
    for (div, (subjectvars, _)) in divvarval_maps[:subject]
        # generate pairs of all possible subject teacher combinations
        # and store them in a matrix
        numpairs = sum(length.(values(subjectteacher_maps[div])))
        subjectteacherpairs = Matrix{Int}(undef, numpairs, 2)
        i = 0
        for (subject, teachers) in subjectteacher_maps[div]
            for teacher in teachers
                i += 1
                subjectteacherpairs[i, 1:2] = [subject, teacher]
            end
        end
        # add using CS.TableSet(table)
        teachervars = divvarval_maps[:teacher][div][1]
        for (svar, tvar) in zip(subjectvars, teachervars)
            # model vars corresponding to this pair of subject and teacher vars
            (modelsvar, modeltvar) = modelvars[svar], modelvars[tvar]
            @constraint(m, [modelsvar, modeltvar] in CS.TableSet(subjectteacherpairs))
        end
    end

    return nothing
end

"""
    define_teacherconstraints!(m::Model, modelvars::OrderedDict, vardata::VariableData)::Nothing

Define teacher constraints in the model `m` using the `vardata` object and the model variables `modelvars`.

# Notes
- Modify the model `m` in place.
"""
function define_teacherconstraints!(m::Model, modelvars::OrderedDict, vardata::VariableData)::Nothing
    # variables with same row,period should have different teachers
    divvarval_maps = vardata.divvarval_maps

    # following block of code creates a list of vars
    # grouped by the same row and period
    i = 0
    while true
        i += 1
        # store var or 0 if var not found
        rowperiodvars = [get(vars, i, 0) for (_, (vars, _)) in divvarval_maps[:teacher]]
        # remove zeroes from rowperiodvars
        filter!(!iszero, rowperiodvars)
        # if no vars were pushed to list,
        # then we have exceeded the number of rows and periods
        if length(rowperiodvars) == 0
            break
        end
        # now get model vars from these vars
        rowperiodmodelvars = [modelvars[x] for x in rowperiodvars]
        # add alldifferent constraint
        @constraint(m, rowperiodmodelvars in CS.AllDifferent())
    end

    return nothing
end

"""
    define_constraints!(m::Model, modelvars::OrderedDict, vardata::VariableData)::Nothing

Define all constraints in the model `m` using the `vardata` object and the model variables `modelvars`.

# Notes
- Modify the model `m` in place.
- Call [`define_subjectconstraints!`](@ref), [`define_subjectteacherconstraints!`](@ref)
  and [`define_teacherconstraints!`](@ref) internally, in the respective order.
"""
function define_constraints!(m::Model, modelvars::OrderedDict, vardata::VariableData)::Nothing
    define_subjectconstraints!(m, modelvars, vardata)

    define_subjectteacherconstraints!(m, modelvars, vardata)

    define_teacherconstraints!(m, modelvars, vardata)

    return nothing
end

"""
    get_solution(m::Model, modelvars::OrderedDict, vardata::VariableData)::Tuple{OrderedDict{Int,Int},MOI.TerminationStatusCode}

Return a solution and status code after solving the given model `m` and using model variables `modelvars` and variable data `vardata` to access the solution values.

# Notes
- Model is optimally solved only if the status code is [`MOI.OPTIMAL`](https://jump.dev/JuMP.jl/stable/moi/reference/models/#MathOptInterface.TerminationStatusCode).
- Solution is a mapping from variable to its value (both as ints).
"""
function get_solution(
    m::Model, modelvars::OrderedDict, vardata::VariableData)::Tuple{OrderedDict{Int,Int},MOI.TerminationStatusCode}
    # solve model
    optimize!(m)

    # convert model values to var=>sol mapping
    solution = OrderedDict{Int,Int}()

    # model status
    status = JuMP.termination_status(m)

    if status != MOI.OPTIMAL
        error("Model unsolved!\nstatus = $status")
    end

    for (symbol, symbolmap) in vardata.divvarval_maps
        for (div, (vars, vals)) in symbolmap
            for var in vars
                # get the corresponding var defined in the model
                mvar = modelvars[var]
                # get its value
                mval = Int(JuMP.value(mvar))
                # add var val mapping to solution
                solution[var] = mval
            end
        end
    end

    # return solution
    return solution, status
end

"""
    convertsolution(rawsolution::OrderedDict{Int,Int}, vardata::VariableData)::OrderedDict{String,String}

Return a solution in the form of a mapping from variable to its value (both as strings).

# Notes
- Convert an `Int=>Int` solution to a `String=>String` one using the `inverse*` mappings defined in `vardata`.
"""
function convertsolution(
    rawsolution::OrderedDict{Int,Int}, vardata::VariableData)::OrderedDict{String,String}
    solution = OrderedDict{String,String}()
    # iterate over every var and val in solution
    for (var, val) in rawsolution
        # get var val strs using mappings
        var_str = vardata.varmap[var]
        val_str = vardata.valmap[val]
        # add it to solution
        solution[var_str] = val_str
    end
    return solution
end

"""
    applysolution!(schedule::Schedule, solution::OrderedDict{String,String})::Nothing

Apply a solution (of strings) to the given schedule.

# Notes
- Modify the schedule in place.
"""
function applysolution!(schedule::Schedule, solution::OrderedDict{String,String})::Nothing
    # iterate over every var and val in solution
    for (var, val) in solution
        # modify schedule with var and val
        modify!(schedule, var, val)
    end
end

"""
    solve!(schedule::Schedule, timeout::Float64 = Inf, all_solutions::Bool = false)::Tuple{OrderedDict{String,String},MOI.TerminationStatusCode}

Solve a given schedule in place and return a `String=>String` solution and status code.

# Arguments
- `schedule::Schedule`: Schedule object representing the problem.
- `timeout::Float64=Inf`: Float representing the maximum time to run the model.
- `all_solutions::Bool=false`: Bool representing whether to return all solutions or just the first.

# Notes
- Schedule is optimally solved only if the status code is [`MOI.OPTIMAL`](https://jump.dev/JuMP.jl/stable/moi/reference/models/#MathOptInterface.TerminationStatusCode).
"""
function solve!(
    schedule::Schedule, timeout::Float64 = Inf, all_solutions::Bool = false
)::Tuple{OrderedDict{String,String},MOI.TerminationStatusCode}
    # get vardata
    vardata = VariableData(schedule)

    # get model and model vars
    m, modelvars = get_model(vardata, timeout, all_solutions)

    # get solution
    rawsolution, status = get_solution(m, modelvars, vardata)

    # convert solution from ints to strs
    solution = convertsolution(rawsolution, vardata)

    # apply solution
    applysolution!(schedule, solution)

    # return solution
    return solution, status
end


text:
# # Instructions for Assignment04
#
# ## Introduction
#
# Note: this file is runnable in its current state,
# but is incomplete.
# You can run the file from the command line,
# or use the VS Code julia extension to run individual lines.

## To view additional instructions for this assignment, #jl
## see https://wellesley-bisc195.github.io/assignments/Assignment04/ #jl

## To view a rendered version of this document, #jl
## see https://wellesley-bisc195.github.io/assignments/Assignment04/#assignment04_code #jl

# ## Base Composition
#
# We've already seen how to calculate GC content of a DNA sequence
# but sometimes it might be useful to have more information
# about the overall base composition.
# To do this, we can loop through a DNA sequence
# and count each type of base as it occurs.
#
# But first, let's write a "helper" function
# that does some work to make sure that sequences
# are normalized into a particular form.
# In particular, we'll make sure they're always uppercase,
# are `String`s rather than `Char`s,
# and don't contain any invalid bases.
# Note: you can convert a `Char` to a `String` with `string(c)`.
#
# ### Question 1

"""
    normalizeDNA(sequence)

Normalizes DNA sequences or bases to uppercase `String`s.
Throws an error if invalid bases are encountered.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> normalizeDNA("aaatg")
    "AAATG"

    julia> normalizeDNA("XTG")
    ERROR: Invalid base X encountered

    julia> normalizeDNA(5)
    ERROR: Invalid base 5 encountered

    julia> normalizeDNA('G')
    "G"

    julia> b = normalizeDNA('G');

    julia> typeof(b)
    String
"""
function normalizeDNA(sequence)
    ## 1. make sequence into a String
    ## 2. make sequence uppercase

    for base in sequence
        if !occursin(base, "ACGT")
            error("Invalid base $base encountered")
        end
    end

    ## return the normalized sequence
end

# ### Question 2
#
# Now, we should be able to set `sequence = normalizeDNA(sequence)`
# in all of our other functions, and never worry about the form of the input.
#
# The following function returns a `Tuple`,
# which is an immutable structure that contains multiple values.

"""
    basecomposition(sequence)

Counts the number of each type of base
in a DNA sequence and returns a tuple those counts
in the order A, C, G, T

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> basecomposition("AATCGGG")
    (2, 1, 3, 1)

    julia> basecomposition('C')
    (0, 1, 0, 0)

    julia> A,C,G,T = basecomposition("accgggtttt")
    (1, 2, 3, 4)

    julia> A
    1

    julia> T
    4

    julia> basecomposition("BACCGGGTTTT")
    ERROR: Invalid base B encountered
"""
function basecomposition(sequence)
    sequence = normalizeDNA(sequence) # make uppercase string, check invalid bases
    a = c = g = t = 0 # sets all 4 variables to `0`

    for base in sequence
        ## add 1 to each base as it occurs
    end

    return a,c,g,t
end

# Hint: If all of your base counts are zeros,
# What's the type of `base` as you loop through a `String`?

# ### Question 3
#
# Now let's rewrite our `gc_content()` function from Assignment03
# using our `basecomposition()` function.
# This one has the same docstring,
# but should require substantially less code.

"""
    gc_content(sequence)

Calculates the GC ratio of a DNA sequence.
The GC ratio is the total number of G and C bases divided by the total length of the sequence.
For more info about GC content, see here:

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> gc_content("AATG")
    0.25

    julia> gc_content("cccggg") * 100
    100.0

    julia> gc_content("ATta")
    0.0

    julia> gc_content("ATty")
    Error: Invalid base Y encountered
"""
function gc_content(sequence)
    ## Be sure to use `basecomposition()` in your answer.
    ## Note: Since `basecomposition()` already calls `normalizeDNA`,
    ## there's no need to call it here.
end

# ## Kmer Composition
#
# Be sure that you've read the sections of Lesson 5
# on kmers and dictionaries before attempting the following questions.
#
# What we'd like to do is to get a count of all of the unique kmers in a DNA sequence
# But we can't simply generate all possible kmers
# (there are too many for even small-ish values of k).
#
# ### Question 4
#
# Let's think of the process in pseudocode,
# a plain-language high-level description of what we want to happen.
# Here, `sequence` represents some DNA sequence,
# and `k` is an integer representing the length of each kmer.
#
# ```
# make a dictionary to store kmers
# for each base in sequence
#   get a kmer of the base and the next (k-1) bases
#   make sure the kmer is a valid DNA sequence
#   if this kmer is a key in the dictionary
#       add 1 to the value referenced by that kmer
#   otherwise
#       make a new entry in the dictionary with a value of 1
# ```

"""
   kmercount(sequence, k)

Finds all kmers in a sequence,
returning a dictionary of those kmers
and the number of times they appear in the sequence.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> kmercount("ggg", 3)
    Dict{Any,Any} with 1 entry:
    "GGG" => 1

    julia> kmercount("ATATATATA", 4)
    Dict{Any,Any} with 2 entries:
    "TATA" => 3
    "ATAT" => 3

    julia> kmercount("ATATATATAx", 4)
    ERROR: Invalid base X encountered

    julia> kmercount("A", 2)
    ERROR: k must be a positive integer less than the length of the sequence
"""
function  kmercount(sequence, k)
    1 <= k <= length(sequence) || error("k must be a positive integer less than the length of the sequence")
    kmers = Dict() # initialize dictionary

    ## We're going to loop through the string with numerical index,
    ## each time grabbing the bases at position i through i+k-1.
    ## What is the last index that we should search?
    stopindex = 0

    for i in 1:stopindex
        kmer = "" # Change to index the sequence from i to i+k-1
        kmer = normalizeDNA(kmer)
        ##   if this kmer is a key the dictionary
        ##       add 1 to the value referenced by that kmer
        ##   otherwise
        ##       make a new entry in the dictionary with a value of 1

    end
    return kmers
end

ext: jl


text:
function get_parameters_drivers(obj::IonosphericModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function path_delay(obj::IonosphericModel, arg0::FieldSpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{RealFieldElement})
    return jcall(obj, "pathDelay", RealFieldElement, (FieldSpacecraftState, TopocentricFrame, jdouble, Vector{RealFieldElement}), arg0, arg1, arg2, arg3)
end

function path_delay(obj::IonosphericModel, arg0::SpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{jdouble})
    return jcall(obj, "pathDelay", jdouble, (SpacecraftState, TopocentricFrame, jdouble, Vector{jdouble}), arg0, arg1, arg2, arg3)
end


text:
using ViscousFlow
using Test
##using TestSetExtensions
using Literate
using Suppressor

const GROUP = get(ENV, "GROUP", "All")

ENV["GKSwstype"] = "nul" # removes GKS warnings during plotting

notebookdir = "../examples"
docdir = "../docs/src/manual"
litdir = "./literate"

if GROUP == "All" || GROUP == "Auxiliary"
  include("pointforce.jl")
end

if GROUP == "Literate"
  for (root, dirs, files) in walkdir(litdir)
    for file in files
      endswith(file,".jl") && @testset "$file" begin include(joinpath(root,file)) end
    end
  end
end

if GROUP == "All" || GROUP == "Notebooks"
  for (root, dirs, files) in walkdir(litdir)
    for file in files
      #endswith(file,".jl") && startswith(file,"8") && Literate.notebook(joinpath(root, file),notebookdir)
      endswith(file,".jl") && Literate.notebook(joinpath(root, file),notebookdir)
    end
  end
end

if GROUP == "Documentation"
  for (root, dirs, files) in walkdir(litdir)
    for file in files
      endswith(file,".jl") && Literate.markdown(joinpath(root, file),docdir)
    end
  end
end


text:
function InverseLUSolve(A::Array,B::Vector)
      #A = the matrix of coefficients 'A' must be squared
  #B = The Vector 'B' right hand side vector Vertical
    m,n =size(A)
    v = size(B,1)
    if m != n
        error("Matrix must be squared")
    elseif m != v
        error("Vector Length must be equal to Matrix")
    end

    N,P=InverseLU(A)

    X=N*(P*B)

    return X
end


text:
#!/usr/bin/julia
module dpd_functions
include("dpd_io.jl")

function wR(r, rc=1.0)
    nr = norm(r)
    wr = 0.0
    if nr <= rc
        wr = 1 - nr / rc
    end
    wr
end

function tot_PE(X, A, bl, L, rc=1.0)
    N = size(X, 1)
    box = L * eye(3)
    inv_box = pinv(box)
    g = zeros(3); dr = zeros(3)
    E = 0.0

    for i in 1:N
        for j in 1:i-1
            dr = X[i, :] - X[j, :]
            g = inv_box * dr
            g = g - round(g)
            dr = box * g
            E += A[bl[i], bl[j]] * rc * wR(norm(dr), rc)^2 / 2.0
        end
    end
    E
end

function tot_KE(V)
    KE = sum(V .* V) / 2.0
end

function init_pos(N, L)
    X = rand(N, 3) .* L
end

function init_vel(N, kT)
    V = randn(N, 3) * kT
    V - ones(N, 3) .* mean(V, 1)
end

function F_tot(r, v, a, gamma, kT, dt, rc)
    nr = norm(r)
    FC = a * wR(r, rc) * r / nr
    FD = -gamma * wR(r, rc)^2 * dot(r, v) * r / nr^2
    FR = sqrt(2 * gamma * kT) * wR(r, rc) * randn() / sqrt(dt) * r / nr
    ftot = FC + FD + FR
    ftot
end

function force_mat(X, V, A, bl, L, gamma, kT, dt, rc)
    """Force matrix. Output: (N, N, 3) matrix"""
    N = size(X, 1)
    fcube = zeros(N, N, 3)
    box = L * eye(3)
    inv_box = pinv(box)
    g = zeros(3); rij = zeros(3); vij = zeros(3)
    F = zeros(N, 3)

    for i in 1:N
        for j in 1:i-1
            rij = X[i, :] - X[j, :]
            g = inv_box * rij
            g = g - round(g)
            rij = box * g
            vij = V[i, :] - V[j, :]
            fcube[i, j, :] = F_tot(rij, vij, A[bl[i], bl[j]], gamma, kT, dt, rc)
            fcube[j, i, :] = -fcube[i, j, :]
        end
    end
    F = sum(fcube, 2)[:, 1, :] # sum rows
end

function integrate_euler(X, V, A, bl, L, gamma, kT, dt, rc, Nt, thermo)
    N = size(X, 1)
    T = zeros(Nt+1)
    KE = zeros(Nt+1)
    PE = zeros(Nt+1)
    F = zeros(size(X))
    KE[1] = tot_KE(V)
    PE[1] = tot_PE(X, A, bl, L, rc)
    @printf "Initial data | KE: %.6e | PE: %.6e\n" KE[1] PE[1]

    ti = time()
    for i in 2:Nt+1
        F = force_mat(X, V, A, bl, L, gamma, kT, dt, rc)
        V = V + F * dt
        X = X + V * dt

        X = mod(X, L)
        KE[i] = tot_KE(V)
        PE[i] = tot_PE(X, A, bl, L, rc)
        T[i] = KE[i] / (3.0 * N / 2.0)
        tf = time()
        if i % thermo == 0
            writedlm("Dump/dump_$(i).xyz", [bl X])   # rewrite
#            write_xyz("Dump/dump_$(i).xyz", bl, X)
            @printf "Step %i | t: %.3f | T: %.3f | PE: %.3e | Time: %.2f\n" i  i * dt  T[i]  PE[i]  tf - ti
        end
    end
    T, KE, PE
end

function integrate_verlet(X, V, A, bl, L, gamma, kT, dt, rc, Nt, thermo)
    N = size(X, 1)
    T = zeros(Nt+1)
    KE = zeros(Nt+1)
    PE = zeros(Nt+1)
    Fnew = zeros(size(X))

    writedlm("init_force.out", F)
    KE[1] = tot_KE(V)
    T[1] = KE[1] / (3.0 * N / 2)
    PE[1] = tot_PE(X, A, bl, L, rc)
    write_xyz("Dump/dump_1.xyz", bl, X)
    @printf "Initial data | KE: %.6e | PE: %.6e\n" KE[1] PE[1]

    ti = time()
    F = force_mat(X, V, A, bl, L, gamma, kT, dt, rc)
    for i in 2:Nt+1
        X = X + V * dt + 0.5 * F * dt^2
        Fnew = force_mat(X, V, A, bl, L, gamma, kT, dt, rc)
        V = V + 0.5 * (F + Fnew) * dt
        F = Fnew

        X = mod(X, L)
        KE[i] = tot_KE(V)
        PE[i] = tot_PE(X, A, bl, L, rc)
        T[i] = KE[i] / (3.0 * N / 2.0)
        tf = time()
        if i % thermo == 0
            write_xyz("Dump/dump_$(i).xyz", bl, X)
            @printf "Step: %i | t: %.3f | T: %.3f | PE: %.3e | Time: %.2f\n" i  i * dt  T[i]  PE[i]  tf - ti
        end
    end
    T, KE, PE
end

end

max_stars_repo_name: petervanya/DPDsim

max_issues_repo_name: petervanya/DPDsim

max_forks_repo_name: petervanya/DPDsim


text:
function _x_in_y(x, y::AbstractArray)::Bool
    for i = 1:length(y)
        if isassigned(y, i) && _compare(x, y[i])
            return true
        end
    end
    return false
end

_x_in_y(x::Missing, y::AbstractArray)::Bool = false

function _how_many_elements_occur_in_this_array(elements::_DataElements{T}, arr::AbstractArray; kwargs...)::Int where T
    temp::Vector{Bool} = Vector{Bool}(undef, length(elements.v))
    for i = 1:length(elements.v)
        temp[i] = _x_in_y(elements.v[i], arr)
    end
    result::Int = sum(temp)
    @debug("$(result) element(s) occured in the array.")
    return result
end

function _sanitize!(arr::AbstractArray, data::Vector{Data}, elements::_DataElements; required_matches::Integer = 5, kwargs...)
    for i = 1:length(arr)
        if isassigned(arr, i)
            _sanitize!(arr[i], data::Vector{Data}, elements::_DataElements; kwargs...)
        end
    end
    if _how_many_elements_occur_in_this_array(elements, arr; kwargs...) >= required_matches
        zero!(arr)
    end
    return arr
end

function _elements!(all_elements::Vector{Any}, arr::AbstractArray; kwargs...)
    # push!(all_elements, arr)
    for object in arr
        _elements!(all_elements, object; kwargs...)
    end
    return all_elements
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libfreenect_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libfreenect")
JLLWrappers.@generate_main_file("libfreenect", UUID("05d52523-c636-5a6c-8a6a-526edab90d71"))
end  # module libfreenect_jll


text:
struct FixedAdmittance <: ElectricLoad
    name::String
    available::Bool
    bus::Bus
    Y::Complex{Float64} # [Z]
end


text:
__precompile__()
"""
This module checks a given safety property without explicitly computing the reachable states.
For example, to check whether x_42 is less than 5, we just need to track the maximum values of x_42, i.e., evaluate the support of x_42 in the positive infinity direction.
This is cheaper than computing the reachable states also for the negative direction (and also for the other variable in the respective block).
"""
module Properties

import LazySets.Approximations:decompose,
                               overapproximate
import Reachability.tocc

using LazySets, MathematicalSystems, ..Utils, Expokit, ProgressMeter

# ==============================
# Property struct and evaluation
# ==============================
include("Property.jl")
export Property,
       inout_map_property

include("LinearConstraintProperty.jl")
export LinearConstraintProperty,
       Clause

include("IntersectionProperty.jl")
export IntersectionProperty

include("SubsetProperty.jl")
export SubsetProperty

# dictionary of registered algorithms
available_algorithms = Dict{String, Dict{String, Any}}()

# "explicit" backends
include("check_blocks.jl")
push!(available_algorithms, "explicit_blocks"=>Dict("func"=>check_blocks,
                                                    "is_full"=>false,
                                                    "is_explicit"=>true))

include("check_property.jl")
export check_property

# ====================================================
# Algorithms to find a threshold for property checking
# ====================================================
include("tune.jl")
export tune_δ

end #module Properties


text:
module DLagrangianFESpacesTests

using Test
using Gridap

model = CartesianDiscreteModel(partition=(2,2))

grid = Grid(model)

facelabels = FaceLabels(model)

node_to_label = labels_on_dim(facelabels,0)

tag_to_labels = facelabels.tag_to_labels

diritags = [1,2,5]
dirimasks = [true,true,true]

fespace = DLagrangianFESpace(
  Float64,grid,node_to_label,tag_to_labels,diritags,dirimasks)

r = [[-1, -2, 1, 2], [-3, -4, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]

@test collect(fespace.cell_to_dofs) == r

trian = Triangulation(grid)
quad = CellQuadrature(trian,degree=2)

bh = FEBasis(fespace)

uh = zero(fespace)

cellmat = integrate(inner(bh,bh),trian,quad)
cellvec = integrate(inner(bh,uh),trian,quad)

ufun(x) = x[1]

nfree = 12
ndiri = 4

test_fe_space(fespace, nfree, ndiri, cellmat, cellvec, ufun)

uh = interpolate(fespace,ufun)
u = CellField(trian,ufun)
e = u - uh

tol = 1.0e-8
@test sum(integrate(inner(e,e),trian,quad)) < tol

ufun1(x) = 1.0
ufun2(x) = 2.0
ufun5(x) = 5.0

dv = interpolate_diri_values(fespace,[ufun1,ufun2,ufun5])

@test dv == [1.0, 5.0, 5.0, 2.0]

# DLagrangianFESpace for vector values

T = VectorValue{2,Float64}

diritags = [1,2,5]
dirimasks = [(true,false),(false,true),(true,true)]

fespace = DLagrangianFESpace(
  T,grid,node_to_label,tag_to_labels,diritags,dirimasks)

bh = FEBasis(fespace)

uh = zero(fespace)

cellmat = integrate(inner(bh,bh),trian,quad)
cellvec = integrate(inner(bh,uh),trian,quad)

ufun(x) = x

nfree = 26
ndiri = 6

test_fe_space(fespace, nfree, ndiri, cellmat, cellvec, ufun)

uh = interpolate(fespace,ufun)
u = CellField(trian,ufun)
e = u - uh

tol = 1.0e-8
@test sum(integrate(inner(e,e),trian,quad)) < tol

ufun1(x) = VectorValue(1.0,-1.0)
ufun2(x) = VectorValue(2.0,-2.0)
ufun5(x) = VectorValue(5.0,-5.0)

dv = interpolate_diri_values(fespace,[ufun1,ufun2,ufun5])

@test dv == [1.0, 5.0, -5.0, 5.0, -5.0, -2.0]

order = 1
fespace = DLagrangianFESpace(T,model,order,diritags,dirimasks)
test_fe_space(fespace, nfree, ndiri, cellmat, cellvec, ufun)

order = 2
T = VectorValue{2,Float64}
diritags = [1,2,5]
dirimasks = [(true,false),(false,true),(true,true)]

fespace = DLagrangianFESpace(T,model,order,diritags,dirimasks)

@test fespace.dirimasks === dirimasks

fespace = DLagrangianFESpace(T,model,order,diritags)

@test fespace.dirimasks == [(true, true), (true, true), (true, true)]

order = 2
T = Float64
diritags = ["physical_tag_1","physical_tag_2","physical_tag_5"]
dirimasks = [false,false,true]
fespace = DLagrangianFESpace(T,model,order,diritags,dirimasks)
@test fespace.dirimasks === dirimasks
@test fespace.diritags == [1,2,5]

fespace = DLagrangianFESpace(T,model,order,diritags)
@test fespace.dirimasks == [(true,), (true,), (true,)]

fespace = DLagrangianFESpace(T,model,order)
@test fespace.diritags == Int[]
@test fespace.dirimasks == Tuple{Bool}[]
@test fespace.ndiridofs == 0

labels = FaceLabels(model)
labels = FaceLabels(
  labels.dim_to_nface_to_label,copy(labels.tag_to_labels),copy(labels.tag_to_name))
diritag = "diri_boundary"
add_tag_from_tags!(labels,diritag,diritags)
dirimask = [true,]
fespace = DLagrangianFESpace(T,model,labels,order,diritag,dirimask)

end # module


text:
module KhepriRevit
using KhepriBase
using Sockets
using EzXML

# functions that need specialization
include(khepribase_interface_file())
include("Revit.jl")

function __init__()
  add_current_backend(revit)
end
end


text:
using Test
using NBodyIPs
using NBodyIPs.Polys
using NBodyIPs: BondLengthDesc, AnalyticTransform, OneBody, BASIS
using JuLIP: save_json, load_json, decode_dict

##
println("Check (De-)Dictionisation of `BondLengthDesc`")
D = BondLengthDesc(PolyTransform(3, 2.0), CosCut(6.0, 9.0))
Ds = Dict(D)
D1 = BondLengthDesc(Ds)
println(@test D1 == D)
println(@test hash(BASIS(), D) == hash(BASIS(), D1))
##
println("generate some basis functions")
rcuts = [9.2, 6.2, 4.5]
TRANSFORM = PolyTransform(3, 2.9)
CUTOFF = [CosCut(0.66*rcut, rcut) for rcut in rcuts]
B1 = [OneBody(1.0)]
B2 = blpolys(2, TRANSFORM, CUTOFF[1], 12)
B3 = blpolys(3, TRANSFORM, CUTOFF[2], 10)
B4 = blpolys(4, TRANSFORM, CUTOFF[3], 6)
B = [B1; B2; B3; B4]
c = rand(length(B))
IP = NBodyIP(B, c)
println("Check (De-)Dictionisation of `NBodyIP`")
D_IP = Dict(IP)
IP1 = NBodyIP(D_IP)
println(@test Dict(IP1) == D_IP)
println(@test IP1 == IP)

##
# store it as a JSON file and load it again
fname = tempname() * ".json"
save_ip(fname, IP)
IP2, _ = load_ip(fname)
println(@test IP2 == IP)
run(`rm $fname`)

##
# We can also see what happens if we store a vector of basis functions
# and then try to recombine them?
# (and this also tests a bit the storing and loading of basis sets)
B_dict = Dict.(B)
fname = tempname() * ".json"
save_json(fname, Dict("B" => B_dict))
B_dict2 = load_json(fname)["B"]
run(`rm $fname`)
Ba = decode_dict.(B_dict2)
H = hash.(Ref(BASIS()), B)
Ha = hash.(Ref(BASIS()), Ba)
println(@test(H == Ha))
# check that the hashs reproduce the right number of basis functions
Hu = unique(H)
println(@test(length(Hu) == 4))
lens = [length(findall(H .== Hu[n])) for n = 1:4]
println(@test(sort(lens) == length.([B1, B2, B3, B4])))

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
import Base:  ==, hash
"""
    Location()
The Agent or Patch location lying on a Field.
"""
mutable struct Location{P<:Position}
    pos::P
    object::Union{Agent,Patch}
end

function ==(x::Location{P}, y::Location{P}) where {P<:Position}
    x.pos==y.pos && x.object==y.object
end

function hash(x::Location{P}, h::UInt) where {P<:Position}
    hash((x.pos, x.object), h)
end


text:
struct CircularArray{T,N} <: AbstractArray{T,N}
    x::AbstractArray{T,N}
    function CircularArray(x::AbstractArray{T,N}) where {T,N}
        return new{T,N}(x)
    end
end

Base.size(A::CircularArray) = size(A.x) #size of the vector

Base.length(A::CircularArray)=length(A.x)

function Base.getindex(A::CircularArray, I::Vararg{Int, N}) where N
I2 = size(A)
return Base.getindex(A.x,(mod1.(I,I2))...)
end

Base.getindex(A::CircularArray, I) = [A[i] for i in I]

function Base.setindex!(A::CircularArray,value,I...)
    I2 = size(A)
    return Base.setindex!(A.x,value,(mod1.(I,I2))...)
end

Base.iterate(S::CircularArray, state=1) = Base.iterate(S.x)
Base.IndexStyle(::Type{CircularArray}) = IndexCartesian()
Base.checkbounds(A::CircularArray, I...) = nothing


text:
"""`care(A, B, Q, R)`

Compute 'X', the solution to the continuous-time algebraic Riccati equation,
defined as A'X + XA - (XB)R^-1(B'X) + Q = 0, where R is non-singular.

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf
"""
function care(A, B, Q, R)
    G = try
        B*inv(R)*B'
    catch y
        if y isa SingularException
            error("R must be non-singular in care.")
        else
            throw(y)
        end
    end

    Z = [A  -G;
        -Q  -A']

    S = schur(Z)
    S = ordschur(S, real(S.values).<0)
    U = S.Z

    (m, n) = size(U)
    U11 = U[1:div(m, 2), 1:div(n,2)]
    U21 = U[div(m,2)+1:m, 1:div(n,2)]
    return U21/U11
end

"""`dare(A, B, Q, R)`

Compute `X`, the solution to the discrete-time algebraic Riccati equation,
defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where Q>=0 and R>0

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf
"""
function dare(A, B, Q, R)
    if (!ishermitian(Q) || minimum(eigvals(real(Q))) < 0)
        error("Q must be positive-semidefinite.");
    end
    if (!isposdef(R))
        error("R must be positive definite.");
    end

    n = size(A, 1);

    E = [
        Matrix{Float64}(I, n, n) B/R*B';
        zeros(size(A)) A'
    ];
    F = [
        A zeros(size(A));
        -Q Matrix{Float64}(I, n, n)
    ];

    QZ = schur(F, E);
    QZ = ordschur(QZ, abs.(QZ.alpha./QZ.beta) .< 1);

    return QZ.Z[(n+1):end, 1:n]/QZ.Z[1:n, 1:n];
end

"""`dlyap(A, Q)`

Compute the solution `X` to the discrete Lyapunov equation
`AXA' - X + Q = 0`.
"""
function dlyap(A::AbstractMatrix, Q)
    lhs = kron(A, conj(A))
    lhs = I - lhs
    x = lhs\reshape(Q, prod(size(Q)), 1)
    return reshape(x, size(Q))
end

"""`gram(sys, opt)`

Compute the grammian of system `sys`. If `opt` is `:c`, computes the
controllability grammian. If `opt` is `:o`, computes the observability
grammian."""
function gram(sys::AbstractStateSpace, opt::Symbol)
    if !isstable(sys)
        error("gram only valid for stable A")
    end
    func = iscontinuous(sys) ? lyap : dlyap
    if opt == :c
        # TODO probably remove type check in julia 0.7.0
        return func(sys.A, sys.B*sys.B')#::Array{numeric_type(sys),2} # lyap is type-unstable
    elseif opt == :o
        return func(Matrix(sys.A'), sys.C'*sys.C)#::Array{numeric_type(sys),2} # lyap is type-unstable
    else
        error("opt must be either :c for controllability grammian, or :o for
                observability grammian")
    end
end

"""`obsv(A, C)` or `obsv(sys)`

Compute the observability matrix for the system described by `(A, C)` or `sys`.

Note that checking for observability by computing the rank from `obsv` is
not the most numerically accurate way, a better method is checking if
`gram(sys, :o)` is positive definite."""
function obsv(A::AbstractMatrix, C::AbstractMatrix)
    T = promote_type(eltype(A), eltype(C))
    n = size(A, 1)
    ny = size(C, 1)
    if n != size(C, 2)
        throw(ArgumentError("C must have the same number of columns as A"))
    end
    res = fill(zero(T), n*ny, n)
    res[1:ny, :] = C
    for i=1:n-1
        res[(1 + i*ny):(1 + i)*ny, :] = res[((i - 1)*ny + 1):i*ny, :] * A
    end
    return res
end
obsv(sys::StateSpace) = obsv(sys.A, sys.C)

"""`ctrb(A, B)` or `ctrb(sys)`

Compute the controllability matrix for the system described by `(A, B)` or
`sys`.

Note that checking for controllability by computing the rank from
`ctrb` is not the most numerically accurate way, a better method is
checking if `gram(sys, :c)` is positive definite."""
function ctrb(A::AbstractMatrix, B::AbstractMatrix)
    T = promote_type(eltype(A), eltype(B))
    n = size(A, 1)
    nu = size(B, 2)
    if n != size(B, 1)
        throw(ArgumentError("B must have the same number of rows as A"))
    end
    res = fill(zero(T), n, n*nu)
    res[:, 1:nu] = B
    for i=1:n-1
        res[:, (1 + i*nu):(1 + i)*nu] = A * res[:, ((i - 1)*nu + 1):i*nu]
    end
    return res
end
ctrb(sys::StateSpace) = ctrb(sys.A, sys.B)

"""`P = covar(sys, W)`

Calculate the stationary covariance `P = E[y(t)y(t)']` of an lti-model `sys`, driven by gaussian
white noise 'w' of covariance `E[w(t)w(τ)]=W*δ(t-τ)` where δ is the dirac delta.

The ouput is if Inf if the system is unstable. Passing white noise directly to
the output will result in infinite covariance in the corresponding outputs
(D*W*D' .!= 0) for contunuous systems."""
function covar(sys::AbstractStateSpace, W)
    (A, B, C, D) = (sys.A, sys.B, sys.C, sys.D)
    if !isa(W, UniformScaling) && (size(B,2) != size(W, 1) || size(W, 1) != size(W, 2))
        error("W must be a square matrix the same size as `sys.B` columns")
    end
    if !isstable(sys)
        return fill(Inf,(size(C,1),size(C,1)))
    end
    func = iscontinuous(sys) ? lyap : dlyap
    Q = try
        func(A, B*W*B')
    catch
        error("No solution to the Lyapunov equation was found in covar")
    end
    P = C*Q*C'
    if iscontinuous(sys)
        #Variance and covariance infinite for direct terms
        directNoise = D*W*D'
        for i in 1:size(C,1)
            if directNoise[i,i] != 0
                P[i,:] .= Inf
                P[:,i] .= Inf
            end
        end
    else
        P += D*W*D'
    end
    return P
end

covar(sys::TransferFunction, W) = covar(ss(sys), W)

"""
    covar(C,W)
If `C` is a matrix, return CWC'
"""
covar(C::Union{AbstractMatrix,UniformScaling}, R) = C*R*C'

# Note: the H∞ norm computation is probably not as accurate as with SLICOT,
# but this seems to be still reasonably ok as a first step
"""
`..  norm(sys, p=2; tol=1e-6)`

`norm(sys)` or `norm(sys,2)` computes the H2 norm of the LTI system `sys`.

`norm(sys, Inf)` computes the L∞ norm of the LTI system `sys`.
The H∞ norm is the same as the L∞ for stable systems, and Inf for unstable systems.
If the peak gain frequency is required as well, use the function `norminf` instead.

`tol` is an optional keyword argument, used only for the computation of L∞ norms.
It represents the desired relative accuracy for the computed L∞ norm
(this is not an absolute certificate however).

sys is first converted to a state space model if needed.

The L∞ norm computation implements the 'two-step algorithm' in:
N.A. Bruinsma and M. Steinbuch, 'A fast algorithm to compute the H∞-norm
of a transfer function matrix', Systems and Control Letters 14 (1990), pp. 287-293.
For the discrete-time version, see, e.g.,: P. Bongers, O. Bosgra, M. Steinbuch, 'L∞-norm
calculation for generalized state space systems in continuous and discrete time',
American Control Conference, 1991.
"""
function LinearAlgebra.norm(sys::AbstractStateSpace, p::Real=2; tol=1e-6)
    if p == 2
        return sqrt(tr(covar(sys, I)))
    elseif p == Inf
        if sys.Ts == 0
            return normLinf_twoSteps_ct(sys,tol)[1]
        else
            return normLinf_twoSteps_dt(sys,tol)[1]
        end
    else
        error("`p` must be either `2` or `Inf`")
    end
end

function LinearAlgebra.norm(sys::TransferFunction, p::Real=2; tol=1e-6)
    return norm(ss(sys), p, tol=tol)
end

"""
`.. (peakgain, peakgainfrequency) = norminf(sys; tol=1e-6)`

Compute the L∞ norm of the LTI system `sys`, together with the frequency
`peakgainfrequency` (in rad/TimeUnit) at which the gain achieves its peak value `peakgain`.
The H∞ norm is the same as the L∞ for stable systems, and Inf for unstable systems.

`tol` is an optional keyword argument representing the desired relative accuracy for
the computed L∞ norm (this is not an absolute certificate however).

sys is first converted to a state space model if needed.

The L∞ norm computation implements the 'two-step algorithm' in:
N.A. Bruinsma and M. Steinbuch, 'A fast algorithm to compute the H∞-norm
of a transfer function matrix', Systems and Control Letters 14 (1990), pp. 287-293.
For the discrete-time version, see, e.g.,: P. Bongers, O. Bosgra, M. Steinbuch, 'L∞-norm
calculation for generalized state space systems in continuous and discrete time',
American Control Conference, 1991.
"""
function norminf(sys::AbstractStateSpace; tol=1e-6)
    if sys.Ts == 0
        return normLinf_twoSteps_ct(sys,tol)
    else
        return normLinf_twoSteps_dt(sys,tol)
    end
end

function norminf(sys::TransferFunction, ; tol=1e-6)
    return norminf(ss(sys), tol=tol)
end

function normLinf_twoSteps_ct(sys::AbstractStateSpace, tol=1e-6, maxIters=1000, approximag=1e-10)
    # `maxIters`: the maximum  number of iterations allowed in the algorithm (default 1000)
    # approximag is a tuning parameter: what does it mean for a number to be on the imaginary axis
    # Because of this tuning for example, the relative precision that we provide on the norm computation
    # is not a true guarantee, more an order of magnitude
    # outputs: pair of Float64, namely L∞ norm approximation and frequency fpeak at which it is achieved
    T = promote_type(numeric_type(sys), Float64)
    if sys.nx == 0  # static gain
        return (svdvals(sys.D)[1], T(0))
    end
    p = pole(sys)
    # Check if there is a pole on the imaginary axis
    pidx = findfirst(map(x->isapprox(x,0.0),real(p)))
    if !(pidx isa Nothing)
        return (T(Inf), imag(p[pidx]))
        # note: in case of cancellation, for s/s for example, we return Inf, whereas Matlab returns 1
    else
        # Initialization: computation of a lower bound from 3 terms
        lb = maximum(svdvals(sys.D))
        fpeak = T(Inf)
        (lb, idx) = findmax([lb, T(maximum((svdvals(evalfr(sys,0)))))]) #TODO remove T() in julia 0.7.0
        if idx == 2
            fpeak = T(0)
        end
        if isreal(p)  # only real poles
            omegap = minimum(abs.(p))
        else  # at least one pair of complex poles
            tmp = maximum(abs.(imag.(p)./(real.(p).*abs.(p))))
            omegap = abs(p[argmax(tmp)])    # TODO This is highly suspicious
        end
        (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys, omegap*1im))))]) #TODO remove T() in julia 0.7.0
        if idx == 2
            fpeak = omegap
        end

        # Iterations
        iter = 1;
        while iter <= maxIters
            res = (1+2*T(tol))*lb
            R = sys.D'*sys.D - res^2*I
            S = sys.D*sys.D' - res^2*I
            M = sys.A-sys.B*(R\sys.D')*sys.C
            H = [         M              -res*sys.B*(R\sys.B') ;
                   res*sys.C'*(S\sys.C)            -M'            ]
            omegas = eigvals(H) .+ 0im # To make type stable
            omegaps = imag.(omegas[ (abs.(real.(omegas)).<=approximag) .& (imag.(omegas).>=0) ])
            sort!(omegaps)
            if isempty(omegaps)
                return (1+T(tol))*lb, fpeak
            else  # if not empty, omegaps contains at least two values
                ms = [(x+y)/2 for x=omegaps[1:end-1], y=omegaps[2:end]]
                for mval in ms
                    (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys,mval*1im))))]) #TODO remove T() in julia 0.7.0
                    if idx == 2
                        fpeak = mval
                    end
                end
            end
            iter += 1
        end
        error("In norminf: The computation of the H-infinity norm did not converge in $maxIters iterations")
    end
end

# discrete-time version of normHinf_twoSteps_ct above
# The value fpeak returned by the function is in the range [0,pi)/sys.Ts (in rad/s)
function normLinf_twoSteps_dt(sys::AbstractStateSpace,tol=1e-6, maxIters=1000, approxcirc=1e-8)
    T = promote_type(numeric_type(sys), Float64)
    if sys.nx == 0  # static gain
        return (svdvals(sys.D)[1], T(0))
    end
    p = pole(sys)
    # Check first if there is a pole on the unit circle
    pidx = findfirst(map(x->isapprox(x,1.0),abs.(p)))
    if !(pidx isa Nothing)
        return (T(Inf), angle(p[pidx])/abs(T(sys.Ts)))
    else
        # Initialization: computation of a lower bound from 3 terms
        lb = T(maximum(svdvals(evalfr(sys,1))))  #TODO remove T() in julia 0.7.0
        fpeak = T(0)
        (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys,-1))))]) #TODO remove T() in julia 0.7.0
        if idx == 2
            fpeak = T(pi)
        end

        p = p[imag(p).>0]
        if ~isempty(p)  # not just real poles
            # find frequency of pôle closest to unit circle
            omegap = angle(p[findmin(abs.(abs.(p).-1))[2]])
        else
            omegap = T(pi)/2
        end
        (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys, exp(omegap*1im)))))]) #TODO remove T() in julia 0.7.0
        if idx == 2
            fpeak = omegap
        end

        # Iterations
        iter = 1;
        while iter <= maxIters
            res = (1+2*T(tol))*lb
            R = res^2*I - sys.D'*sys.D
            RinvDt = R\sys.D'
            L = [ sys.A+sys.B*RinvDt*sys.C  sys.B*(R\sys.B');
                  zeros(T, sys.nx,sys.nx)      I]
            M = [ I                                 zeros(T, sys.nx,sys.nx);
                  sys.C'*(I+sys.D*RinvDt)*sys.C     L[1:sys.nx,1:sys.nx]']
            # +0im to make type stable
            zs = eigvals(L,M) .+ 0im # generalized eigenvalues
            # are there eigenvalues on the unit circle?
            omegaps = angle.(zs[ (abs.(abs.(zs).-1) .<= approxcirc) .& (imag(zs).>=0)])
            sort!(omegaps)
            if isempty(omegaps)
                return (1+T(tol))*lb, fpeak/T(sys.Ts)
            else  # if not empty, omegaps contains at least two values
                ms = [(x+y)/2 for x=omegaps[1:end-1], y=omegaps[2:end]]
                for mval in ms
                    (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys,exp(mval*1im)))))]) #TODO remove T() in julia 0.7.0
                    if idx == 2
                        fpeak = mval
                    end
                end
            end
            iter += 1
        end
        error("In norminf: The computation of the H-infinity norm did not converge in $maxIters iterations")
    end
end

"""`S, P, B = balance(A[, perm=true])`

Compute a similarity transform `T` resulting in `B = T\\A*T` such that the row
and column norms of `B` are approximately equivalent. If `perm=false`, the
transformation will only scale `A` using diagonal `S`, and not permute `A` (i.e., set `P=I`)."""
function balance(A, perm::Bool=true)
    n = LinearAlgebra.checksquare(A)
    B = copy(A)
    job = perm ? 'B' : 'S'
    ilo, ihi, scaling = LAPACK.gebal!(job, B)

    S = diagm(0 => scaling)
    for j = 1:(ilo-1)   S[j,j] = 1 end
    for j = (ihi+1):n   S[j,j] = 1 end

    P = Matrix{Int}(I,n,n)
    if perm
        if ilo > 1
            for j = (ilo-1):-1:1 cswap!(j, round(Int, scaling[j]), P) end
        end
        if ihi < n
            for j = (ihi+1):n    cswap!(j, round(Int, scaling[j]), P) end
        end
    end
    return S, P, B
end

function cswap!(i::Integer, j::Integer, X::StridedMatrix)
    for k = 1:size(X,1)
        X[i, k], X[j, k] = X[j, k], X[i, k]
    end
end

"""
`sysr, G = balreal(sys::StateSpace)`

Calculates a balanced realization of the system sys, such that the observability and reachability gramians of the balanced system are equal and diagonal `G`

See also `gram`, `baltrunc`

Glad, Ljung, Reglerteori: Flervariabla och Olinjära metoder
"""
function balreal(sys::ST) where ST <: AbstractStateSpace
    P = gram(sys, :c)
    Q = gram(sys, :o)

    Q1 = try
        cholesky(Hermitian(Q)).U
    catch
        throw(ArgumentError("Balanced realization failed: Observability grammian not positive definite, system needs to be observable"))
    end
    U,Σ,V = svd(Q1*P*Q1')
    Σ .= sqrt.(Σ)
    Σ1 = diagm(0 => sqrt.(Σ))
    T = Σ1\(U'Q1)

    Pz = T*P*T'
    Qz = inv(T')*Q*inv(T)
    if norm(Pz-Qz) > sqrt(eps())
        @warn("balreal: Result may be inaccurate")
        println("Controllability gramian before transform")
        display(P)
        println("Controllability gramian after transform")
        display(Pz)
        println("Observability gramian before transform")
        display(Q)
        println("Observability gramian after transform")
        display(Qz)
        println("Singular values of PQ")
        display(Σ)
    end

    sysr = ST(T*sys.A/T, T*sys.B, sys.C/T, sys.D, sys.Ts), diagm(0 => Σ)
end

"""
`sysr, G = baltrunc(sys::StateSpace, atol = √ϵ, rtol=1e-3, unitgain=true)`

Reduces the state dimension by calculating a balanced realization of the system sys, such that the observability and reachability gramians of the balanced system are equal and diagonal `G`, and truncating it such that all states corresponding to singular values less than `atol` and less that `rtol σmax` are removed. If `unitgain=true`, the matrix `D` is chosen such that unit static gain is achieved.

See also `gram`, `balreal`

Glad, Ljung, Reglerteori: Flervariabla och Olinjära metoder
"""
function baltrunc(sys::ST; atol = sqrt(eps()), rtol = 1e-3, unitgain = true) where ST <: AbstractStateSpace
    sysbal, S = balreal(sys)
    S = diag(S)
    S = S[S .>= atol]
    S = S[S .>= S[1]*rtol]
    n = length(S)
    A = sysbal.A[1:n,1:n]
    B = sysbal.B[1:n,:]
    C = sysbal.C[:,1:n]
    D = sysbal.D
    if unitgain
        D = D/(C*inv(-A)*B)
    end

    return ST(A,B,C,D,sys.Ts), diagm(0 => S)
end

"""
    syst = similarity_transform(sys, T)
Perform a similarity transform `T : Tx̃ = x` on `sys` such that
```
Ã = T⁻¹AT
B̃ = T⁻¹ B
C̃ = CT
D̃ = D
```
"""
function similarity_transform(sys::ST, T) where ST <: AbstractStateSpace
    Tf = factorize(T)
    A = Tf\sys.A*T
    B = Tf\sys.B
    C = sys.C*T
    D = sys.D
    ST(A,B,C,D,sys.Ts)
end

"""
sysi = innovation_form(sys, R1, R2)
sysi = innovation_form(sys; sysw=I, syse=I, R1=I, R2=I)

Takes a system
```
x' = Ax + Bu + w ~ R1
y  = Cx + e ~ R2
```
and returns the system
```
x' = Ax + Kv
y  = Cx + v
```
where `v` is the innovation sequence.

If `sysw` (`syse`) is given, the covariance resulting in filtering noise with `R1` (`R2`) through `sysw` (`syse`) is used as covariance.

See Stochastic Control, Chapter 4, Åström
"""
function innovation_form(sys::ST, R1, R2) where ST <: AbstractStateSpace
    K = kalman(sys, R1, R2)
    ST(sys.A, K, sys.C, Matrix{eltype(sys.A)}(I, sys.ny, sys.ny), sys.Ts)
end
# Set D = I to get transfer function H = I + C(sI-A)\ K
function innovation_form(sys::ST; sysw=I, syse=I, R1=I, R2=I) where ST <: AbstractStateSpace
	K = kalman(sys, covar(sysw,R1), covar(syse, R2))
	ST(sys.A, K, sys.C, Matrix{eltype(sys.A)}(I, sys.ny, sys.ny), sys.Ts)
end


text:
using Test

include("namebadge.jl")

@testset "Full data" begin
    @test badge(455, "Mary M. Brown", "MARKETING") == "[455] - Mary M. Brown - MARKETING"
end

@testset "Full data with mixed case department" begin
    @test badge(7, "Eneus", "Wizardry & Witchcraft") == "[7] - Eneus - WIZARDRY & WITCHCRAFT"
end

@testset "No ID" begin
    @test badge(nothing, "Amália Valéria Kóbor", "Security") == "Amália Valéria Kóbor - SECURITY"
end

@testset "Owner's Badge" begin
    @test badge(1, "Anna Johnson", nothing) == "[1] - Anna Johnson - OWNER"
end

@testset "Owner's Badge without ID" begin
    @test badge(nothing, "Elvira McFee", nothing) == "Elvira McFee - OWNER"
end


text:
# API 581 Corrosion and regression
# Linear Interpolation

using CSV
using Statistics
using Missings
using Plots
using StatsBase

y1 = 0.005
x1 = 225
y2 = 0.01
x2 = 160
x = 200
y = y1 + ((x - x1) * ((y2-y1) / (x2 - x1)))
y = 0.0041 ipy # 240 operating temperature

operating_temperature = [10,18,43,90,160,225,275,325,350]
marine_corrosion_rate = [0,1,5,5,10,5,2,1,0]
marine_corrosion_rate = marine_corrosion_rate ./ 1000
plot(operating_temperature,marine_corrosion_rate, seriestype=:scatter,title="API 581 Table 16.2 - Corrosion Rates (ipy) - CUI", label ="Estimated Corrosion Rate API 581")

# First set of data
y = marine_corrosion_rate[1:5]
x = operating_temperature[1:5]
# Regression coefficients (y = b0 .+ (b1 .* x))
b1 = cov(x,y) / var(x) # b (slope)
b0 = mean(y) - b1 * mean(x) # a (intercept)

# Apply the regression formula to the training data
y_hat_train = fill(0.0,size(y,1)) # initialize array
y_hat_train = b0 .+ (b1 .* x)
training_set_residuals = y .- y_hat_train
r_2 = cor(x,y) ^2
n = size(train_set,1)
standard_error = sqrt((sum(y_hat_train .- y)^2) / (n-2))

plot!(x,y_hat_train, smooth=true, label = "Least Squares Regression")
plot_annotation = "y = ",b0," + (",b1," * x)"

#annotate!(30, 0.0018, text(plot_annotation, :red, :left, 10))

# Second set of data
y = marine_corrosion_rate[5:9]
x = operating_temperature[5:9]
# Regression coefficients (y = b0 .+ (b1 .* x))
b1 = cov(x,y) / var(x) # b (slope)
b0 = mean(y) - b1 * mean(x) # a (intercept)

# Apply the regression formula to the training data
y_hat_train = fill(0.0,size(y,1)) # initialize array
y_hat_train = b0 .+ (b1 .* x)
training_set_residuals = y .- y_hat_train
r_2 = cor(x,y) ^2
plot!(x,y_hat_train, smooth=true, label = "Least Squares Regression")

# margin of error
upper_bound = y + (.5 * y)
lower_bound = y - (.5 * y)
range = collect(lower_bound:.001:upper_bound)
difference_to_tmin = .280 - .120
time_to_governing_tmin = zeros(size(range,1))
for i in 1:length(range)
    time_to_governing_tmin[i] = difference_to_tmin / range[i]
end

plot(range,time_to_governing_tmin, seriestype=:scatter,title="Estimating +/-50% incorrect corrosion rate-Interpolated = 0.0041ipy", label ="Time To Governing Tmin")
plot!([y], seriestype="vline", label = "Interpolated Rate")
plot!([23.17], seriestype="hline", label = "Years MARS in-service")

#  Calculate B31.3 Pressure Tmin
# ASME B31.3 - 304.1.2 Straight Pipe Under Internal Pressure
P = 740  # Shell specification I, Design Pressure 740psig @ 100F
D = 6.625 # Outside Diameter
Y = .4 # Values of Coefficient Y for t < D/6 Table 304.1.1
E = 1.0 # Weld joint efficiency, 1.0 seamless
W = 1.0 # Weld Joint Strength Reduction Factor, W, Table 302.3.5
S = 20000 # ASTM A106 Grade B, Seamless @100F
c = 0.05 # Shell specification D, Corrosion allowance

t = (P * D) / (2*((S * E * W) + (P * Y)))
t = 0.120
tm = t + c # Note t = exclusive of corrosion allowance while tm is inclusive of corrosion allowance
tm = 0.170

max_stars_repo_name: JuliaTagBot/Machine-Learning

max_issues_repo_name: JuliaTagBot/Machine-Learning

max_forks_repo_name: JuliaTagBot/Machine-Learning


text:
@testset "get_removed_quantifier_form" begin
    @customtest String(get_removed_quantifier_form(Formula("(Q(x) | ?{x}*{y}(P(f(x),z) & Q(a))) | *{z}R(x,z,g(x)))"))) ==
        "((Q(x) | (P(f(x), z) & Q(a))) | R(x, z, g(x)))"
    @customtest String(get_removed_quantifier_form(Formula("(?{x}*{y}Foo(bar, baz)) & *{abc}~B(abc)"))) ==
        "(Foo(bar, baz) & ~B(abc))"
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Core.IR
const Compiler = Core.Compiler

# TODO: this test is broken
#let code = Any[
#        Expr(:gotoifnot, SlotNumber(2), 4),
#        Expr(:(=), SlotNumber(3), 2),
#        # Test a SlotNumber as a value of a PhiNode
#        PhiNode(Any[2,3], Any[1, SlotNumber(3)]),
#        Expr(:return, SSAValue(3))
#    ]
#
#    ci = eval(Expr(:new, CodeInfo,
#        code,
#        nothing,
#        Any[Any, Any, Any],
#        Any[Any],
#        UInt8[0, 0, 0],
#        Any[Symbol("Self"), :arg, :slot],
#        false, false, false, false
#    ))
#
#    NullLineInfo = Core.LineInfoNode(Symbol(""), Symbol(""), 0, 0)
#    Compiler.run_passes(ci, 1, [NullLineInfo])
#    # XXX: missing @test
#end

# Issue #31121
using .Compiler: CFG, BasicBlock

# We have the following CFG and corresponding DFS numbering:
#
#     CFG     DFS
#
#      A       1
#      | \     | \
#      B C     2 5
#     /|/     /|/
#    | D     | 3
#     \|      \|
#      E       4
#
# In the bug `E` got the wrong dominator (`B` instead of `A`), because the DFS
# tree had the wrong parent (i.e. we recorded the parent of `4` as `2` rather
# than `3`, so the idom search missed that `1` is `3`'s semi-dominator). Here
# we manually construct that CFG and verify that the DFS records the correct
# parent.
make_bb(preds, succs) = BasicBlock(Compiler.StmtRange(0, 0), preds, succs)
let cfg = CFG(BasicBlock[
    make_bb([]     , [2, 3]),
    make_bb([1]    , [4, 5]),
    make_bb([1]    , [4]   ),
    make_bb([2, 3] , [5]   ),
    make_bb([2, 4] , []    ),
], Int[])
    dfs = Compiler.DFS(cfg, Compiler.BBNumber(1))
    @test dfs.numbering[dfs.parents[dfs.reverse[5]]] == 4
    let correct_idoms = Compiler.naive_idoms(cfg)
        @test Compiler.SNCA(cfg) == correct_idoms
        # For completeness, reverse the order of pred/succ in the CFG and verify
        # the answer doesn't change (it does change the which node is chosen
        # as the semi-dominator, since it changes the DFS numbering).
        for (a, b, c, d) in Iterators.product(((true, false) for _ = 1:4)...)
            let cfg′ = Compiler.copy(cfg)
                a && reverse!(cfg′.blocks[1].succs)
                b && reverse!(cfg′.blocks[2].succs)
                c && reverse!(cfg′.blocks[4].preds)
                d && reverse!(cfg′.blocks[5].preds)
                @test Compiler.SNCA(cfg′) == correct_idoms
            end
        end
    end
end

# test >:
let
    f(a, b) = a >: b
    code_typed(f, Tuple{Any, Any})
    # XXX: missing @test
end

for compile in ("min", "yes")
    cmd = `$(Base.julia_cmd()) --compile=$compile interpreter_exec.jl`
    if !success(pipeline(Cmd(cmd, dir=@__DIR__); stdout=stdout, stderr=stderr))
        error("Interpreter test failed, cmd : $cmd")
    end
end

# Issue #27104
# Test whether meta nodes are still present after code optimization.
let
    @noinline f(x, y) = x + y
    @test any(code_typed(f)[1][1].code) do ex
        Meta.isexpr(ex, :meta)
    end
end

# PR #32145
# Make sure IncrementalCompact can handle blocks with predecessors of index 0
# while removing blocks with no predecessors.
let cfg = CFG(BasicBlock[
    make_bb([]        , [2, 4]),
    make_bb([1]       , [4, 5]),
    make_bb([]        , [4]   ), # should be removed
    make_bb([0, 1, 2] , [5]   ), # 0 predecessor should be preserved
    make_bb([2, 3]    , []    ),
], Int[])
    code = Compiler.IRCode(
        [], [], Int32[], UInt8[], cfg, LineInfoNode[], [], [], [])
    compact = Compiler.IncrementalCompact(code, true)
    @test length(compact.result_bbs) == 4 && 0 in compact.result_bbs[3].preds
end

# Issue #32579 - Optimizer bug involving type constraints
function f32579(x::Int, b::Bool)
    if b
        x = nothing
    end
    if isa(x, Int)
        y = x
    else
        y = x
    end
    if isa(y, Nothing)
        z = y
    else
        z = y
    end
    return z === nothing
end
@test f32579(0, true) === true
@test f32579(0, false) === false

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
using PtFEM
using Test: @test

data = Dict(
  # Rod(nels, np_types, nip, fin_el(nod, nodof))
  :struc_el => Rod(4, 1, 1, Line(2, 1)),
  :properties => [1.0e5;],
  :x_coords => range(0, stop=1, length=5),
  :support => [
    (5, [0])
    ],
  :loaded_nodes => [
    (1, [-0.625]),
    (2, [-1.25]),
    (3, [-1.25]),
    (4, [-1.25]),
    (5, [-0.625])
    ]
)

@time m, dis_fm, fm_df = p41(data)

println("\nTesting: round.(m.displacements, digits=8) == [-2.5e-5 -2.344e-5 -1.875e-5 -1.094e-5 0.0]'\n")

@test round.(m.displacements, digits=8) == [-2.5e-5 -2.344e-5 -1.875e-5 -1.094e-5 0.0]'


text:
# Compare with tabulated values,  classical hypergeometric formulas

@testset "Bessel" begin

    T = Float64

    for α in  (3/2, 1/2, 1, 2)
        P = Bessel{α,T}
        β =  2

        p0,p1,p2,p3,p4 = basis.(P, 0:4)

        # Krall and Frink
        for x in range(0, 1, length=5)
            @test p0(x) ≈ 1
            @test p1(x) ≈ 1 + α*(x/β)
            @test p2(x) ≈ 1 + 2(α+1)*(x/β) + (α+1)*(α+2)*(x/β)^2
            @test p3(x) ≈ 1 + 3(α+2)*(x/β) + 3(α+2)*(α+3)*(x/β)^2 +  (α+2)*(α+3)*(α+4)*(x/β)^3
            @test p4(x) ≈ 1 + 4(α+3)*(x/β) + 6(α+3)*(α+4)*(x/β)^2 +  4(α+3)*(α+4)*(α+5)*(x/β)^3 + (α+3)*(α+4)*(α+5)*(α+6)*(x/β)^4
        end
    end

    x = variable(Polynomial)
    for  α  ∈ (1/4,1/2,3/4, 1, 2)
        P =  Bessel{α}
        for n in 0:5
            p  = basis(P,n)
            @test  p(x) ≈ SP.classical_hypergeometric(P, n,  x)
        end
    end

end

@testset "ChebyshevT" begin

    p0 = Chebyshev([1])  # 1
    p1 = Chebyshev([0,1]) # x
    p2 = Chebyshev([0,0,1]) #
    p3 = Chebyshev([0,0,0,1]) #
    p4 = Chebyshev([0,0,0,0,1]) #

    for x in range(-1, 1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ x
        @test p2(x) ≈ 2x^2 - 1
        @test p3(x) ≈ 4x^3  - 3x
        @test p4(x) ≈ 8x^4 - 8x^2 + 1
    end

end

@testset "ChebyshevU" begin

    p0 = ChebyshevU([1])  # 1
    p1 = ChebyshevU([0,1]) # x
    p2 = ChebyshevU([0,0,1]) #
    p3 = ChebyshevU([0,0,0,1]) #
    p4 = ChebyshevU([0,0,0,0,1]) #

    for x in range(-1, 1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ 2x
        @test p2(x) ≈ 4x^2 - 1
        @test p3(x) ≈ 8x^3 - 4x
        @test p4(x) ≈ 16x^4 - 12x^2 +  1
    end

end

@testset "Hermite" begin

    p0 = Hermite([1])  # 1
    p1 = Hermite([0,1]) # x
    p2 = Hermite([0,0,1]) #
    p3 = Hermite([0,0,0,1]) #
    p4 = Hermite([0,0,0,0,1]) #

    for x in range(-1, 1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ 2x
        @test p2(x) ≈ 4x^2 - 2
        @test p3(x) ≈ 8x^3 - 12x
        @test p4(x) ≈ 16x^4 - 48x^2 + 12
    end

    x = variable(LaurentPolynomial)
    P =  Hermite
    for n in 2:5
        p  = basis(P,n)
        @test  p(x) ≈ SP.classical_hypergeometric(P, n,  x)
    end

end

@testset "ChebyshevHermite" begin

    p0 = ChebyshevHermite([1])  # 1
    p1 = ChebyshevHermite([0,1]) # x
    p2 = ChebyshevHermite([0,0,1]) #
    p3 = ChebyshevHermite([0,0,0,1]) #
    p4 = ChebyshevHermite([0,0,0,0,1]) #

    for x in range(-1, 1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ x
        @test p2(x) ≈ x^2 - 1
        @test p3(x) ≈ x^3 - 3x
        @test p4(x) ≈ x^4 - 6x^2 + 3
    end

    x = variable(LaurentPolynomial)
    P =  ChebyshevHermite
    for n in 2:5
        p  = basis(P,n)
        @test  p(x) ≈ SP.classical_hypergeometric(P, n,  x)
    end

end

@testset "Legendre" begin

    p0 = Legendre([1])  # 1
    p1 = Legendre([0,1]) # x
    p2 = Legendre([0,0,1]) #
    p3 = Legendre([0,0,0,1]) #
    p4 = Legendre([0,0,0,0,1]) #

    for x in range(-1, stop=1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ x
        @test p2(x) ≈ 1/2*(3x^2-1)
        @test p3(x) ≈ 1/2*(5x^3-3x)
        @test p4(x) ≈ 1/8*(35x^4 - 30x^2 + 3)
    end

end

@testset "Shifted Legendre" begin

    T = Float64
    P = ShiftedLegendre{T}

    p0 = P([1])  # 1
    p1 = P([0,1]) # x
    p2 = P([0,0,1]) #
    p3 = P([0,0,0,1]) #
    p4 = P([0,0,0,0,1]) #

    for x in range(0, stop=1, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ 2x - 1
        @test p2(x) ≈ 6x^2 - 6x + 1
        @test p3(x) ≈ 20x^3 - 30x^2+ 12x - 1
        @test p4(x) ≈ 70x^4 - 140x^3 + 90x^2 - 20x + 1
    end

end

@testset "Laguerre" begin

    P = Laguerre{0}
    p0 = P([1])  # 1
    p1 = P([0,1]) # x
    p2 = P([0,0,1]) #
    p3 = P([0,0,0,1]) #
    p4 = P([0,0,0,0,1]) #

    for x in range(0, stop=5, length=5)
        @test p0(x) ≈ 1
        @test p1(x) ≈ -x + 1
        @test p2(x) ≈ 1/2*(x^2 - 4x + 2)
        @test p3(x) ≈ 1/6*(-x^3 + 9x^2 -18x + 6)
        @test p4(x) ≈ 1/24 * (x^4 - 16x^3 + 72x^2 - 96x + 24)
    end

    x = variable(Polynomial)
    for  α  ∈ (0, 1/4,1/2,3/4)
        P =  Laguerre{α}
        for n in 2:5
            p  = basis(P,n)
            @test  p(x) ≈ SP.classical_hypergeometric(P, n,  x)
        end
    end

end

@testset "Gegenbauer" begin

    for α in (1/4, 1/3, 1/2, 1, 2)
        P = Gegenbauer{α, Float64}
        p0 = P([1])  # 1
        p1 = P([0,1]) # x
        p2 = P([0,0,1]) #
        p3 = P([0,0,0,1]) #
        p4 = P([0,0,0,0,1]) #

        for x in range(-1, stop=1, length=5)
            @test p0(x) ≈ 1
            @test p1(x) ≈ 2α*x
            @test p2(x) ≈ -α + 2α*(1+α)*x^2
            @test p3(x) ≈ -2α * (1 + α) *x + 4/3*α*(1+α) * (2 + α)*x^3
        end
    end

    x = variable(Polynomial)
    for  α  ∈ (1/4,1/2,3/4, 1)
        P =  Gegenbauer{α}
        for n in 2:5
            p  = basis(P,n)
            @test  p(x) ≈ SP.classical_hypergeometric(P, n,  x)
        end
    end

end

@testset "Jacobi" begin

    T = Float64
    x = variable()

    # ChebyshevT = J(-1/2, -1/2)
    for (alpha_beta, Q) in  (((-1/2, -1/2), Chebyshev{T}),
                             ((1/2,1/2), ChebyshevU{T}),
                             ((0,0), Legendre{T})
                             )

        P = Jacobi{alpha_beta..., T}
        for i in  2:6
            p = basis(P,i)
            q = basis(Q,i)
            @test SP.monic(p)(x) ≈ SP.monic(q)(x)
        end

    end

    x = variable(Polynomial)
    for  (α,β)  ∈ ((1/2,1/2),  (-1/2,1/2), (1/2,-1/2), (-1/2, -1/2))
        P =  Jacobi{α,β}
        for n in 2:5
            p  = basis(P,n)
            @test  p(x) ≈ SP.classical_hypergeometric(P, n,  x)
        end
    end

    ## We have `jacobi_eval` for α and β outside [-1, ∞)
    α,β = 1/2, -1/2
    for n in 2:5
        x0 = 1/2
        @test SP.jacobi_eval(α, β, n, x0) ≈ basis(Jacobi{α, β},n)(x0)
    end

end

## CDOP

@testset "Charlier" begin

    x = variable()
    for μ ∈ (1/4, 1/2, 1,2)
        P = Charlier{μ}
        @test basis(P,0)(x) ≈ one(x)
        @test basis(P,1)(x) ≈ 1 - x/μ
        @test basis(P,2)(x) ≈ (x^2  +  μ^2 - x *(1 + 2μ))/μ^2

    end

    for μ ∈ (1/4, 1/2, 1, 2)
        P = Charlier{μ}
        for i in 0:5
            @test basis(P,i)(x) ≈  SP.classical_hypergeometric(P,i,x)
        end
    end

end

@testset  "Hahn"  begin

    x  = variable()
    for  α ∈  (1/4, 1/2)
        for β  ∈ (1/4, 1/2)
            𝐍  = 5

            P = Hahn{α,β,𝐍}
            for i in 0:𝐍
                @test basis(P,i)(x) ≈  SP.classical_hypergeometric(P,i,x)
            end
        end
    end

end

@testset "Krawchouk" begin

    x = variable()
    for p ∈ (1/4, 1/2, 1, 2)
        for N ∈ (3,6)
            P = Krawchouk{p,N}
            @test basis(P,0)(x) ≈ one(x)
            @test basis(P,1)(x) ≈ -N*p + x
            @test basis(P,2)(x) ≈ 1/2*(N^2*p^2 + x*(2p + x - 1) - N*p*(p+2x))
        end

    end

    for p ∈ (1/4, 1/2, 1, 2)
        for N ∈ (3,6)
            P = Krawchouk{p,N}
            for i in 0:N
                @test basis(P,i)(x) ≈  SP.classical_hypergeometric(P,i,x)
            end
        end
    end

end

@testset  "Meixner" begin

    x = variable()
    for γ ∈ (1/4, 1/2, 3/4)
        for μ ∈ (1/4, 1/2, 3/4)
            P = Meixner{γ,μ}
            for i in 0:4
                @test basis(P,i)(x) ≈  SP.classical_hypergeometric(P,i,x)
            end
        end
    end

end


text:
struct InterfaceQuadratures
    quads::Any
    normals::Any
    celltoquad::Any
    ncells::Any
    function InterfaceQuadratures(quads, normals, celltoquad)
        ncells = length(celltoquad)
        nphase, nquads = size(quads)

        @assert nphase == 2
        @assert length(normals) == nquads
        @assert all(celltoquad .>= 0)
        @assert all(celltoquad .<= nquads)

        new(quads, normals, celltoquad, ncells)
    end
end

function InterfaceQuadratures(
    cellsign,
    levelset,
    levelsetcoeffs,
    nodalconnectivity,
    numqp,
    cellmap,
)

    numcells = length(cellsign)
    @assert size(nodalconnectivity)[2] == numcells

    xL, xR = [-1.0, -1.0], [1.0, 1.0]
    invjac = inverse_jacobian(cellmap)

    hasinterface = cellsign .== 0
    numinterfaces = count(hasinterface)
    quads = Matrix(undef, 2, numinterfaces)
    normals = Vector(undef, numinterfaces)
    celltoquad = zeros(Int, numcells)

    counter = 1
    for cellid = 1:numcells
        if cellsign[cellid] == 0
            nodeids = nodalconnectivity[:, cellid]
            update!(levelset, levelsetcoeffs[nodeids])

            try
                squad = surface_quadrature(levelset, xL, xR, numqp, numsplits = 2)
                n = levelset_normal(levelset, squad.points, invjac)
                quads[1, counter] = squad
                quads[2, counter] = squad
                normals[counter] = n
            catch e
                squad = surface_quadrature(levelset,xL,xR,numqp,numsplits = 3)
                n = levelset_normal(levelset, squad.points, invjac)
                quads[1, counter] = squad
                quads[2, counter] = squad
                normals[counter] = n
            end

            celltoquad[cellid] = counter
            counter += 1
        end
    end
    return InterfaceQuadratures(quads, normals, celltoquad)
end

function InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh::CutMesh, numqp)
    cellsign = cell_sign(cutmesh)
    nodalconnectivity = nodal_connectivity(cutmesh.mesh)
    cellmap = cell_map(cutmesh, 1)
    return InterfaceQuadratures(
        cellsign,
        levelset,
        levelsetcoeffs,
        nodalconnectivity,
        numqp,
        cellmap,
    )
end

function Base.getindex(iquads::InterfaceQuadratures, s, cellid)
    row = cell_sign_to_row(s)
    idx = iquads.celltoquad[cellid]
    idx > 0 || error("Cell $cellid does not have an interface quadrature rule")
    return iquads.quads[row, idx]
end

function Base.show(io::IO, interfacequads::InterfaceQuadratures)
    ncells = interfacequads.ncells
    numinterfaces = length(interfacequads.normals)
    str = "InterfaceQuadratures\n\tNum. Cells: $ncells\n\tNum. Interfaces: $numinterfaces"
    print(io, str)
end

function interface_normals(iquads::InterfaceQuadratures, cellid)
    idx = iquads.celltoquad[cellid]
    return iquads.normals[idx]
end


text:
    using LowRankIntegrators, LaTeXStrings, GLMakie
    # Reaction Network
    # ∅ -> A -> ∅, [k1/(1+(B/θ)^3), k2*A]
    # ∅ -> B -> ∅, [k3/(1+A), k4*B]

    # computing the propensities
    function props(k1, k2, k3, k4, θ, N)
        ax = [zeros(N) for i in 1:4]
        ay = [zeros(N) for i in 1:4]
        for x in 1:N
            ax[1][x] = 1
            ax[2][x] = k2*(x-1)
            ax[3][x] = k3/(1+(x-1))
            ax[4][x] = 1
        end
        for y in 1:N
            ay[1][y] = k1/(1+((y-1)/θ)^3)
            ay[2][y] = 1
            ay[3][y] = 1
            ay[4][y] = k4*(y-1)
        end
        return ax, ay
    end

    # shift operators:
    # let A be an n×n matrix. Then,
    # shift(x,n)*A shifts rows by x rows down and padds the first x rows with 0s
    # A*shift(x,n)' shifts cols by x cols right and padds the x leftmost cols with 0s
    function shift(x,n)
        S = zeros(n,n)
        if x > 0
            S[1+x:end,1:end-x] = Matrix(I,n-x, n-x)
        end
        if x < 0
            S[1:end+x,1-x:end] = Matrix(I,n+x,n+x)
        end
        if x == 0
            S = I
        end
        return S
    end

    ## simulation parameters
    const N = 600
    const Nrxn = 4
    const ν = [(1,0), (-1,0), (0,1), (0,-1)]
    k1 = 30.0
    k2 = 1
    k3 = 10
    k4 = 1
    θ = 1

    ax, ay = props(k1,k2,k3,k4,θ,N) # propensity factors
    Srows = [shift(ν[r][1],N) for r in 1:Nrxn] # shift operators
    Scols = [shift(ν[r][2],N)' for r in 1:Nrxn] # shift operators
    A = [Srows[r]*TwoFactorRepresentation(ax[r],ay[r])*Scols[r] for r in 1:Nrxn] # coefficient matrices
    Asum = TwoFactorRepresentation(hcat(ax...), hcat(ay...)) # coefficient matrix

    # chemical master equation
    function cme(P, (A,Asum,Srows,Scols), t)
        return sum(A[r] .* (Srows[r]*P*Scols[r]) for r in 1:Nrxn) - Asum .* P
    end

    # Gaussian initial condition
    x0 = 20
    y0 = 20
    Q = qr(randn(2,2)).Q
    D = Q*Diagonal([rand()/20, rand()/20])*Q'
    P0(x,y) = exp(-[x-x0; y-y0]'*D*[x-x0; y-y0])
    P_init = [P0(x,y) for x in 1:N, y in 1:N]
    P_init ./= sum(P_init)
    Tf = 10.0

    dt = 2e-2
    prob = ODEProblem(cme, P_init, (0, Tf), (Matrix.(A), Matrix(Asum), Srows, Scols))
    t_direct = @elapsed sol = DifferentialEquations.solve(prob, Tsit5(), saveat = dt)

    P_lr = truncated_svd(P_init, tol = 1e-6)
    @show rank(P_lr)
    lr_prob = MatrixDEProblem((x,t) -> cme(x, (A, Asum, Srows, Scols), t), P_lr, (0.0, 10.0))
    t_lr = @elapsed lr_sol = LowRankIntegrators.solve(lr_prob, UnconventionalAlgorithm(), dt)

    @show t_direct/t_lr

    # comparison movie
    fig = Figure()
    ax_lra = Axis(fig[1,1], xlabel = "A", ylabel = "B", title = "full")
    ax_full = Axis(fig[1,2], xlabel = "A", ylabel = "B", title = "low rank")
    ax_error = Axis(fig[2,1:2], yscale = log10, xlabel="time", ylabel = "error")
    limits!(ax_error, 0, Tf, 1e-8, 1e-0)
    display(fig)

    N_relevant = 40
    sol_lra = Observable(Matrix(lr_sol.Y[1][1:N_relevant, 1:N_relevant]))
    sol_full = Observable(sol.u[1][1:N_relevant, 1:N_relevant])
    error = Observable([norm(sol.u[1] - Matrix(lr_sol.Y[1]))])
    ts = Observable([0.0])
    Makie.heatmap!(ax_lra, sol_lra)
    Makie.heatmap!(ax_full, sol_full)
    Makie.lines!(ax_error, ts, error, linewidth = 2, color = :black)

    record(fig, "markov_chain.mp4", enumerate(lr_sol.t); framerate = 24) do (i,t)
        sol_lra[] = Matrix(lr_sol.Y[i][1:N_relevant,1:N_relevant])
        sol_full[] = sol(t)[1:N_relevant,1:N_relevant]
        push!(error.val, maximum(abs.(Matrix(lr_sol.Y[i][1:N_relevant, 1:N_relevant]) - sol(t)[1:N_relevant, 1:N_relevant])))
        push!(ts.val, t)
        ts[] = ts[]
        error[] = error[]
    end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
"""
$(TYPEDEF)

Epigraph of vector relative entropy cone of dimension `dim`.

    $(FUNCTIONNAME){T}(dim::Int, use_dual::Bool = false)
"""
mutable struct EpiRelEntropy{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    w_dim::Int
    v_idxs::UnitRange{Int64}
    w_idxs::UnitRange{Int64}

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    inv_hess_aux_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, SparseMatrixCSC{T, Int}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    lwv::Vector{T}
    tau::Vector{T}
    sigma::Vector{T}
    z::T
    Hiuu::T
    Hiuv::Vector{T}
    Hiuw::Vector{T}
    Hivv::Vector{T}
    Hivw::Vector{T}
    Hiww::Vector{T}
    temp1::Vector{T}
    temp2::Vector{T}

    function EpiRelEntropy{T}(
        dim::Int;
        use_dual::Bool = false,
        ) where {T <: Real}
        @assert dim >= 3
        @assert isodd(dim)
        cone = new{T}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        w_dim = cone.w_dim = div(dim - 1, 2)
        cone.v_idxs = 1 .+ (1:w_dim)
        cone.w_idxs = w_dim .+ cone.v_idxs
        return cone
    end
end

reset_data(cone::EpiRelEntropy) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.inv_hess_aux_updated =
    cone.hess_fact_updated = false)

function setup_extra_data!(cone::EpiRelEntropy{T}) where {T <: Real}
    w_dim = cone.w_dim
    cone.lwv = zeros(T, w_dim)
    cone.tau = zeros(T, w_dim)
    cone.sigma = zeros(T, w_dim)
    cone.Hiuv = zeros(T, w_dim)
    cone.Hiuw = zeros(T, w_dim)
    cone.Hivv = zeros(T, w_dim)
    cone.Hivw = zeros(T, w_dim)
    cone.Hiww = zeros(T, w_dim)
    cone.temp1 = zeros(T, w_dim)
    cone.temp2 = zeros(T, w_dim)
    return cone
end

get_nu(cone::EpiRelEntropy) = cone.dim

function set_initial_point!(arr::AbstractVector, cone::EpiRelEntropy)
    (arr[1], v, w) = get_central_ray_epirelentropy(cone.w_dim)
    @views arr[cone.v_idxs] .= v
    @views arr[cone.w_idxs] .= w
    return arr
end

function update_feas(cone::EpiRelEntropy{T}) where T
    @assert !cone.feas_updated
    u = cone.point[1]
    @views v = cone.point[cone.v_idxs]
    @views w = cone.point[cone.w_idxs]

    if all(>(eps(T)), v) && all(>(eps(T)), w)
        @. cone.lwv = log(w / v)
        cone.z = u - dot(w, cone.lwv)
        cone.is_feas = (cone.z > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::EpiRelEntropy{T}) where T
    u = cone.dual_point[1]
    @views v = cone.dual_point[cone.v_idxs]
    @views w = cone.dual_point[cone.w_idxs]

    if all(>(eps(T)), v) && u > eps(T)
        return all(u * (1 + log(vi / u)) + wi > eps(T) for (vi, wi) in zip(v, w))
    end

    return false
end

function update_grad(cone::EpiRelEntropy)
    @assert cone.is_feas
    @views v = cone.point[cone.v_idxs]
    @views w = cone.point[cone.w_idxs]
    z = cone.z
    sigma = cone.sigma
    tau = cone.tau
    g = cone.grad

    @. sigma = w / v / z
    @. tau = (cone.lwv + 1) / -z
    g[1] = -inv(z)
    @. @views g[cone.v_idxs] = -sigma - inv(v)
    @. @views g[cone.w_idxs] = -tau - inv(w)

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::EpiRelEntropy{T}) where T
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data
    u = cone.point[1]
    v_idxs = cone.v_idxs
    w_idxs = cone.w_idxs
    point = cone.point
    @views v = point[v_idxs]
    @views w = point[w_idxs]
    z = cone.z
    sigma = cone.sigma
    tau = cone.tau
    g = cone.grad

    # H_u_u, H_u_v, H_u_w parts
    H[1, 1] = abs2(g[1])
    @. @views H[1, v_idxs] = sigma / z
    @. @views H[1, w_idxs] = tau / z

    # H_v_v, H_v_w, H_w_w parts
    zinv = inv(z)
    @inbounds for (i, v_idx, w_idx) in zip(1:cone.w_dim, v_idxs, w_idxs)
        vi = point[v_idx]
        wi = point[w_idx]
        taui = tau[i]
        sigmai = sigma[i]

        H[v_idx, v_idx] = abs2(sigmai) - g[v_idx] / vi
        H[w_idx, w_idx] = abs2(taui) + (zinv + inv(wi)) / wi

        @. H[v_idx, w_idxs] = sigmai * tau
        @. H[w_idx, v_idxs] = sigma * taui
        H[v_idx, w_idx] -= zinv / vi

        @inbounds for j in (i + 1):cone.w_dim
            H[v_idx, v_idxs[j]] = sigmai * sigma[j]
            H[w_idx, w_idxs[j]] = taui * tau[j]
        end
    end

    cone.hess_updated = true
    return cone.hess
end

# auxiliary calculations for inverse Hessian
function update_inv_hess_aux(cone::EpiRelEntropy{T}) where T
    @assert !cone.inv_hess_aux_updated
    point = cone.point
    @views v = point[cone.v_idxs]
    @views w = point[cone.w_idxs]
    z = cone.z

    HiuHu = zero(T)
    @inbounds for i in 1:cone.w_dim
        wi = w[i]
        vi = v[i]
        lwvi = cone.lwv[i]
        zwi = z + wi
        z2wi = zwi + wi
        wz2wi = wi / z2wi
        vz2wi = vi / z2wi
        uvvi = wi * (wi * lwvi - z)
        uwwi = wi * (z + lwvi * zwi) * wz2wi
        HiuHu += wz2wi * uvvi - uwwi * (lwvi + 1)
        cone.Hiuv[i] = vz2wi * uvvi
        cone.Hiuw[i] = uwwi
        cone.Hivw[i] = wi * vi * wz2wi
        cone.Hiww[i] = wi * zwi * wz2wi
        cone.Hivv[i] = vi * zwi * vz2wi
    end
    cone.Hiuu = abs2(z) - HiuHu
    if cone.Hiuu < zero(T)
        @warn("bad Hiuu $(cone.Hiuu)")
    end

    cone.inv_hess_aux_updated = true
    return
end

function alloc_inv_hess!(cone::EpiRelEntropy{T}) where T
    # initialize sparse idxs for upper triangle of inverse Hessian
    dim = cone.dim
    w_dim = cone.w_dim
    nnz_tri = 2 * dim - 1 + w_dim
    I = Vector{Int}(undef, nnz_tri)
    J = Vector{Int}(undef, nnz_tri)
    idxs1 = 1:dim
    @views I[idxs1] .= 1
    @views J[idxs1] .= idxs1
    idxs2 = (dim + 1):(2 * dim - 1)
    @views I[idxs2] .= 2:dim
    @views J[idxs2] .= 2:dim
    idxs3 = (2 * dim):nnz_tri
    @views I[idxs3] .= 1 .+ (1:w_dim)
    @views J[idxs3] .= (1 + w_dim) .+ (1:w_dim)
    V = ones(T, nnz_tri)
    cone.inv_hess = Symmetric(sparse(I, J, V, dim, dim), :U)
    return
end

# updates for nonzero values in the inverse Hessian
function update_inv_hess(cone::EpiRelEntropy{T}) where T
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    w_dim = cone.w_dim

    # modify nonzeros of upper triangle of inverse Hessian
    nzval = cone.inv_hess.data.nzval
    nzval[1] = cone.Hiuu
    nz_idx = 2
    @inbounds for i in 1:cone.w_dim
        nzval[nz_idx] = cone.Hiuv[i]
        nzval[nz_idx + 1] = cone.Hivv[i]
        nz_idx += 2
    end
    @inbounds for i in 1:cone.w_dim
        nzval[nz_idx] = cone.Hiuw[i]
        nzval[nz_idx + 1] = cone.Hivw[i]
        nzval[nz_idx + 2] = cone.Hiww[i]
        nz_idx += 3
    end

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiRelEntropy,
    )
    @assert cone.grad_updated
    v_idxs = cone.v_idxs
    w_idxs = cone.w_idxs
    @views v = cone.point[v_idxs]
    @views w = cone.point[w_idxs]
    z = cone.z
    sigma = cone.sigma
    tau = cone.tau

    @inbounds @views begin
        u_arr = arr[1, :]
        u_prod = prod[1, :]
        v_arr = arr[v_idxs, :]
        v_prod = prod[v_idxs, :]
        w_arr = arr[w_idxs, :]
        w_prod = prod[w_idxs, :]
        mul!(u_prod, v_arr', sigma)
        mul!(u_prod, w_arr', tau, true, true)
        @. u_prod += u_arr / z
        @. v_prod = sigma * u_prod' + (sigma * v_arr + v_arr / v - w_arr / z) / v
        @. w_prod = tau * u_prod' + (w_arr / z + w_arr / w) / w - v_arr / v / z
        @. u_prod /= z
    end

    return prod
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiRelEntropy,
    )
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)

    @inbounds @views begin
        u_arr = arr[1, :]
        u_prod = prod[1, :]
        v_arr = arr[cone.v_idxs, :]
        v_prod = prod[cone.v_idxs, :]
        w_arr = arr[cone.w_idxs, :]
        w_prod = prod[cone.w_idxs, :]
        @. u_prod = cone.Hiuu * u_arr
        mul!(u_prod, v_arr', cone.Hiuv, true, true)
        mul!(u_prod, w_arr', cone.Hiuw, true, true)
        mul!(v_prod, cone.Hiuv, u_arr')
        mul!(w_prod, cone.Hiuw, u_arr')
        @. v_prod += cone.Hivv * v_arr + cone.Hivw * w_arr
        @. w_prod += cone.Hivw * v_arr + cone.Hiww * w_arr
    end

    return prod
end

function dder3(cone::EpiRelEntropy, dir::AbstractVector)
    @assert cone.grad_updated
    tau = cone.tau
    z = cone.z
    v_idxs = cone.v_idxs
    w_idxs = cone.w_idxs
    @views v = cone.point[v_idxs]
    @views w = cone.point[w_idxs]
    u_dir = dir[1]
    @views v_dir = dir[v_idxs]
    @views w_dir = dir[w_idxs]
    dder3 = cone.dder3
    @views v_dder3 = dder3[v_idxs]
    @views w_dder3 = dder3[w_idxs]
    wdw = cone.temp1
    vdv = cone.temp2

    i2z = inv(2 * z)
    @. wdw = w_dir / w
    @. vdv = v_dir / v
    const0 = (u_dir + dot(w, vdv)) / z + dot(tau, w_dir)
    const1 = abs2(const0) + sum(w[i] * abs2(vdv[i]) + w_dir[i] *
        (wdw[i] - 2 * vdv[i]) for i in eachindex(w)) / (2 * z)
    dder3[1] = const1 / z

    # v
    v_dder3 .= const1
    @. v_dder3 += (const0 + vdv) * vdv - i2z * wdw * w_dir
    @. v_dder3 *= w
    @. v_dder3 += (z * vdv - w_dir) * vdv + (-const0 + i2z * w_dir) * w_dir
    @. v_dder3 /= v
    @. v_dder3 /= z

    # w
    @. w_dder3 = const1 * tau
    @. w_dder3 += ((const0 - w * vdv / z) / z + (inv(w) + i2z) * wdw) * wdw
    @. w_dder3 += (-const0 + w_dir / z - vdv / 2) / z * vdv

    return dder3
end

# TODO remove this in favor of new hess_nz_count etc functions that directly use uu, uw, ww etc
inv_hess_nz_count(cone::EpiRelEntropy) =
    3 * cone.dim - 2 + 2 * cone.w_dim

inv_hess_nz_count_tril(cone::EpiRelEntropy) =
    2 * cone.dim - 1 + cone.w_dim

inv_hess_nz_idxs_col(cone::EpiRelEntropy, j::Int) =
    (j == 1 ? (1:cone.dim) : (j <= (1 + cone.w_dim) ?
    [1, j, j + cone.w_dim] : [1, j - cone.w_dim, j]))

inv_hess_nz_idxs_col_tril(cone::EpiRelEntropy, j::Int) =
    (j == 1 ? (1:cone.dim) : (j <= (1 + cone.w_dim) ? [j, j + cone.w_dim] : [j]))

# see analysis in
# https://github.com/lkapelevich/HypatiaSupplements.jl/tree/master/centralpoints
function get_central_ray_epirelentropy(w_dim::Int)
    if w_dim <= 10
        return central_rays_epirelentropy[w_dim, :]
    end
    # use nonlinear fit for higher dimensions
    rtw_dim = sqrt(w_dim)
    if w_dim <= 20
        u = 1.2023 / rtw_dim - 0.015
        v = 0.432 / rtw_dim + 1.0125
        w = -0.3057 / rtw_dim + 0.972
    else
        u = 1.1513 / rtw_dim - 0.0069
        v = 0.4873 / rtw_dim + 1.0008
        w = -0.4247 / rtw_dim + 0.9961
    end
    return [u, v, w]
end

const central_rays_epirelentropy = [
    0.827838399	1.290927714	0.805102005;
    0.708612491	1.256859155	0.818070438;
    0.622618845	1.231401008	0.829317079;
    0.558111266	1.211710888	0.838978357;
    0.508038611	1.196018952	0.847300431;
    0.468039614	1.183194753	0.854521307;
    0.435316653	1.172492397	0.860840992;
    0.408009282	1.163403374	0.866420017;
    0.38483862	1.155570329	0.871385499;
    0.364899122	1.148735192	0.875838068;
    ]


text:
abstract type vdW1fRuleModel <: MixingRule end

struct vdW1fRuleParam <: EoSParam
end

@newmodelsimple vdW1fRule vdW1fRuleModel vdW1fRuleParam

export vdW1fRule
function vdW1fRule(components::Vector{String}; activity=nothing, userlocations::Vector{String}=String[], activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    acentricfactor = SingleParam(params["w"],"acentric factor")
    packagedparams = vdW1fRuleParam()
    model = vdW1fRule(packagedparams, verbose=verbose)
    return model
end

function mixing_rule(model::ABCubicModel,V,T,z,mixing_model::vdW1fRuleModel,α,a,b,c)
    n = sum(z)
    invn = (one(n)/n)
    invn2 = invn^2
    #b̄ = dot(z,Symmetric(b),z) * invn2
    ā = zero(T+first(z))
    b̄ = zero(first(z))
    for i in 1:length(z)
        zi = z[i]
        αi = α[i]
        zi2 = zi^2
        b̄ += b[i,i]*zi2
        ā += a[i,i]*αi*zi^2
        for j in 1:(i-1)
            zij = zi*z[j]
            ā += 2*a[i,j]*sqrt(αi*α[j])*zij
            b̄ += 2*b[i,j]*zij
        end
    end
    ā *= invn2
    b̄ *= invn2
    c̄ = dot(z,c)*invn
    #dot(z,Symmetric(a .* sqrt.(α*α')),z) * invn2
    return ā,b̄,c̄
end

is_splittable(::vdW1fRule) = false