text:
# This files stores constraints or concepts that are not yet available within usual_constraints
# The list is based on https://arxiv.org/abs/2009.00514

# SECTION - Generic constraints
# intension: predicate over variables (DONE)
# extension: supports or conflicts with a set of configurations

# SECTION - Constraints defined from Languages
# regular: recognize a regular automaton
# mdd: follow a multi-values decision diagram

# SECTION - Comparison-based Constraints
# all_different: done in Constraints.jl
# candidate: all_different_except (DONE)
# all_equal: done in Constraints.jl
# all_equal_param: done (not in XCSP3)
# ordered: increasing (current ordered in Constraints.jl), strictly_increasing (DONE), decreasing, strictly_decreasing

# SECTION - Counting and Summing Constraints
# sum or linear or linear_sum: sum_equal_param (DONE)
# candidate: more generic linear equation (DONE)
# count: count, among, at_least, at_most
# n_values: n_values, at_least_n_values, at_most_n_values
# candidate: n_values_except (ignores the values in except)
# cardinality: (global_)cardinality

# SECTION - Connection Constraints
# maximum: all values lower than param (DONE)

# minimum
concept_minimum(x; param) = minimum(x) < param

const _minimum = Constraint(
    concept = concept_minimum,
    error = make_error(:minimum),
)

push!(BENCHED_CONSTRAINTS, :minimum => _minimum)

# element: param ∈ values(variables)
# channel: if xᵢ = paramⱼ then xⱼ = paramᵢ

# SECTION - Packing and Scheduling Constraints
# no_overlap: disjunctive (no_overlap_1D), diffn (kD)
# cummulative: check pdf :p

# SECTION - Constraints on Graphs
# circuit: verify that the sequence of vertices forms a circuit (where xᵢ = j => (i,j)∈E)

# SECTION - Elementary Constraints
# instantiation: x == param

## no overlap 1D and 2D (commented)
function concept_no_overlap(x; param)
    for i in 1:(length(x)-1), j in i:length(x)
        x[i] + param > x[j] && x[j] + param > x[i] && return false
    end
    return true
end

# function concept_no_overlap(x; param::AbstractVector)
#     for i in 1:(length(x)-1), j in i:length(x)
#         x[i] + param[i] > x[j] && x[j] + param[j] > x[i] && return false
#     end
#     return true
# end

const no_overlap = Constraint(
    concept = concept_no_overlap,
    error = make_error(:no_overlap),
)

push!(BENCHED_CONSTRAINTS, :no_overlap => no_overlap)

## maximum
concept_maximum(x; param) = maximum(x) < param

const _maximum = Constraint(
    concept = concept_maximum,
    error = make_error(:maximum),
)

push!(BENCHED_CONSTRAINTS, :maximum => _maximum)

## intension
concept_intension(x; param::Base.Function) = param(x)

const _intension = Constraint(
    concept = concept_intension,
    error = make_error(:intension),
)

push!(BENCHED_CONSTRAINTS, :intension => _intension)

## linear
concept_linear(x1, op1::Base.Function, x2, op2::Base.Function, param) = op2(op1(x1,x2),param)

const _linear = Constraint(
    concept = concept_linear,
    error = make_error(:intension),
)

push!(BENCHED_CONSTRAINTS, :linear => _linear)

## all_different_except

concept_all_different_except(x, s) = allunique(filter(element -> element ∉ s , x))

const _all_different_except = Constraint(
    concept = concept_all_different_except,
    error = make_error(:all_different_except),
)

push!(BENCHED_CONSTRAINTS, :all_different_except => _all_different_except)

## strictly_increasing

concept_strictly_increasing(x; param=nothing) = issorted(x) && length(x) == length(Set(x))

const _strictly_increasing = Constraint(
    concept = concept_strictly_increasing,
    error = make_error(:strictly_increasing),
)

push!(BENCHED_CONSTRAINTS, :strictly_increasing => _strictly_increasing)


text:
@testset "hyperPA" begin
    hyper_degrees = AliasTable([2345, 30991, 41226, 29829, 15690, 8247, 3983, 2279, 1669, 1158, 927, 1114, 384, 324, 292, 629])
    degrees = AliasTable([65, 32, 28, 2, 1, 9, 4, 6, 0, 29, 10, 11, 2, 27, 3, 0, 120, 4, 114, 3])
    @test hyper_pa(degrees, hyper_degrees, 16, Int32(1000)) isa Vector{Vector{Int32}}
end

ext: jl

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

# Porting note:
# There is no need for a `Thing` type because it's an abstract thing that could just be
# an entity, relation, of attribute. See `instantiate` function below.

function instantiate(t::Proto.Thing)
    if t._type.encoding == Proto.Type_Encoding.ENTITY_TYPE
        return Entity(t)
    elseif t._type.encoding == Proto.Type_Encoding.RELATION_TYPE
        return Relation(t)
    elseif t._type.encoding == Proto.Type_Encoding.ATTRIBUTE_TYPE
        return Attribute(t)
    else
        throw(TypeDBClientException(CONCEPT_BAD_ENCODING))
    end
end


text:
doc"""
  `Ubound{ESS,FSS}` creates a ubound, which represents an interval on the real line
  which cannot be expressed using a single unum.  Several constructors are available:

  `Ubound{ESS,FSS}(::Ubound)` copies the ubound.

  `Ubound{ESS,FSS}(left::Unum, right::Unum)`
  `Ubound{ESS,FSS}(left::Ubound, right::Unum)`
  `Ubound{ESS,FSS}(left::Unum, right::Ubound)`
  `Ubound{ESS,FSS}(left::Ubound, right::Ubound)`

  create a Ubound with the unums representing the outer hull of the two values.
  Left must be strictly less than right.

  In all cases  Ubound-passed parameters will copy their unums, but unum passed
  parameters will become "owned" by the Ubound.

  These constructors are unsafe when options[:devmode] is unset.  Guaranteed safe
  constructors are provided using the equivalent ubound() functions, these constructors
  will copy their passed unums.
"""
abstract Ubound{ESS,FSS} <: Real

type UboundSmall{ESS,FSS} <: Ubound{ESS,FSS}
  lower::UnumSmall{ESS,FSS}
  upper::UnumSmall{ESS,FSS}
end

#overload the call dispatches for naked ubound types.
@universal function (::Type{Ubound{ESS,FSS}})(x::Unum, y::Unum)
  B(x, y)
end
#an empty constructor defaults to the extended real line.
function (::Type{Ubound{ESS,FSS}}){ESS,FSS}()
  if FSS < 7
    UboundSmall{ESS,FSS}(neg_inf(UnumSmall{ESS,FSS}), pos_inf(UnumSmall{ESS,FSS}))
  else
    UboundLarge{ESS,FSS}(neg_inf(UnumLarge{ESS,FSS}), pos_inf(UnumLarge{ESS,FSS}))
  end
end

type UboundLarge{ESS,FSS} <: Ubound{ESS,FSS}
  lower::UnumLarge{ESS,FSS}
  upper::UnumLarge{ESS,FSS}
end

#to make Ubound construction easier, you don't have to specify the {ESS,FSS} pair when calling it.
(::Type{Ubound}){ESS,FSS}(lower::Unum{ESS,FSS}, upper::Unum{ESS,FSS}) = Ubound{ESS,FSS}(lower, upper)

@universal Base.copy(x::Ubound) = B(copy(x.lower), copy(x.upper))

export Ubound

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
mutable struct Ports <: DynDevice
    state::Vector{Symbol}
    inner::Vector
end

function _make_state_mapping(components::Vector{DT}, states::Vector{Symbol}) where DT <: DynamicComponent

    local_state_mapping = Dict{DT, Vector{Int64}}()
    local_state_space_ix = 0
    for c in components
        component_state_index= Vector{Int64}(undef, length(c.states))
        for (ix,i) in enumerate(c.states)
            local_state_space_ix += 1
            component_state_index[ix] = local_state_space_ix
        end
        local_state_mapping[c] = component_state_index
    end

    return local_state_mapping

end

function _make_port_mapping(components::Vector{DT}, states::Vector{Symbol}) where DT <: DynamicComponent

    input_port_mapping = Dict{DT, Vector{Int64}}()

    for c in components
        index_component_in = Vector{Int64}()
        for i in c.ports.state
            tmp = [(ix,var) for (ix,var) in enumerate(states) if var == i]
            isempty(tmp) && continue
            push!(index_component_in, tmp[1][1])
        end
        input_port_mapping[c] = index_component_in
    end

    return input_port_mapping

end


text:
# ---
# title: 104. Maximum Depth of Binary Tree
# id: problem104
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/maximum-depth-of-binary-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _its maximum depth_.
#
# A binary tree's **maximum depth**  is the number of nodes along the longest
# path from the root node down to the farthest leaf node.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/26/tmp-tree.jpg)
#
#
#
#     Input: root = [3,9,20,null,null,15,7]
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: root = [1,null,2]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: root = []
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: root = [0]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 104]`.
#   * `-100 <= Node.val <= 100`
#
#
## @lc code=start
using LeetCode

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


text:
#####
##### Momentum advection operators
#####

"""
    div_Uu(i, j, k, grid, advection, U, u)

Calculate the advection of momentum in the x-direction using the conservative form, ∇·(Uu)

    1/Vᵘ * [δxᶠᵃᵃ(ℑxᶜᵃᵃ(Ax * u) * ℑxᶜᵃᵃ(u)) + δy_fca(ℑxᶠᵃᵃ(Ay * v) * ℑyᵃᶠᵃ(u)) + δz_fac(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(u))]

which will end up at the location `fcc`.
"""
@inline function div_Uu(i, j, k, grid, advection, U, u)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, momentum_flux_uu, advection, U.u, u) +
                                    δyᵃᶜᵃ(i, j, k, grid, momentum_flux_uv, advection, U.v, u) +
                                    δzᵃᵃᶜ(i, j, k, grid, momentum_flux_uw, advection, U.w, u))
end

"""
    div_Uv(i, j, k, grid, advection, U, v)

Calculate the advection of momentum in the y-direction using the conservative form, ∇·(Uv)

    1/Vʸ * [δx_cfa(ℑyᵃᶠᵃ(Ax * u) * ℑxᶠᵃᵃ(v)) + δyᵃᶠᵃ(ℑyᵃᶜᵃ(Ay * v) * ℑyᵃᶜᵃ(v)) + δz_afc(ℑxᶠᵃᵃ(Az * w) * ℑzᵃᵃᶠ(w))]

which will end up at the location `cfc`.
"""
@inline function div_Uv(i, j, k, grid, advection, U, v)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, momentum_flux_vu, advection, U.u, v) +
                                    δyᵃᶠᵃ(i, j, k, grid, momentum_flux_vv, advection, U.v, v)    +
                                    δzᵃᵃᶜ(i, j, k, grid, momentum_flux_vw, advection, U.w, v))
end

"""
    div_Uw(i, j, k, grid, advection, U, w)

Calculate the advection of momentum in the z-direction using the conservative form, ∇·(Uw)

    1/Vʷ * [δx_caf(ℑzᵃᵃᶠ(Ax * u) * ℑxᶠᵃᵃ(w)) + δy_acf(ℑzᵃᵃᶠ(Ay * v) * ℑyᵃᶠᵃ(w)) + δzᵃᵃᶠ(ℑzᵃᵃᶜ(Az * w) * ℑzᵃᵃᶜ(w))]

which will end up at the location `ccf`.
"""
@inline function div_Uw(i, j, k, grid, advection, U, w)
    return 1/Vᵃᵃᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, momentum_flux_wu, advection, U.u, w) +
                                    δyᵃᶜᵃ(i, j, k, grid, momentum_flux_wv, advection, U.v, w) +
                                    δzᵃᵃᶠ(i, j, k, grid, momentum_flux_ww, advection, U.w, w))
end


text:
using SlurmHyperopt, DataFrames
using Test

# not a good test yet (it really only checks if the code is running, not if it is doing anything useful)
@testset "SlurmHyperopt.jl" begin

    extra_calls = "echo \"------------------------------------------------------------\"
    echo \"SLURM JOB ID: \$SLURM_JOBID\"
    echo \"\$SLURM_NTASKS tasks\"
    echo \"------------------------------------------------------------\"

    module load julia/1.7.0"

    julia_call = "julia test.jl \$SLURM_JOB_NAME \$SLURM_TASK_ID"

    slurm_file = "test_script.sh"

    params = SlurmParams(qos="short",
                    job_name="test",
                    account="test-account",
                    nodes=1,
                    ntasks_per_node=1,
                    extra_calls=extra_calls,
                    julia_call=julia_call,
                    file_path=slurm_file)

    N_jobs = 10
    samp = RandomSampler(a = LinRange(1,5,100), b = [true, false], c = exp10.(LinRange(-1,3,100)))
    ho = SlurmHyperoptimizer(N_jobs, samp, params)

    # fill the results with some values
    ress = rand(N_jobs)
    for i ∈ 1:N_jobs

        pars = ho[i]
        res = ress[i]

        SlurmHyperopt.save_result(ho, HyperoptResults(pars=pars, res=res), i)
    end

    merge_results!(ho)

    pars = get_params(ho) # get all the parameters
    res = get_results(ho) # get only the results

    @test res == ress
    df = DataFrame(ho)

    @test typeof(df) <: DataFrame
    @test size(df) == (10,4)
    @test df[!,"results"] == res
end


text:
#=
types.jl
Types and constructors for LinearMapAA objects.
2026-06-27 Jeff Fessler, University of Michigan
=#

using LinearMaps: LinearMap

LMAAkeys = (:_lmap, :_prop, :_idim, :_odim) # reserved

#=
Note: this old way may not properly allow `setindex!` to work as desired
because it may change the type of the lmap and of the prop:
`struct LinearMapAA{T, M <: LinearMap, P <: NamedTuple} <: AbstractMatrix{T}`
{T, M <: LinearMap, P <: NamedTuple}
    _lmap::M
    _prop::P
    function LinearMapAA{T}(L::M, p::P) where {T, M <: LinearMap, P <: NamedTuple}
    function LinearMapAA(L::LinearMap, p::NamedTuple) # where {T, M <: LinearMap, P <: NamedTuple}
   #    new{T,M,P}(L, p)
        new(L, p)
    end
=#

"""
    mutable struct LinearMapAM{T,Do,Di} <: AbstractMatrix{T}
"matrix" version that is quite akin to a matrix in its behavior
"""
mutable struct LinearMapAM{T,Do,Di} <: AbstractMatrix{T}
    _lmap::LinearMap # "L"
    _prop::NamedTuple # user-defined "named properties" accessible via A.name
    _idim::Dims{Di} # "input" dimensions, always (size(L,2),)
    _odim::Dims{Do} # "output" dimensions, always (size(L,1),)
end

"""
    struct LinearMapAO{T,Do,Di}
"tensor" version that can map from arrays to arrays
(it is not a subtype of `AbstractArray` and `setindex` is unsupported)
"""
struct LinearMapAO{T,Do,Di}
    _lmap::LinearMap # "L"
    _prop::NamedTuple # user-defined "named properties" accessible via A.name
    _idim::Dims{Di} # "input" dimensions, often (size(L,2),)
    _odim::Dims{Do} # "output" dimensions, often (size(L,1),)
end

# most operations apply to both AM and AO types:
LinearMapAX{T,Do,Di} =
    Union{ LinearMapAM{T,Do,Di}, LinearMapAO{T,Do,Di} } where {T,Do,Di}

# constructors

"""
    B = LinearMapAO{T,Do,Di}(A::LinearMapAX)
Make an AO from an AM, despite `idim` and `odim` being 1D,
for expert users who want `B*X` to be an `Array`.
Somewhat an opposite of `undim`.
"""
LinearMapAO(A::LinearMapAX{T,Do,Di}) where {T,Do,Di} =
    LinearMapAO{T,Do,Di}(A._lmap, A._prop, A._idim, A._odim)

"""
    A = LinearMapAA(L::LinearMap ; ...)

Constructor

options
- `prop::NamedTuple = NamedTuple()`
- `T = eltype(L)`
- `idim::Dims = (size(L,2),)`
- `odim::Dims = (size(L,1),)`
- `operator::Bool` by default: `false` if both `idim` & `odim` are 1D.

`prop` cannot include the fields `_lmap`, `_prop`, `_idim`, `_odim`

Output `A` is `LinearMapAO` if `operator` is `true`, else `LinearMapAM`.
"""
function LinearMapAA(L::LinearMap ;
    prop::NamedTuple = NamedTuple(),
    T::Type = eltype(L),
    idim::Dims{Di} = (size(L,2),),
    odim::Dims{Do} = (size(L,1),),
    operator::Bool = length(idim) > 1 || length(odim) > 1,
) where {Di,Do}

    size(L,2) == prod(idim) ||
        throw(DimensionMismatch("size2=$(size(L,2)) vs idim=$idim"))
    size(L,1) == prod(odim) ||
        throw(DimensionMismatch("size1=$(size(L,1)) vs odim=$odim"))
    length(intersect(propertynames(prop), LMAAkeys)) > 0 &&
        throw("invalid property field among $(propertynames(prop))")

    return operator ?
         LinearMapAO{T,Do,Di}(L, prop, idim, odim) :
         LinearMapAM{T,Do,Di}(L, prop, idim, odim)
end

# for backwards compatibility:
LinearMapAA(L, prop::NamedTuple ; kwargs...) =
    LinearMapAA(L::LinearMap ; prop=prop, kwargs...)

"""
    A = LinearMapAA(L::AbstractMatrix ; ...)
Constructor
"""
LinearMapAA(L::AbstractMatrix ; kwargs...) =
    LinearMapAA(LinearMap(L) ; kwargs...)

_ismutating(f) = first(methods(f)).nargs == 3

"""
    A = LinearMapAA(f::Function, fc::Function, D::Dims{2} [, prop::NamedTuple)]
    ; T::DataType = Float32, idim::Dims, odim::Dims)
Constructor
"""
function LinearMapAA(f::Function, fc::Function, D::Dims{2} ;
    T::DataType = Float32,
    idim::Dims = (D[2],),
    odim::Dims = (D[1],),
    kwargs...,
)

    if idim == (D[2],) && odim == (D[1],)
        fnew = f
        fcnew = fc
    else
        fnew = _ismutating(f) ?
            (y,x) -> f(reshape(y,odim), reshape(x,idim)) :
            x -> reshape(f(reshape(x,idim)), odim)
        fcnew = _ismutating(fc) ?
            (x,y) -> fc(reshape(x,idim), reshape(y,odim)) :
            y -> reshape(fc(reshape(y,odim)), idim)
    end
    LinearMapAA(LinearMap{T}(fnew, fcnew, D[1], D[2]) ;
        idim=idim, odim=odim, kwargs...)
end

LinearMapAA(f::Function, fc::Function, D::Dims{2}, prop::NamedTuple ; kwargs...) =
    LinearMapAA(f, fc, D ; prop=prop, kwargs...)

"""
    A = LinearMapAA(f::Function, D::Dims{2} [, prop::NamedTuple]; kwargs...)
Constructor
"""
function LinearMapAA(f::Function, D::Dims{2} ;
    T::DataType = Float32,
    idim::Dims = (D[2],),
    odim::Dims = (D[1],),
    kwargs...,
)

    if idim == (D[2],) && odim == (D[1],)
        fnew = f
    else
        fnew = _ismutating(f) ?
            (y,x) -> f(reshape(y,odim), reshape(x,idim)) :
            x -> reshape(f(reshape(x,idim)), odim)
    end
    LinearMapAA(LinearMap{T}(fnew, D[1], D[2]) ;
        idim=idim, odim=odim, kwargs...)
end

LinearMapAA(f::Function, D::Dims{2}, prop::NamedTuple ; kwargs...) =
    LinearMapAA(f, D ; prop=prop, kwargs...)


text:
# ---
# title: 840. Magic Squares In Grid
# id: problem840
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/magic-squares-in-grid/description/>
# hidden: true
# ---
#
# A `3 x 3` magic square is a `3 x 3` grid filled with distinct numbers
# **from**`1` **to**`9` such that each row, column, and both diagonals all have
# the same sum.
#
# Given a `row x col` `grid` of integers, how many `3 x 3` "magic square"
# subgrids are there?  (Each subgrid is contiguous).
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/11/magic_main.jpg)
#
#
#
#     Input: grid = [[4,3,8,4],[9,5,1,9],[2,7,6,2]]
#     Output: 1
#     Explanation:
#     The following subgrid is a 3 x 3 magic square:
#     ![](https://assets.leetcode.com/uploads/2020/09/11/magic_valid.jpg)
#     while this one is not:
#     ![](https://assets.leetcode.com/uploads/2020/09/11/magic_invalid.jpg)
#     In total, there is only one magic square inside the given grid.
#
#
# **Example 2:**
#
#
#
#     Input: grid = [[8]]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: grid = [[4,4],[3,3]]
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: grid = [[4,7,8],[9,5,1],[2,3,6]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `row == grid.length`
#   * `col == grid[i].length`
#   * `1 <= row, col <= 10`
#   * `0 <= grid[i][j] <= 15`
#
#
## @lc code=start
using LeetCode

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


text:
export nullspacex

"""
`nullspacex(A)` returns an exact basis for the matrix `A`
"""
function nullspacex(A::AbstractMatrix{T}) where {T}
    r, c = size(A)
    B = rrefx(A)

    leads = Int[]
    # in each row, find first 1
    for i = 1:r
        row = B[i, :]
        if all(row .== 0)
            continue
        end
        k = findfirst(row .!= 0)
        append!(leads, k)
    end

    frees = setdiff(collect(1:c), leads)

    result = Matrix{T}(undef, c, 0)
    for f in frees
        v = zeros(T, c)
        v[f] = T(1)
        for i = 1:length(leads)
            l = leads[i]  # (i,l) is a leading 1
            if l < f
                v[l] = -B[i, f]
            end
        end
        result = [result v]
    end
    return result
end

function nullspacex(A::AbstractMatrix{T}) where {T<:IntegerX}
    return nullspacex(big.(A) // 1)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Miscellaneous functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _get_significand_and_base(number::Number[, new_decimal_base::Integer])

Return the significand and the base of the `number` in base 10. If the parameter
`new_decimal_base` is passed, then the output base is converted to this value.
"""
function _get_significand_and_base(number::Number)
    if abs(number) > 0
        base        = floor(log10(abs(number)))
        significand = number / 10^base
        return significand, base
    else
        return number, 0
    end
end

function _get_significand_and_base(number::Number, new_decimal_base::Integer)
    if abs(number) > 0
        significand, base = _get_significand_and_base(number)
        fact = 10^(new_decimal_base - base)
        significand /=  fact
        return significand, new_decimal_base
    else
        return number, new_decimal_base
    end
end


text:
pprint(it...;sep=" ",ed="\n")=length(it)>0 ? (@inbounds for i in 1:length(it)-1;print(it[i],sep);end;print(last(it),ed)) : print(ed)

function main()
    a,b,c=parse.(Int64, split(readline()))
    count = 0
    while true
        if count>32
            pprint(-1)
            break
        end
        if a%2+b%2+c%2>0
            pprint(count)
            break
        end
        a,b,c=b/2+c/2,a/2+c/2,b/2+c/2
        count+=1
    end
end
main()

max_stars_repo_name: sugitanishi/competitive-programming

max_issues_repo_name: sugitanishi/competitive-programming

max_forks_repo_name: sugitanishi/competitive-programming


text:
using Test
using WaveProp
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Element iterator 1d" begin
    d = ReferenceLine()
    h = 0.1
    mesh = CartesianMesh(d;h=h)
    els = ElementIterator(mesh)
    @test length(els) == 10
    @test eltype(els) == HyperRectangle{1,Float64}
    @test first(els) == HyperRectangle((0,),(h,))
    for (n,el) in enumerate(els)
        @test el ≈ HyperRectangle((n-1)*h,n*h)
    end
end

@testset "Element iterator 2d" begin
    # d    = ReferenceSquare()
    # h    = 0.1
    # mesh = CartesianMesh(d;h=(h,h),n=(1,1))
    # els = ElementIterator(mesh)
    @test_broken length(els) == 100
    @test_broken eltype(els) == HyperRectangle{2,Float64}
    @test_broken first(els) == HyperRectangle((0,0),(h,h))
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
module Day9

using Advent: save_input

function _solve_a(nums, preamble=25)
    for i in (preamble + 1):lastindex(nums)
        found = any(Iterators.product(i - 25:i, i - 25:i)) do (j, k)
            nums[i] == nums[j] + nums[k] && j != k
        end
        found || return nums[i]
    end
    big(-1)
end

function solve_a(filename)
    nums = map(x -> parse(BigInt, x), eachline(filename))
    _solve_a(nums)
end

function solve_b(filename)
    nums = map(x -> parse(BigInt, x), eachline(filename))
    target = _solve_a(nums)
    for i in eachindex(nums), j in i + 2:lastindex(nums)
        @views sum(nums[i:j]) == target && return sum(extrema(nums[i:j]))
    end
    big(-1)
end

solve(filename = save_input(2020, 9)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
Base.Vector(v::Sym, n) = [v[i] for i in 1:n]
Base.Matrix(A::Sym, n, m) = [A[i,j] for i in 1:n, j in 1:m]


text:
# diff_map
SUITE["diff_map"] = BenchmarkGroup()
let sz=(2^7, 2^7+1)
    x = randn(prod(sz))

    T2d = diff2d_map_old(sz...)
    Tnd = MIRT.diff_map(sz)
    d = T2d * x

    SUITE["diff_map"]["diff_forward"] = BenchmarkGroup()
    SUITE["diff_map"]["diff_forward"]["baseline"] = @benchmarkable $T2d * $x
    SUITE["diff_map"]["diff_forward"]["MIRT"] = @benchmarkable $Tnd * $x

    SUITE["diff_map"]["diff_adjoint"] = BenchmarkGroup()
    SUITE["diff_map"]["diff_adjoint"]["baseline"] = @benchmarkable $T2d' * $d
    SUITE["diff_map"]["diff_adjoint"]["MIRT"] = @benchmarkable $Tnd' * $d
end

# finite difference
SUITE["finite difference"] = BenchmarkGroup()
let sz = (2^10, 2^9), dims=2
    x = rand(sz...);
    h = diff_forw(x, dims=dims)
    g = diffl(x, dims)
    out = Array{Float64, length(sz)}(undef, sz...)

    SUITE["finite difference"]["diff"] = @benchmarkable diff($x, dims=$dims)
    SUITE["finite difference"]["diffl"] = @benchmarkable diffl($x, $dims)
    SUITE["finite difference"]["diffl!"] = @benchmarkable diffl!($out, $x, $dims)
    SUITE["finite difference"]["diff_forw"] = @benchmarkable diff_forw($x, dims=$dims)

    SUITE["finite difference"]["diff_adj"] = @benchmarkable diff_adj($h, $sz; dims=$dims)
    SUITE["finite difference"]["diffl_adj"] = @benchmarkable diffl_adj($g, $dims)
    SUITE["finite difference"]["diffl_adj!"] = @benchmarkable diffl_adj!($out, $g, $dims)
end


text:
""" Bipartite graphs, directed and undirected, as C-sets.

A graph theorist might call these "bipartitioned graphs," as in a graph
*equipped with* a bipartition, as opposed to "bipartite graphs," as in a graph
that *can be* bipartitioned. Here we use the former notion, which is more
natural from the structuralist perspective, but the latter terminology, which is
shorter and more familiar.
"""
module BipartiteGraphs
export HasBipartiteVertices, HasBipartiteGraph,
  AbstractUndirectedBipartiteGraph, UndirectedBipartiteGraph,
  AbstractBipartiteGraph, BipartiteGraph, nv₁, nv₂, vertices₁, vertices₂,
  ne₁₂, ne₂₁, edges₁₂, edges₂₁, src₁, src₂, tgt₁, tgt₂,
  add_vertex₁!, add_vertex₂!, add_vertices₁!, add_vertices₂!,
  rem_vertex₁!, rem_vertex₂!, rem_vertices₁!, rem_vertices₂!,
  add_edge₁₂!, add_edge₂₁!, add_edges₁₂!, add_edges₂₁!,
  rem_edge₁₂!, rem_edge₂₁!, rem_edges₁₂!, rem_edges₂₁!

using ...Present, ...CSetDataStructures
using ..BasicGraphs: flatten
using ...Theories
import ..BasicGraphs: nv, ne, vertices, edges, src, tgt, add_edge!, add_edges!,
  rem_edge!, rem_edges!, inneighbors, outneighbors

# Base types
############

""" Abstract type for C-sets that contain a bipartition of vertices.
"""
@abstract_acset_type HasBipartiteVertices

""" Abstract type for C-sets that contain a (directed) bipartite graph.
"""
@abstract_acset_type HasBipartiteGraph <: HasBipartiteVertices

# Undirected bipartite graphs
#############################

@present TheoryUndirectedBipartiteGraph(FreeSchema) begin
  (V₁, V₂)::Ob
  E::Ob
  src::Hom(E, V₁)
  tgt::Hom(E, V₂)
end

""" Abstract type for undirected bipartite graphs.
"""
@abstract_acset_type AbstractUndirectedBipartiteGraph <: HasBipartiteVertices

""" An undirected bipartite graph.

It is a matter of perspective whether to consider such graphs "undirected," in
the sense that the edges have no orientation, or "unidirected," in the sense
that all edges go from vertices of type 1 to vertices of type 2.
"""
@acset_type UndirectedBipartiteGraph(TheoryUndirectedBipartiteGraph, index=[:src, :tgt]) <:
  AbstractUndirectedBipartiteGraph

function (::Type{T})(nv₁::Int, nv₂::Int) where T <: HasBipartiteVertices
  g = T()
  add_parts!(g, :V₁, nv₁)
  add_parts!(g, :V₂, nv₂)
  return g
end

""" Number of vertices of type 1 in a bipartite graph.
"""
nv₁(g::HasBipartiteVertices) = nparts(g, :V₁)

""" Number of vertices of type 2 in a bipartite graph.
"""
nv₂(g::HasBipartiteVertices) = nparts(g, :V₂)

""" Vertices of type 1 in a bipartite graph.
"""
vertices₁(g::HasBipartiteVertices) = parts(g, :V₁)

""" Vertices of type 2 in a bipartite graph.
"""
vertices₂(g::HasBipartiteVertices) = parts(g, :V₂)

nv(g::HasBipartiteVertices) = (nv₁(g), nv₂(g))
vertices(g::HasBipartiteVertices) = (vertices₁(g), vertices₂(g))

""" Add a vertex of type 1 to a bipartite graph.
"""
add_vertex₁!(g::HasBipartiteVertices; kw...) = add_part!(g, :V₁; kw...)

""" Add a vertex of type 2 to a bipartite graph.
"""
add_vertex₂!(g::HasBipartiteVertices; kw...) = add_part!(g, :V₂; kw...)

""" Add vertices of type 1 to a bipartite graph.
"""
add_vertices₁!(g::HasBipartiteVertices, n::Int; kw...) =
  add_parts!(g, :V₁, n; kw...)

""" Add vertices of type 2 to a bipartite graph.
"""
add_vertices₂!(g::HasBipartiteVertices, n::Int; kw...) =
  add_parts!(g, :V₂, n; kw...)

""" Remove vertex of type 1 from a bipartite graph.
"""
rem_vertex₁!(g::HasBipartiteVertices, v::Int; kw...) =
  rem_vertices₁!(g, v:v; kw...)

""" Remove vertex of type 2 from a bipartite graph.
"""
rem_vertex₂!(g::HasBipartiteVertices, v::Int; kw...) =
  rem_vertices₂!(g, v:v; kw...)

""" Remove vertices of type 1 from a bipartite graph.
"""
function rem_vertices₁!(g::AbstractUndirectedBipartiteGraph, vs;
                        keep_edges::Bool=false)
  if !keep_edges
    rem_parts!(g, :E, unique!(sort!(flatten(incident(g, vs, :src)))))
  end
  rem_parts!(g, :V₁, vs)
end

""" Remove vertices of type 2 from a bipartite graph.
"""
function rem_vertices₂!(g::AbstractUndirectedBipartiteGraph, vs;
                        keep_edges::Bool=false)
  if !keep_edges
    rem_parts!(g, :E, unique!(sort!(flatten(incident(g, vs, :tgt)))))
  end
  rem_parts!(g, :V₂, vs)
end

ne(g::AbstractUndirectedBipartiteGraph) = nparts(g, :E)
edges(g::AbstractUndirectedBipartiteGraph) = parts(g, :E)
edges(g::AbstractUndirectedBipartiteGraph, src::Int, tgt::Int) =
  (e for e in incident(g, src, :src) if subpart(g, e, :tgt) == tgt)
src(g::AbstractUndirectedBipartiteGraph, args...) = subpart(g, args..., :src)
tgt(g::AbstractUndirectedBipartiteGraph, args...) = subpart(g, args..., :tgt)

add_edge!(g::AbstractUndirectedBipartiteGraph, src::Int, tgt::Int; kw...) =
  add_part!(g, :E; src=src, tgt=tgt, kw...)

function add_edges!(g::AbstractUndirectedBipartiteGraph,
                    srcs::AbstractVector{Int}, tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  add_parts!(g, :E, n; src=srcs, tgt=tgts, kw...)
end

rem_edge!(g::AbstractUndirectedBipartiteGraph, e::Int) = rem_part!(g, :E, e)
rem_edge!(g::AbstractUndirectedBipartiteGraph, src::Int, tgt::Int) =
  rem_edge!(g, first(edges(g, src, tgt)))
rem_edges!(g::AbstractUndirectedBipartiteGraph, es) = rem_parts!(g, :E, es)

inneighbors(g::AbstractUndirectedBipartiteGraph, v::Int) =
  subpart(g, incident(g, v, :tgt), :src)
outneighbors(g::AbstractUndirectedBipartiteGraph, v::Int) =
  subpart(g, incident(g, v, :src), :tgt)

# Bipartite graphs
##################

@present TheoryBipartiteGraph(FreeSchema) begin
  (V₁, V₂)::Ob
  (E₁₂, E₂₁)::Ob
  src₁::Hom(E₁₂, V₁)
  tgt₂::Hom(E₁₂, V₂)
  src₂::Hom(E₂₁, V₂)
  tgt₁::Hom(E₂₁, V₁)
end

""" Abstract type for bipartite graphs.
"""
@abstract_acset_type AbstractBipartiteGraph <: HasBipartiteGraph

""" A bipartite graph, better known as a [bipartite directed
multigraph](https://cs.stackexchange.com/a/91521).

Directed bipartite graphs are isomorphic to port hypergraphs and to whole grain
Petri nets.
"""
@acset_type BipartiteGraph(TheoryBipartiteGraph, index=[:src₁, :src₂, :tgt₁, :tgt₂]) <:
  AbstractBipartiteGraph

""" Number of edges from V₁ to V₂ in a bipartite graph.
"""
ne₁₂(g::HasBipartiteGraph) = nparts(g, :E₁₂)
ne₁₂(g::HasBipartiteGraph, src::Int, tgt::Int) =
  count(subpart(g, e, :tgt₂) == tgt for e in incident(g, src, :src₁))

""" Number of edges from V₂ to V₁ in a bipartite graph.
"""
ne₂₁(g::HasBipartiteGraph) = nparts(g, :E₂₁)
ne₂₁(g::HasBipartiteGraph, src::Int, tgt::Int) =
  count(subpart(g, e, :tgt₁) == tgt for e in incident(g, src, :src₂))

""" Edges from V₁ to V₂ in a bipartite graph.
"""
edges₁₂(g::HasBipartiteGraph) = parts(g, :E₁₂)
edges₁₂(g::HasBipartiteGraph, src::Int, tgt::Int) =
  (e for e in incident(g, src, :src₁) if subpart(g, e, :tgt₂) == tgt)

""" Edges from V₂ to V₁ in a bipartite graph.
"""
edges₂₁(g::HasBipartiteGraph) = parts(g, :E₂₁)
edges₂₁(g::HasBipartiteGraph, src::Int, tgt::Int) =
  (e for e in incident(g, src, :src₂) if subpart(g, e, :tgt₁) == tgt)

ne(g::HasBipartiteGraph) = (ne₁₂(g), ne₂₁(g))
edges(g::HasBipartiteGraph) = (edges₁₂(g), edges₂₁(g))

""" Source vertex of edge from V₁ to V₂ in a bipartite graph.
"""
src₁(g::HasBipartiteGraph, args...) = subpart(g, args..., :src₁)

""" Target vertex of edge from V₁ to V₂ in a bipartite graph.
"""
tgt₂(g::HasBipartiteGraph, args...) = subpart(g, args..., :tgt₂)

""" Source vertex of edge from V₂ to V₁ in a bipartite graph.
"""
src₂(g::HasBipartiteGraph, args...) = subpart(g, args..., :src₂)

""" Target vertex of edge from V₂ to V₁ in a bipartite graph.
"""
tgt₁(g::HasBipartiteGraph, args...) = subpart(g, args..., :tgt₁)

""" Add edge from V₁ to V₂ in a bipartite graph.
"""
add_edge₁₂!(g::HasBipartiteGraph, src::Int, tgt::Int; kw...) =
  add_part!(g, :E₁₂; src₁=src, tgt₂=tgt, kw...)

""" Add edge from V₂ to V₁ in a bipartite graph.
"""
add_edge₂₁!(g::HasBipartiteGraph, src::Int, tgt::Int; kw...) =
  add_part!(g, :E₂₁; src₂=src, tgt₁=tgt, kw...)

""" Add edges from V₁ to V₂ in a bipartite graph.
"""
function add_edges₁₂!(g::HasBipartiteGraph, srcs::AbstractVector{Int},
                      tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  add_parts!(g, :E₁₂, n; src₁=srcs, tgt₂=tgts, kw...)
end

""" Add edges from V₂ to V₁ in a bipartite graph.
"""
function add_edges₂₁!(g::HasBipartiteGraph, srcs::AbstractVector{Int},
                      tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  add_parts!(g, :E₂₁, n; src₂=srcs, tgt₁=tgts, kw...)
end

function rem_vertices₁!(g::AbstractBipartiteGraph, vs; keep_edges::Bool=false)
  if !keep_edges
    rem_parts!(g, :E₁₂, unique!(sort!(flatten(incident(g, vs, :src₁)))))
    rem_parts!(g, :E₂₁, unique!(sort!(flatten(incident(g, vs, :tgt₁)))))
  end
  rem_parts!(g, :V₁, vs)
end

function rem_vertices₂!(g::AbstractBipartiteGraph, vs; keep_edges::Bool=false)
  if !keep_edges
    rem_parts!(g, :E₂₁, unique!(sort!(flatten(incident(g, vs, :src₂)))))
    rem_parts!(g, :E₁₂, unique!(sort!(flatten(incident(g, vs, :tgt₂)))))
  end
  rem_parts!(g, :V₂, vs)
end

""" Remove edge from V₁ to V₂ in a bipartite graph.
"""
rem_edge₁₂!(g::AbstractBipartiteGraph, e::Int) = rem_part!(g, :E₁₂, e)
rem_edge₁₂!(g::AbstractBipartiteGraph, src::Int, tgt::Int) =
  rem_edge₁₂!(g, first(edges₁₂(g, src, tgt)))

""" Remove edge from V₁ to V₂ in a bipartite graph.
"""
rem_edge₂₁!(g::AbstractBipartiteGraph, e::Int) = rem_part!(g, :E₂₁, e)
rem_edge₂₁!(g::AbstractBipartiteGraph, src::Int, tgt::Int) =
  rem_edge₂₁!(g, first(edges₂₁(g, src, tgt)))

""" Remove edges from V₁ to V₂ in a bipartite graph.
"""
rem_edges₁₂!(g::AbstractBipartiteGraph, es) = rem_parts!(g, :E₁₂, es)

""" Remove edges from V₂ to V₁ in a bipartite graph.
"""
rem_edges₂₁!(g::AbstractBipartiteGraph, es) = rem_parts!(g, :E₂₁, es)

end


text:
"""
    AbstractFilterRepresentation{S}

Abstract type for representation of the conditional distribution over the hidden state of type `S`.
"""
abstract type AbstractFilterRepresentation{S} end

"""
    represented_type(rep::AbstractFilterRepresentation)

Return the type which is represented by `rep`.
"""
represented_type(rep::AbstractFilterRepresentation{S}) where S = S

"""
    dim(rep::AbstractFilterRepresentation{S})

Return the dimensionality of the filter representation `rep`.
"""
function dim(rep::AbstractFilterRepresentation{S}) where S end

# How to add new filter representations:
# * Add a struct which is a subtype of AbstractFilterRepresentation{S}, where S is the type to be represented
# * Implement a method for dim which returns the dimensionality of the representation.

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
is_interval_constraint(con_ref::ConstraintRef) = constraint_object(con_ref).set isa MOI.Interval
is_interval_constraint(con_ref::NonlinearConstraintRef) = count(i -> i == :(<=), Meta.parse(string(con_ref)).args) == 2
is_equality_constraint(con_ref::ConstraintRef) = constraint_object(con_ref).set isa MOI.EqualTo
is_equality_constraint(con_ref::NonlinearConstraintRef) = Meta.parse(string(con_ref)).args[1] == :(==)
JuMP.name(con_ref::NonlinearConstraintRef) = ""

"""
    check_constraint!(m::Model, constr::Tuple)

Check constraints in a disjunction Tuple.

    check_constraint!(m::Model, constr_j, constr_list::Vector)

Check nested constraint and update `constr_list`.

    check_constraint!(m::Model, constr)

Check constraint in a Model.

    check_constraint!(m::Model, constr::Nothing)

Return nothing for an empty disjunct.
"""
function check_constraint!(m::Model, constr::Tuple)
    constr_list = []
    map(constr_j -> check_constraint!(m, constr_j, constr_list), constr)
    return Tuple(constr_list)
end
function check_constraint!(m::Model, constr_j::Tuple, constr_list::Vector)
    map(constr_jk -> check_constraint!(m, constr_jk, constr_list), constr_j)
end
function check_constraint!(m::Model, constr_j::AbstractArray{<:ConstraintRef}, constr_list::Vector)
    push!(constr_list, check_constraint!(m, constr_j))
end
function check_constraint!(m::Model, constr_j::ConstraintRef, constr_list::Vector)
    push!(constr_list, check_constraint!(m, constr_j))
end
function check_constraint!(m::Model, constr::ConstraintRef)
    @assert all(is_valid(m, constr)) "$constr is not a valid constraint."
    new_constr = split_constraint(m, constr)
    if isnothing(new_constr)
        new_constr = constr
    else
        constr_name = gen_constraint_name(constr)
        m[constr_name] = new_constr
        constr_str = split(string(constr),"}")[end]
        @warn "$constr_str uses the `MOI.Interval` or `MOI.EqualTo` set. Each instance of the interval set has been split into two constraints, one for each bound."
        delete_original_constraint!(m, constr)
    end
    return new_constr
end
function check_constraint!(m::Model, constr::AbstractArray{<:ConstraintRef})
    @assert all(is_valid.(m, constr)) "$constr is not a valid constraint."
    if !any(is_interval_constraint.(constr)) && !any(is_equality_constraint.(constr))
        new_constr = constr
    else
        idxs = get_indices(constr)
        constr_dict = Dict(union(
            [
                split_constraint(m, constr[idx...]) |>
                    i -> isnothing(i) ?
                        (idx...,"") => constr[idx...] :
                        [(idx...,"lb") => i[1], (idx...,"ub") => i[2]]
                for idx in idxs
            ]...
        ))
        new_constr = Containers.SparseAxisArray(constr_dict)
        constr_name = gen_constraint_name(constr)
        m[constr_name] = new_constr
        constr_str = split(string(constr),"}")[end]
        @warn "$constr_str uses the `MOI.Interval` or `MOI.EqualTo` set. Each instance of the interval set has been split into two constraints, one for each bound."
        delete_original_constraint!(m, constr)
    end
    return new_constr
end
check_constraint!(m::Model, constr::Nothing) = nothing

"""
    split_constraint(m::Model, constr::NonlinearConstraintRef)

Split a nonlinear constraint that is an Interval or EqualTo constraint.

    split_constraint(m::Model, constr::ConstraintRef, constr_name::String = name(constr))

Split a linear or quadratic constraint.

    split_constraint(m::Model, constr_obj::ScalarConstraint, lb_name::String, ub_name::String)

Split a constraint that is a MOI.EqualTo or MOI.Interval.

    split_constraint(m::Model, func::Union{AffExpr,QuadExpr}, lb::Float64, ub::Float64, lb_name::String, ub_name::String)

Create split constraint for linear or quadratic constraint.

    split_constraint(m::Model, constr::ConstraintRef, constr_func_expr::Expr, lb::Float64, ub::Float64)

Split a nonlinear constraint.

    split_constraint(args...)

Return nothing for an empty disjunct.
"""
function split_constraint(m::Model, constr::NonlinearConstraintRef)
    constr_expr = Meta.parse(string(constr))
    if constr_expr.args[1] == :(==) #replace == for lb <= expr <= ub and split
        lb, ub = 0, 0#rhs is always 0, but could get obtained from: constr_expr.args[3]
        constr_func_expr = copy(constr_expr.args[2])
        new_constraints = split_constraint(m, constr, constr_func_expr, lb, ub)
        return new_constraints
    elseif count(x -> x == :(<=), constr_expr.args) == 2 #split lb <= expr <= ub
        lb = constr_expr.args[1]
        ub = constr_expr.args[5]
        constr_func_expr = copy(constr_expr.args[3]) #get func part of constraint
        new_constraints = split_constraint(m, constr, constr_func_expr, lb, ub)
        return new_constraints
    else
        return nothing
    end
end
function split_constraint(m::Model, constr::ConstraintRef, constr_name::String = name(constr))
    if isempty(constr_name)
        constr_name = "[$constr]"
    end
    lb_name = name_split_constraint(constr_name, :lb)
    ub_name = name_split_constraint(constr_name, :ub)
    constr_obj = constraint_object(constr)
    new_constraints = split_constraint(m, constr_obj, lb_name, ub_name)

    return new_constraints
end
function split_constraint(m::Model, func::Union{AffExpr,QuadExpr}, lb::Float64, ub::Float64, lb_name::String, ub_name::String)
    return [
        @constraint(m, lb <= func, base_name = lb_name),
        @constraint(m, func <= ub, base_name = ub_name)
    ]
end
function split_constraint(m::Model, constr_obj::ScalarConstraint{T,<:MOI.EqualTo}, lb_name::String, ub_name::String) where T
    split_constraint(m, constr_obj.func, constr_obj.set.value, constr_obj.set.value, lb_name, ub_name)
end
function split_constraint(m::Model, constr_obj::ScalarConstraint{T,<:MOI.Interval}, lb_name::String, ub_name::String) where T
    split_constraint(m, constr_obj.func, constr_obj.set.lower, constr_obj.set.upper, lb_name, ub_name)
end
function split_constraint(m::Model, constr::ConstraintRef, constr_func_expr::Expr, lb::Float64, ub::Float64)
    replace_JuMPvars!(constr_func_expr, m) #replace Expr with JuMP vars
    #replace original constraint with lb <= func
    lb_constr = JuMP._NonlinearConstraint(
        JuMP._NonlinearExprData(m, constr_func_expr),
        lb,
        Inf
    )
    m.nlp_data.nlconstr[constr.index.value] = lb_constr
    #create new constraint for func <= ub
    constr_expr_ub = Expr(:call, :(<=), constr_func_expr, ub)
    ub_constr = add_nonlinear_constraint(m, constr_expr_ub)

    return [constr, ub_constr]
end
split_constraint(args...) = nothing

delete_original_constraint!(m::Model, constr::ConstraintRef) = delete(m,constr)
delete_original_constraint!(m::Model, constr::NonlinearConstraintRef) = nothing
delete_original_constraint!(m::Model, constr::AbstractArray{<:ConstraintRef}) = map(c -> delete_original_constraint!(m,c), constr)

"""
    parse_constraint(constr::ConstraintRef)

Extract constraint operator symbol (op), constraint function expression (LHS), and constraint RHS.
"""
function parse_constraint(constr::ConstraintRef)
    #check function has a single comparrison operator (<=, >=, ==)
    constr_str = replace(split(string(constr),": ")[end], " " => "", "²" => "^2") #remove name, blanks and power notation
    constr_expr = Meta.parse(constr_str).args
    op = constr_expr[1] #comparrison operator
    lhs = constr_expr[2] #LHS of the constraint
    rhs = constr_expr[3] #RHS of constraint
    return op, lhs, rhs
end

"""

"""
function replace_constraint(constr::NonlinearConstraintRef, sym_expr, op, rhs)
    #convert symbolic function to expression
    expr = Base.remove_linenums!(build_function(sym_expr)).args[2].args[1]
    #replace symbolic variables by their JuMP variables and math operators with their symbols
    m = constr.model
    replace_JuMPvars!(expr, m)
    replace_operators!(expr)
    # determine bounds of original constraint (if op is ==, both bounds are set to rhs)
    upper_b = (op == :(>=)) ? Inf : rhs
    lower_b = (op == :(<=)) ? -Inf : rhs
    # replace NL constraint currently in the model with the reformulated one
    m.nlp_data.nlconstr[constr.index.value] = JuMP._NonlinearConstraint(JuMP._NonlinearExprData(m, expr), lower_b, upper_b)
end
function replace_constraint(constr::ConstraintRef, sym_expr, op, rhs)
    #convert symbolic function to expression
    op = eval(op)
    expr = Base.remove_linenums!(build_function(op(sym_expr,rhs))).args[2].args[1]
    #replace symbolic variables by their JuMP variables and math operators with their symbols
    m = constr.model
    replace_JuMPvars!(expr, m)
    replace_operators!(expr)
    #add a nonlinear constraint with the perspective function and delete old constraint
    add_nonlinear_constraint(m, expr)
    delete(m, constr)
end


text:
type Lungs
    Rp::Float64
    Ra::Vector{Float64}
    Rv::Vector{Float64}
    Ca::Vector{Float64}
    Cv::Vector{Float64}
    La::Float64
    Lv::Float64
    V0::Vector{Float64}
    Pp::Vector{Float64}
    Pa::Any
    Va::Any
    Qa::Any
    Pv::Any
    Vv::Any
    Qv::Any

    function Lungs()
        this = new()
        this.Rp = 0.0251*mmHgToPa/cm3Tom3
        this.Ra = [0.0227*mmHgToPa/cm3Tom3];
        push!(this.Ra,0.03*mmHgToPa/cm3Tom3);
        push!(this.Ra,0.021*mmHgToPa/cm3Tom3);
        this.Rv = [0.010*mmHgToPa/cm3Tom3];
        push!(this.Rv,0.010*mmHgToPa/cm3Tom3);
        this.Ca = [2.222*cm3Tom3/mmHgToPa];
        push!(this.Ca,1.481*cm3Tom3/mmHgToPa);
        push!(this.Ca,1.778*cm3Tom3/mmHgToPa);
        this.Cv = [13*cm3Tom3/mmHgToPa];
        push!(this.Cv,74*cm3Tom3/mmHgToPa);
        this.La = 5e-5*mmHgToPa/cm3Tom3;
        this.Lv = 5e-5*mmHgToPa/cm3Tom3;
        this.V0 = [10*cm3Tom3];
        append!(this.V0,[10,10,10,10]*cm3Tom3);
        this.Pp = Vector{Float64}[];
        this.Pa = Array{Float64,2}[];
        this.Va = Array{Float64,2}[];
        this.Qa = Array{Float64,2}[];
        this.Pv = Array{Float64,2}[];
        this.Vv = Array{Float64,2}[];
        this.Qv = Array{Float64,2}[];
        return this
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
# random functions that dispatch either to CURAND or GPUArrays' generic RNG

using Random

export rand_logn!, rand_poisson!

# the interface is split in two levels:
# - functions that extend the Random standard library, and take an RNG as first argument,
#   will only ever dispatch to CURAND and as a result are limited in the types they support.
# - functions that take an array will dispatch to either CURAND or GPUArrays
# - non-exported functions are provided for constructing GPU arrays from only an eltype

const curand_rng = CURAND.default_rng
gpuarrays_rng() = GPUArrays.default_rng(CuArray)

function seed!(seed=Base.rand(UInt64))
    Random.seed!(curand_rng(), seed)
    Random.seed!(gpuarrays_rng(), seed)
end

# CURAND in-place
Random.rand!(A::CURAND.UniformArray) = Random.rand!(curand_rng(), A)
Random.randn!(A::CURAND.NormalArray; kwargs...) = Random.randn!(curand_rng(), A; kwargs...)
rand_logn!(A::CURAND.LognormalArray; kwargs...) = CURAND.rand_logn!(curand_rng(), A; kwargs...)
rand_poisson!(A::CURAND.PoissonArray; kwargs...) = CURAND.rand_poisson!(curand_rng(), A; kwargs...)

# CURAND out-of-place
rand(T::CURAND.UniformType, dims::Dims) = Random.rand(curand_rng(), T, dims)
randn(T::CURAND.NormalType, dims::Dims; kwargs...) = Random.randn(curand_rng(), T, dims; kwargs...)
rand_logn(T::CURAND.LognormalType, dims::Dims; kwargs...) = CURAND.rand_logn(curand_rng(), T, dims; kwargs...)
rand_poisson(T::CURAND.PoissonType, dims::Dims; kwargs...) = CURAND.rand_poisson(curand_rng(), T, dims; kwargs...)

# support all dimension specifications
rand(T::CURAND.UniformType, dim1::Integer, dims::Integer...) =
    Random.rand(curand_rng(), T, Dims((dim1, dims...)))
randn(T::CURAND.NormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(curand_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_logn(T::CURAND.LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_logn(curand_rng(), T, Dims((dim1, dims...)); kwargs...)
rand_poisson(T::CURAND.PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_poisson(curand_rng(), T, Dims((dim1, dims...)); kwargs...)

# GPUArrays in-place
Random.rand!(A::CuArray) = Random.rand!(gpuarrays_rng(), A)
Random.randn!(A::CuArray; kwargs...) =
    error("CUDA.jl does not support generating normally-distributed random numbers of type $(eltype(A))")
# FIXME: GPUArrays.jl has a randn! nowadays, but it doesn't work with e.g. Cuint
rand_logn!(A::CuArray; kwargs...) =
    error("CUDA.jl does not support generating lognormally-distributed random numbers of type $(eltype(A))")
rand_poisson!(A::CuArray; kwargs...) =
    error("CUDA.jl does not support generating Poisson-distributed random numbers of type $(eltype(A))")

# GPUArrays out-of-place
rand(T::Type, dims::Dims) = Random.rand!(CuArray{T}(undef, dims...))
randn(T::Type, dims::Dims; kwargs...) = Random.randn!(CuArray{T}(undef, dims...); kwargs...)
rand_logn(T::Type, dims::Dims; kwargs...) = rand_logn!(CuArray{T}(undef, dims...); kwargs...)
rand_poisson(T::Type, dims::Dims; kwargs...) = rand_poisson!(CuArray{T}(undef, dims...); kwargs...)

# support all dimension specifications
rand(T::Type, dim1::Integer, dims::Integer...) =
    Random.rand!(CuArray{T}(undef, dim1, dims...))
randn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_logn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_poisson(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson!(CuArray{T}(undef, dim1, dims...); kwargs...)

# untyped out-of-place
rand(dim1::Integer, dims::Integer...) =
    Random.rand(curand_rng(), Dims((dim1, dims...)))
randn(dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(curand_rng(), Dims((dim1, dims...)); kwargs...)
rand_logn(dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_logn(curand_rng(), Dims((dim1, dims...)); kwargs...)
rand_poisson(dim1::Integer, dims::Integer...; kwargs...) =
    CURAND.rand_poisson(curand_rng(), Dims((dim1, dims...)); kwargs...)


text:
using CryptoTools
using Test
using Random

@testset "Shamir" begin
    scheme = ShamirSecretSharing()
    secret = BigInt(1234931)
    shares = split(scheme, secret, 6, 3)
    @test join(scheme, shares) == secret
    @test join(scheme, shares[1:3]) == secret
    @test join(scheme, shares[4:6]) == secret
end


text:
#=
Model: Beam (348 variables, 1 input)
=#
using Reachability, MAT, SparseArrays

beam(o::Pair{Symbol, <:Any}...) = beam(Options(Dict{Symbol, Any}(o)))

function beam(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    file = matopen(@relpath "beam.mat")
    A = read(file, "A")

    # initial set
    # - x1-x300 are 0,
    # - the rest is in [0.002, 0.0015]
    X0 = Hyperrectangle([zeros(300); fill(0.00175, 48)],
                        [zeros(300); fill(0.00025, 48)])

    # input set
    B = read(file, "B")
    U = BallInf([0.5], 0.3)

    # instantiate continuous LTI system
    S = InitialValueProblem(
        ConstrainedLinearControlContinuousSystem(A, B, nothing, U), X0)

    # safety property: x89 < 2100
    property = SafeStatesProperty(HalfSpace(sparsevec([89], [1.0], 348), 2100.))

    # =======================
    # Problem default options
    # =======================
    partition = [(2*i-1:2*i) for i in 1:174] # 2D blocks

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [89],
                                  :partition => partition,
                                  :plot_vars => [0, 89])

    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => [89],
                                  :partition => partition,
                                  :property => property)
    end

    return (S, merge(problem_options, input_options))
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer

import SimpleDirectMediaLayer.LoadBMP

SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLEBUFFERS, 16)
SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLESAMPLES, 16)

SDL2.init()

win = SDL2.CreateWindow("Hello World!", Int32(100), Int32(100), Int32(800), Int32(600),
    UInt32(SDL2.WINDOW_SHOWN))
SDL2.SetWindowResizable(win,true)

renderer = SDL2.CreateRenderer(win, Int32(-1),
    UInt32(SDL2.RENDERER_ACCELERATED | SDL2.RENDERER_PRESENTVSYNC))

import Base.unsafe_convert
unsafe_convert(::Type{Ptr{SDL2.RWops}}, s::String) = SDL2.RWFromFile(s, "rb")

LoadBMP(src::String) = SDL2.LoadBMP_RW(src,Int32(1))

bkg = SDL2.Color(200, 200, 200, 255)

# Create text
font = TTF_OpenFont(joinpath(@__DIR__,"../../assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf"), 14)
#txt = "@BinDeps.install Dict([ (:glib, :libglib) ])"
txt = "Hello, World!"
text = TTF_RenderText_Blended(font, txt, SDL2.Color(20,20,20,255))
tex = SDL2.CreateTextureFromSurface(renderer,text)

fx,fy = Int[1], Int[1]
TTF_SizeText(font, txt, pointer(fx), pointer(fy))
fx,fy = fx[1],fy[1]

#img = SDL2.LoadBMP("LB2951.jpg")
#tex = SDL2.CreateTextureFromSurface(ren, img)
#SDL2.FreeSurface(img)

function pollEvent!()
    #SDL2.Event() = [SDL2.Event(NTuple{56, Uint8}(zeros(56,1)))]
    SDL_Event() = Array{UInt8}(zeros(56))
    e = SDL_Event()
    success = (SDL2.PollEvent(e) != 0)
    return e,success
end
function getEventType(e::Array{UInt8})
    # HAHA This is still pretty janky, but I guess that's all you can do w/ unions.
    bitcat(UInt32, e[4:-1:1])
end
function getEventType(e::SDL2.Event)
    e._Event[1]
end

function bitcat(outType::Type{T}, arr)::T where T<:Number
    out = zero(outType)
    for x in arr
        out = out << (sizeof(x)*8)
        out |= convert(T, x)  # the `convert` prevents signed T from promoting to Int64.
    end
    out
end

gameIsRunning = true

while gameIsRunning
    x,y = Int[1], Int[1]
    SDL2.PumpEvents()
    SDL2.GetMouseState(pointer(x), pointer(y))

    # Set render color to red ( background will be rendered in this color )
    SDL2.SetRenderDrawColor(renderer, 200, 200, 200, 255)
    SDL2.RenderClear(renderer)

    SDL2.SetRenderDrawColor(renderer, 20, 50, 105, 255)
    SDL2.RenderDrawLine(renderer,0,0,800,600)

    rect = SDL2.Rect(1,1,200,200)
    SDL2.RenderFillRect(renderer, pointer_from_objref(rect) )
    SDL2.RenderCopy(renderer, tex, C_NULL, pointer_from_objref(SDL2.Rect(x[1],y[1],fx,fy)))

    SDL2.RenderPresent(renderer)
    println("check")
    e,_ = pollEvent!()
           t = getEventType(e)
           if t == SDL2.QUIT
               println("quit")
               #throw(QuitException())
               gameIsRunning = false
           end
    sleep(0.01)

end
SDL2.Quit()

max_stars_repo_name: Kyjor/julia-game

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: Kyjor/julia-game

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: Kyjor/julia-game

max_forks_repo_licenses: MIT Unlicense


text:
using POMDPModels
using POMDPs
using POMDPModelTools
using StaticArrays
using Random
using DiscreteValueIteration
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using Test

# write your own tests here
include("runtests_versus_discrete_vi.jl")


text:
"""
Iso Periodic kernel
"""
struct IsoPeriodKernel{T, VT<:AbstractVector{T}} <: AbstractKernel
	lp::VT
	ll::VT
	lσ::VT
end
@functor IsoPeriodKernel
function (IsoPK::IsoPeriodKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	n = size(x, 2)
	scaled_x = x .* exp(-IsoPK.lp[1])
	d = norm2_metric(scaled_x)
	σ_square = exp(2*IsoPK.lσ[1])
	(@. σ_square*exp(-2.0*sinpi(d)*sinpi(d)*exp(-2*IsoPK.ll[1]))) + Diagonal(λ*ones(n))
end
function (IsoPK::IsoPeriodKernel{T, VT})(x, xo) where {T, VT}
	scaled_x = x .* exp(-IsoPK.lp[1])
	scaled_xo = xo .* exp(-IsoPK.lp[1])
	d = norm2_metric(scaled_x, scaled_xo)
	σ_square = exp(2*IsoPK.lσ[1])
	@. σ_square*exp(-2.0*sinpi(d)*sinpi(d)*exp(-2*IsoPK.ll[1]))
end

function IsoPeriodKernel()
	lp = 1.0 ./ rand(1)
	ll = rand(1)
	lσ = rand(1)
	IsoPeriodKernel(lp, ll, lσ)
end

reset(IsoPK::IsoPeriodKernel) = IsoPeriodKernel()

function Base.show(io::IO, kernel::IsoPeriodKernel)
	print(io, "IsoPeriodKernel(")
	print(io, "lp=", kernel.lp[1], ", ", "ll=", kernel.ll[1], ", ", "lσ=", kernel.lσ[1])
	print(io, ")")
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
push!(LOAD_PATH, "../src/Utils/");  push!(LOAD_PATH, "../src/Leaf/")

using MathToolsMod, Leaf, LeafPhotosynthesisMod, Plots

tmax = 10000 # Time in seconds here
l = leaf_params{Float32}()
f = LeafPhotosynthesisMod.fluxes{Float32}()
l.dynamic_state = true
l.gstyp = 1

uu = zeros(tmax,10)
uu[1,2] = 0.07
uu[1,1] = 0.5
f.APAR = 100
l.rdleaf = 0.0
for c = 2:1:tmax
   #println(c)a
    if c==1000 f.APAR = 1500 end
    if c==5000 f.APAR = 100 end
       l.gs = uu[c-1,2]
       l.Kn = uu[c-1,1]

       uu[c-1,3]=l.ϕs
       uu[c-1,4]=l.Ci
       uu[c-1,5]=f.φ
       uu[c-1,6]=l.Kp
       uu[c-1,8]=f.an
       uu[c-1,7]=f.APAR
       #l.Kp = 4
       LeafPhotosynthesisMod.LeafPhotosynthesis(f,l,298.0)
       uu[c,1] = uu[c-1,1]+(l.Kn_ss-l.Kn)/10*1/60
       uu[c,2] = uu[c-1,2]+(l.gs_ss-l.gs)/15*1/60
       #println(l.Kn_ss, " ",  l.Ci, " ", f.φ)
end

max_stars_repo_name: ruillercoaster/LSM-SPAM

max_issues_repo_name: ruillercoaster/LSM-SPAM

max_forks_repo_name: ruillercoaster/LSM-SPAM


text:
import Base: *, ∈, isempty

export CartesianProduct,
       CartesianProductArray,
       CartesianProduct!,
       array

"""
    CartesianProduct{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}

Type that represents a Cartesian product of two convex sets.

### Fields

- `X` -- first convex set
- `Y` -- second convex set

### Notes

The Cartesian product of three elements is obtained recursively.
See also `CartesianProductArray` for an implementation of a Cartesian product of
many sets without recursion, instead using an array.

The `EmptySet` is the absorbing element for `CartesianProduct`.

Constructors:

- `CartesianProduct{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}(X1::S1, X2::S2)`
  -- default constructor

- `CartesianProduct(Xarr::Vector{S}) where {S<:LazySet}`
  -- constructor from an array of convex sets
"""
struct CartesianProduct{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
    X::S1
    Y::S2
end

# constructor from an array
CartesianProduct(Xarr::Vector{S}) where {N<:Real, S<:LazySet{N}} =
    (length(Xarr) == 0
        ? EmptySet{N}()
        : length(Xarr) == 1
            ? Xarr[1]
            : length(Xarr) == 2
                ? CartesianProduct(Xarr[1], Xarr[2])
                : CartesianProduct(Xarr[1],
                                   CartesianProduct(Xarr[2:length(Xarr)])))

# EmptySet is the absorbing element for CartesianProduct
@absorbing(CartesianProduct, EmptySet)

"""
```
    *(X::LazySet, Y::LazySet)
```

Alias for the binary Cartesian product.
"""
*(X::LazySet, Y::LazySet) = CartesianProduct(X, Y)

"""
    ×

Alias for the binary Cartesian product.
"""
×(X::LazySet, Y::LazySet) = CartesianProduct(X, Y)

"""
    dim(cp::CartesianProduct)::Int

Return the dimension of a Cartesian product.

### Input

- `cp` -- Cartesian product

### Output

The ambient dimension of the Cartesian product.
"""
function dim(cp::CartesianProduct)::Int
    return dim(cp.X) + dim(cp.Y)
end

"""
    σ(d::AbstractVector{N}, cp::CartesianProduct{N}) where {N<:Real}

Return the support vector of a Cartesian product.

### Input

- `d`  -- direction
- `cp` -- Cartesian product

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the wrapped sets.
"""
function σ(d::AbstractVector{N}, cp::CartesianProduct{N}) where {N<:Real}
    n1 = dim(cp.X)
    return [σ(d[1:n1], cp.X); σ(d[n1+1:length(d)], cp.Y)]
end

"""
    ρ(d::AbstractVector{N}, cp::CartesianProduct{N}) where {N<:Real}

Return the support function of a Cartesian product.

### Input

- `d`  -- direction
- `cp` -- Cartesian product

### Output

The support function in the given direction.
If the direction has norm zero, the result depends on the wrapped sets.
"""
function ρ(d::AbstractVector{N}, cp::CartesianProduct{N}) where {N<:Real}
    n1 = dim(cp.X)
    return ρ(d[1:n1], cp.X) + ρ(d[n1+1:length(d)], cp.Y)
end

"""
    isbounded(cp::CartesianProduct)::Bool

Determine whether a Cartesian product is bounded.

### Input

- `cp` -- Cartesian product

### Output

`true` iff both wrapped sets are bounded.
"""
function isbounded(cp::CartesianProduct)::Bool
    return isbounded(cp.X) && isbounded(cp.Y)
end

"""
    ∈(x::AbstractVector{N}, cp::CartesianProduct{N})::Bool where {N<:Real}

Check whether a given point is contained in a Cartesian product.

### Input

- `x`  -- point/vector
- `cp` -- Cartesian product

### Output

`true` iff ``x ∈ cp``.
"""
function ∈(x::AbstractVector{N}, cp::CartesianProduct{N})::Bool where {N<:Real}
    @assert length(x) == dim(cp)

    n1 = dim(cp.X)
    return ∈(view(x, 1:n1), cp.X) &&
           ∈(view(x, n1+1:length(x)), cp.Y)
end

"""
    isempty(cp::CartesianProduct)::Bool

Return if a Cartesian product is empty or not.

### Input

- `cp` -- Cartesian product

### Output

`true` iff any of the sub-blocks is empty.
"""
function isempty(cp::CartesianProduct)::Bool
    return isempty(cp.X) || isempty(cp.Y)
end

"""
    constraints_list(cp::CartesianProduct{N}) where {N<:Real}

Return the list of constraints of a (polytopic) Cartesian product.

### Input

- `cp` -- Cartesian product

### Output

A list of constraints.
"""
function constraints_list(cp::CartesianProduct{N}) where {N<:Real}
    return constraints_list(CartesianProductArray([cp.X, cp.Y]))
end

"""
    vertices_list(cp::CartesianProduct{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a (polytopic) Cartesian product.

### Input

- `cp` -- Cartesian product

### Output

A list of vertices.

### Algorithm

We assume that the underlying sets are polytopic.
Then the high-dimensional set of vertices is just the Cartesian product of the
low-dimensional sets of vertices.
"""
function vertices_list(cp::CartesianProduct{N}
                      )::Vector{Vector{N}} where {N<:Real}
    # collect low-dimensional vertices lists
    vlist_low = (vertices_list(cp.X), vertices_list(cp.Y))

    # create high-dimensional vertices list
    vlist = Vector{Vector{N}}()
    m = length(vlist_low[1]) * length(vlist_low[2])
    sizehint!(vlist, m)
    for v1 in vlist_low[1]
        for v2 in vlist_low[2]
            push!(vlist, vcat(v1, v2))
        end
    end

    return vlist
end

# ======================================
#  Cartesian product of an array of sets
# ======================================
"""
    CartesianProductArray{N<:Real, S<:LazySet{N}} <: LazySet{N}

Type that represents the Cartesian product of a finite number of convex sets.

### Fields

- `array` -- array of sets

### Notes

The `EmptySet` is the absorbing element for `CartesianProductArray`.

Constructors:

- `CartesianProductArray(array::Vector{<:LazySet})` -- default constructor

- `CartesianProductArray([n]::Int=0, [N]::Type=Float64)`
  -- constructor for an empty product with optional size hint and numeric type
"""
struct CartesianProductArray{N<:Real, S<:LazySet{N}} <: LazySet{N}
    array::Vector{S}
end

# constructor for an empty product with optional size hint and numeric type
function CartesianProductArray(n::Int=0, N::Type=Float64)::CartesianProductArray
    arr = Vector{LazySet{N}}()
    sizehint!(arr, n)
    return CartesianProductArray(arr)
end

# EmptySet is the absorbing element for CartesianProductArray
@absorbing(CartesianProductArray, EmptySet)

# add functions connecting CartesianProduct and CartesianProductArray
@declare_array_version(CartesianProduct, CartesianProductArray)

"""
    array(cpa::CartesianProductArray{N, S}
         )::Vector{S} where {N<:Real, S<:LazySet{N}}

Return the array of a Cartesian product of a finite number of convex sets.

### Input

- `cpa` -- Cartesian product array

### Output

The array of a Cartesian product of a finite number of convex sets.
"""
function array(cpa::CartesianProductArray{N, S}
              )::Vector{S} where {N<:Real, S<:LazySet{N}}
    return cpa.array
end

"""
    dim(cpa::CartesianProductArray)::Int

Return the dimension of a Cartesian product of a finite number of convex sets.

### Input

- `cpa` -- Cartesian product array

### Output

The ambient dimension of the Cartesian product of a finite number of convex
sets.
"""
function dim(cpa::CartesianProductArray)::Int
    return length(cpa.array) == 0 ? 0 : sum([dim(Xi) for Xi in cpa.array])
end

"""
    σ(d::AbstractVector{N}, cpa::CartesianProductArray{N}) where {N<:Real}

Support vector of a Cartesian product array.

### Input

- `d`   -- direction
- `cpa` -- Cartesian product array

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the product sets.
"""
function σ(d::AbstractVector{N}, cpa::CartesianProductArray{N}) where {N<:Real}
    svec = similar(d)
    i0 = 1
    for Xi in cpa.array
        i1 = i0 + dim(Xi) - 1
        svec[i0:i1] = σ(d[i0:i1], Xi)
        i0 = i1 + 1
    end
    return svec
end

"""
    ρ(d::AbstractVector{N}, cp::CartesianProductArray{N}) where {N<:Real}

Return the support function of a Cartesian product array.

### Input

- `d`   -- direction
- `cpa` -- Cartesian product array

### Output

The support function in the given direction.
If the direction has norm zero, the result depends on the wrapped sets.
"""
function ρ(d::AbstractVector{N}, cpa::CartesianProductArray{N}) where {N<:Real}
    sfun = zero(N)
    i0 = 1
    for Xi in cpa.array
        i1 = i0 + dim(Xi) - 1
        sfun += ρ(d[i0:i1], Xi)
        i0 = i1 + 1
    end
    return sfun
end

"""
    isbounded(cpa::CartesianProductArray)::Bool

Determine whether a Cartesian product of a finite number of convex sets is
bounded.

### Input

- `cpa` -- Cartesian product of a finite number of convex sets

### Output

`true` iff all wrapped sets are bounded.
"""
function isbounded(cpa::CartesianProductArray)::Bool
    return all(x -> isbounded(x), cpa.array)
end

"""
    ∈(x::AbstractVector{N}, cpa::CartesianProductArray{N}
     )::Bool where {N<:Real}

Check whether a given point is contained in a Cartesian product of a finite
number of sets.

### Input

- `x`   -- point/vector
- `cpa` -- Cartesian product array

### Output

`true` iff ``x ∈ \\text{cpa}``.
"""
function ∈(x::AbstractVector{N}, cpa::CartesianProductArray{N}
          )::Bool where {N<:Real}
    @assert length(x) == dim(cpa)

    i0 = 1
    for Xi in cpa.array
        i1 = i0 + dim(Xi) - 1
        if !∈(x[i0:i1], Xi)
            return false
        end
        i0 = i1 + 1
    end
    return true
end

"""
    isempty(cpa::CartesianProductArray)::Bool

Return if a Cartesian product is empty or not.

### Input

- `cp` -- Cartesian product

### Output

`true` iff any of the sub-blocks is empty.
"""
function isempty(cpa::CartesianProductArray)::Bool
    return any(X -> isempty(X), array(cpa))
end

"""
    constraints_list(cpa::CartesianProductArray{N}) where {N<:Real}

Return the list of constraints of a (polytopic) Cartesian product of a finite
number of sets.

### Input

- `cpa` -- Cartesian product array

### Output

A list of constraints.
"""
function constraints_list(cpa::CartesianProductArray{N}) where {N<:Real}
    clist = Vector{LinearConstraint{N}}()
    n = dim(cpa)
    sizehint!(clist, n)
    prev_step = 1
    # create high-dimensional constraints list
    for c_low in array(cpa)
        c_low_list = constraints_list(c_low)
        if !isempty(c_low_list)
            indices = prev_step : (dim(c_low_list[1]) + prev_step - 1)
        end
        for constr in c_low_list
            new_constr = LinearConstraint(sparsevec(indices, constr.a, n),
                                          constr.b)
            push!(clist, new_constr)
        end
        prev_step += dim(c_low_list[1])
    end

    return clist
end

"""
    vertices_list(cpa::CartesianProductArray{N}
                 )::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a (polytopic) Cartesian product of a finite
number of sets.

### Input

- `cpa` -- Cartesian product array

### Output

A list of vertices.

### Algorithm

We assume that the underlying sets are polytopic.
Then the high-dimensional set of vertices is just the Cartesian product of the
low-dimensional sets of vertices.
"""
function vertices_list(cpa::CartesianProductArray{N}
                      )::Vector{Vector{N}} where {N<:Real}
    # collect low-dimensional vertices lists
    vlist_low = [vertices_list(X) for X in array(cpa)]

    # create high-dimensional vertices list
    indices_max = [length(vl) for vl in vlist_low]
    m = prod(indices_max)
    vlist = Vector{Vector{N}}(undef, m)
    indices = ones(Int, length(vlist_low))
    v = zeros(N, dim(cpa))
    dim_start_j = 1
    for vl in vlist_low
        v_low = vl[1]
        v[dim_start_j:dim_start_j+length(v_low)-1] = v_low
        dim_start_j += length(v_low)
    end
    i = 1
    j = 1
    # iterate through all index combinations
    while true
        indices[1] = 0
        while indices[1] < indices_max[1]
            indices[1] += 1
            v_low = vlist_low[1][indices[1]]
            v[1:length(v_low)] = v_low
            vlist[i] = copy(v)
            i += 1
        end
        if i > m
            break
        end
        j = 1
        dim_start_j = 1
        while indices[j] == indices_max[j]
            indices[j] = 1
            v_low = vlist_low[j][1]
            v[dim_start_j:dim_start_j+length(v_low)-1] = v_low
            dim_start_j += length(v_low)
            j += 1
        end
        indices[j] += 1
        v_low = vlist_low[j][indices[j]]
        v[dim_start_j:dim_start_j+length(v_low)-1] = v_low
    end

    return vlist
end

"""
    same_block_structure(x::AbstractVector{S1}, y::AbstractVector{S2}
                        )::Bool where {S1<:LazySet, S2<:LazySet}

Check whether two vectors of sets have the same block structure, i.e., the
``i``-th entry in the vectors have the same dimension.

### Input

- `x` -- first vector
- `y` -- second vector

### Output

`true` iff the vectors have the same block structure.
"""
function same_block_structure(x::AbstractVector{S1}, y::AbstractVector{S2}
                             )::Bool where {S1<:LazySet, S2<:LazySet}
    if length(x) != length(y)
        return false
    end
    for i in 1:length(x)
        if dim(x[i]) != dim(y[i])
            return false
        end
    end
    return true
end


text:
module PkgHelp

using PkgDev
using YAML
using Coverage

function expanduser_(path...)
    if is_unix()
        expanduser(joinpath(path...))
    else
        if length(path) > 1 && path[1] == "~"
            joinpath(homedir(), path[2:end]...)
        else
            joinpath(path...)
        end
    end
end

function get_config(fun::Function, variable::String)
    pkghelp_conf = expanduser_("~",".julia","pkghelp.yml")
    config = isfile(pkghelp_conf) ? YAML.load(open(pkghelp_conf)) : Dict()

    if variable in keys(config)
        return config[variable]
    end

    config[variable] = fun()

    open(pkghelp_conf, "w") do file
        for k in keys(config)
            write(file, "$(k): $(config[k])\n")
        end
    end

    config[variable]
end

function get_syncdir()
    get_config("dir") do
        f = ""
        while !isdir(f)
            print("Input Julia pkg sync dir: ")
            path = strip(readline(STDIN))
            path = split(path, contains(path, "/") ? "/" : "\\")
            f = expanduser_(path...)
        end
        f
    end
end

function get_readme()
    get_config("readme") do
        print("Change default README format from Markdown to Org? [Y/n] ")
        lowercase(strip(readline())) == "n" ? "markdown" : "org"
    end
end

function get_git_config(variable)
    try
        strip(readstring(`git config --get $(variable)`))
    catch
        ""
    end
end

function get_github_user()
    github_user = get_git_config("github.user")
    if github_user == ""
        print("Enter your Github username: ")
        github_user = strip(readline(STDIN))
        run(`git config --global github.user $(github_user)`)
    end
    github_user
end

function set_git_ssh(dir)
    cd(dir) do
        git_url = get_git_config("remote.origin.url")
        m = match(r"https://github.com/(.+)/(.+)", git_url)
        if m != nothing
            new_git_url = "git@github.com:$(m[1])/$(m[2])"
            println("Remote url $(git_url) -> $(new_git_url)")
            run(`git remote set-url origin $(new_git_url)`)
            run(`git config remote.origin.pushurl $(new_git_url)`)
        end
    end
end

function org_readme(dir, pkg_name)
    cd(dir) do
        open("README.org", "w") do file
            write(file, "#+TITLE: $(pkg_name).jl\n")
            user_name = get_git_config("user.name")
            user_name != "" && write(file, "#+AUTHOR: $(user_name)\n")
            user_email = get_git_config("user.email")
            user_email != "" && write(file, "#+EMAIL: $(user_email)\n")
            github_user = get_github_user()
            write(file, "\n")
            for (url,badge) in
                [("https://travis-ci.org/$(github_user)/$(pkg_name).jl",
                  "https://travis-ci.org/$(github_user)/$(pkg_name).jl.svg?branch=master"),
                 ("https://coveralls.io/github/$(github_user)/$(pkg_name).jl?branch=master",
                  "https://coveralls.io/repos/github/$(github_user)/$(pkg_name).jl/badge.svg?branch=master"),
                 ("http://codecov.io/gh/$(github_user)/$(pkg_name).jl",
                  "http://codecov.io/gh/$(github_user)/$(pkg_name).jl/branch/master/graph/badge.svg")]
                write(file, "[[$(url)][$(badge)]]\n")
            end
            write(file, "\n")
        end
        rm("README.md")
        run(`git add README.org`)
        run(`git commit -a -m "README.md -> README.org"`)
    end
end

function generate(pkg_name, license; org::Bool=false, kwargs...)
    ismatch(r"\.jl$", pkg_name) && (pkg_name = pkg_name[1:end-3])

    PkgDev.generate(pkg_name, license; kwargs...)

    pkg_dir = Pkg.dir(pkg_name)
    pkg_sync_dir = joinpath(get_syncdir(), "$(pkg_name).jl")

    get_github_user()
    set_git_ssh(pkg_dir)
    get_readme() == "org" && org_readme(pkg_dir, pkg_name)

    if org
        open(joinpath(pkg_dir, "src", "$(pkg_name).org"), "w") do file
            write(file, "#+TITLE: $(pkg_name).jl\n")
            user_name = get_git_config("user.name")
            user_name != "" && write(file, "#+AUTHOR: $(user_name)\n")
            user_email = get_git_config("user.email")
            user_email != "" && write(file, "#+EMAIL: $(user_email)\n")
            write(file, "\n#+PROPERTY: header-args:julia :session *julia-$(pkg_name)*\n")
        end

        open(joinpath(pkg_dir, "src", "$(pkg_name).jl"), "w") do file
            write(file, "module $(pkg_name)\n")
            write(file, "\n")
            write(file, "const codefile = joinpath(dirname(@__FILE__), \"literate_org_tangled_code.jl\")\n")
            write(file, "if isfile(codefile)\n")
            write(file, "    include(codefile)\n")
            write(file, "else\n")
            write(file, "    error(\"$(pkg_name) not properly installed. Please run Pkg.build(\\\"$(pkg_name)\\\") then restart Julia.\")\n")
            write(file, "end\n")
            write(file, "\n")
            write(file, "end # module\n")
        end

        open(joinpath(pkg_dir, "test", "runtests.jl"), "w") do file
            write(file, "using $(pkg_name)\n")
            write(file, "using Base.Test\n")
            write(file, "\n")
            write(file, "const testfile = joinpath(dirname(@__FILE__), \"literate_org_tangled_tests.jl\")\n")
            write(file, "if isfile(testfile)\n")
            write(file, "    include(testfile)\n")
            write(file, "else\n")
            write(file, "    error(\"$(pkg_name) not properly installed. Please run Pkg.build(\\\"$(pkg_name)\\\") then restart Julia.\")\n")
            write(file, "end\n")
        end

        open(joinpath(pkg_dir, "REQUIRE"), "r+") do file
            seekend(file)
            write(file, "\nLiterateOrg\n")
        end

        open(joinpath(pkg_dir, ".gitignore"), "r+") do file
            seekend(file)
            write(file, "\n/src/literate_org_tangled_code.jl\n/test/literate_org_tangled_tests.jl\n")
        end

        open(joinpath(pkg_dir, "README.org"), "r+") do file
            seekend(file)
            write(file, "\nThis is a [[https://github.com/jagot/LiterateOrg.jl][LiterateOrg.jl]] project. The documentation is found [[file:src/$(pkg_name).org][within the code]].\n")
        end

        mkpath(joinpath(pkg_dir, "deps"))
        open(joinpath(pkg_dir, "deps", "build.jl"), "w") do file
            write(file, "using LiterateOrg\n")
            write(file, "tangle_package(joinpath(Pkg.dir(\"$(pkg_name)\", \"src\", \"$(pkg_name).org\")), \"$(pkg_name)\")")
        end

        cd(pkg_dir) do
            run(`git add .`)
            run(`git commit -a -m "LiterateOrg project"`)
        end
    end

    println("Moving $(pkg_dir) -> $(pkg_sync_dir)")
    mv(pkg_dir, pkg_sync_dir)

    link(pkg_name)
end

function link(pkg_name)
    ismatch(r"\.jl$", pkg_name) && (pkg_name = pkg_name[1:end-3])

    pkg_dir = Pkg.dir(pkg_name)
    pkg_sync_dir = joinpath(get_syncdir(), "$(pkg_name).jl")

    println("Linking $(pkg_sync_dir) -> $(pkg_dir)")
    symlink(pkg_sync_dir, pkg_dir)
end

function test(pkg_name)
    ismatch(r"\.jl$", pkg_name) && (pkg_name = pkg_name[1:end-3])
    pkg_folder = joinpath(Pkg.dir(pkg_name), "src")
    clean_folder(dirname(pkg_folder))
    Pkg.test(pkg_name, coverage = true)
    coverage = process_folder(pkg_folder)
    coverage_file = joinpath(pkg_folder, "lcov.info")
    LCOV.writefile(coverage_file, coverage)
    try
        real_folder = readlink(dirname(pkg_folder))
        str = readstring(coverage_file)
        open(coverage_file, "w") do file
            write(file, replace(str, dirname(pkg_folder), real_folder))
        end
    catch
    end
    summary = get_summary(coverage)
    @printf("Coverage: %0.2f %%\n", 100summary[1]/summary[2])
end

function clone(github_repo)
    m = match(r"([a-zA-Z0-9]+)/([a-zA-Z0-9]+)$", github_repo)
    user_name = m[1]
    pkg_name = m[2]
    try
        Pkg.installed(pkg_name) != nothing && return
    catch
        Pkg.clone("https://github.com/$(github_repo).jl.git", pkg_name)
    end
    user_name == get_github_user() && set_git_ssh(Pkg.dir(pkg_name))
end

end # module


text:
using TeukolskyQNMFunctions
using Documenter

DocMeta.setdocmeta!(TeukolskyQNMFunctions, :DocTestSetup, :(using TeukolskyQNMFunctions); recursive=true)

makedocs(;
    modules=[TeukolskyQNMFunctions],
    authors="Justin L. Ripley",
    repo="https://github.com/JLRipley314/TeukolskyQNMFunctions.jl.git",
    sitename="TeukolskyQNMFunctions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JLRipley314.github.io/TeukolskyQNMFunctions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JLRipley314/TeukolskyQNMFunctions.jl",
    devbranch="main",
)


text:
travis = get(ENV, "TRAVIS", "") == "true"

import GLFW

if !travis
	include("windowclose.jl")
end


text:
using DReal

begin
reset_ctx!()

# Example
d = -0.45
g = 9.8

# x position
x = DReal.Var(Float64, 10.0, 15.0)

# velocity
v = DReal.Var(Float64, 10.0, 15.0)

# dx/dt =v
dx = DReal.ODE(x, v)
dv = DReal.ODE(v, -g + d*v)
t0 = DReal.Var(Float64, 0.0, 3.0)

x_0 = DReal.Var(Float64, 0.0, 15.0)
v_0 = DReal.Var(Float64, 0.0, 0.0)
x_t = DReal.Var(Float64, 0.0, 15.0)
v_t = DReal.Var(Float64, 0.0, 0.0)

zeroex = convert(DReal.Ex{Float64}, 0.0)

integral_constraint = DReal.integral([dx, dv], [x_t, v_t], zeroex, t0, [x_0, v_0])
println("Made integral constraint")

DReal.add!(integral_constraint)
println("added constraint")

res = DReal.is_satisfiable()
println("result is ", res)

end


text:
using Pkg
using Pkg.Artifacts
using SHA
using Tar

# we have an issue with unpack(), so wrote a custom download
# https://discourse.julialang.org/t/unable-to-automatically-install-artifact/51984/2
function download_brainflow_artifact()

    url = brainflow_url()

    download_path = "$(tempname())-download.tar"
    download(url, download_path)

    brainflow_hash = create_artifact() do artifact_dir
        # Pkg.PlatformEngines.unpack() gives errors for some users on Windows
        Tar.extract(download_path, artifact_dir)
    end

    rm(download_path)

    return brainflow_hash
end

# library path you can optionally use while developing
function dev_library_path()
    return abspath(joinpath(dirname(@__DIR__), "lib"))
end

function get_brainflow_artifact_path()
    artifacts_toml = find_artifacts_toml(@__DIR__) # is there a better way? @__DIR__ makes brainflow non-relocatable I believe.
    brainflow_hash = artifact_hash("brainflow", artifacts_toml)
    if isdir(dev_library_path())
        # developer library takes precedence
        println(string("Using local libraries: ", dev_library_path()))
        return dev_library_path()
    elseif artifact_exists(brainflow_hash)
        # libraries are automatically stored here for users of brainflow
        return artifact_path(brainflow_hash)
    else
        println("Downloading artifact: brainflow")
        brainflow_hash = download_brainflow_artifact()
        return artifact_path(brainflow_hash)
    end
end

function interface_path(library::AbstractString)
    return abspath(INTERFACE_PATH, interface_name(library))
end

function interface_name(library::AbstractString)
    if Sys.iswindows()
        return "$library.dll"
    elseif Sys.isapple()
        return "lib$library.dylib"
    else
        return "lib$library.so"
    end
end

const INTERFACE_PATH = get_brainflow_artifact_path()

const DATA_HANDLER_INTERFACE = interface_path("DataHandler")
const BOARD_CONTROLLER_INTERFACE = interface_path("BoardController")
const ML_MODULE_INTERFACE = interface_path("MLModule")


text:
module TestDet

import ConicHulls.Dets

for n=1:5
    for k=1:100
        A = rand(-8:8, (n,n))
        d0 = iround(det(A))
        d1 = Dets.det([A[:,k] for k=1:n]...)
        d2 = dot(Dets.incomplete_det([A[:,k] for k=1:(n-1)]...), A[:,end])

        @assert d1 == d0
        @assert d2 == d0
    end
end

end # module


text:
# <md>
# strong chimeras in two-cluster networks

# In[]
# loading necessary packages
#using DifferentialEquations
using Plots; pyplot()
using LightGraphs
#using GraphPlot
using LinearAlgebra
using Statistics
using BenchmarkTools
using ProgressMeter
using LaTeXStrings
using DelimitedFiles

# In[]
n = 3     # size of the clusters
c = .2    # intra- and inter-coupling strength ratio
r = 1.2   # parameter of the optoelectronic map
K = 1.2   # overall coupling strength

# In[]
# discrete map model of optoelectronic oscillators
function f(x)
    return sin(x+π/4)^2
end

# In[]
# reduce the two-cluster system to a two-oscillator system
# simulate the two-oscillator reduced system to study the two-cluster coherent state
T = 100000 # total simulation time
x = zeros(2,T) # oscillator states
x[:,1] = rand(2,1) # random initial condition
p = Progress(T, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow) # progress bar
@time for t in 2:T
    x[1,t] = r*f(x[1,t-1]) + c*n*K*(f(x[2,t-1])-f(x[1,t-1]))
    x[2,t] = r*f(x[2,t-1]) + c*n*K*(f(x[1,t-1])-f(x[2,t-1]))
    next!(p)
end

# In[]
# plot showing the states of the two oscillators in the reduced system
scatter(1:T,x[1,:],marker=(:circle,4,.6,:blue,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
scatter!(1:T,x[2,:],marker=(:circle,3,.6,:green,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
title!("Reduced System for Two-Cluster Coherent States")
xaxis!(L"t",font(15),(T-99,T))
yaxis!(L"x_i",font(15),(1,1.1))
#yticks!([-π,0,π],[L"\pi",L"0",L"\pi"])

# In[]
# plot showing the effective input from one cluster to the other
scatter(1:T,c*n*K*(f.(x[2,:])-f.(x[1,:])),marker=(:circle,3,.6,:red,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
title!("Reduced System for Two-Cluster Coherent States")
xaxis!(L"t",font(15),(T-99,T))
yaxis!("Input",font(15),(-.05,.05))

# In[]
# create the quotient network that corresponds to chimera states
function quotient_net(c)
    ring = watts_strogatz(n, 2, 0)
    adj = Matrix(adjacency_matrix(ring))
    return adj = [adj c*n*ones(n,1); c*ones(1,n) 0]
end

# In[]
# simulate the (n+1)-oscillator reduced system to study chimera states
function chimera(adj,T,IC)
    x = zeros(n+1,T) # oscillator states
    driving = zeros(T) # effective input from the incoherent cluster to the coherent cluster
    x[:,1] = IC # initial condition
    for t in 2:T
        for i in 1:n # updating oscillators in the incoherent cluster
            x[i,t] = r*f(x[i,t-1]) + K*adj[i,:]⋅(f.(x[:,t-1]).-f(x[i,t-1]))
        end
        driving[t] = K*adj[n+1,:]⋅(f.(x[:,t-1]).-f(x[n+1,t-1]))
        x[n+1,t] = r*f(x[n+1,t-1]) + driving[t] # updating the oscillator representing the coherent cluster
    end
    # return the mean and the s.t.d. of the effective input from the incoherent cluster to the coherent cluster, and the final oscillator states
    return mean(driving[T-10000:T]), std(driving[T-10000:T]), x[:,T]
end

# In[]
# calculate the mean and the s.t.d. of the effective input for different values of c
m = 100 # number of different values of c
T = 100000 # total simulation time for each c
μ = zeros(m) # mean of the effective input
σ = zeros(m) # s.t.d. of the effective input
IC = r*rand(n+1) # random initial conditions
c_range = range(0,stop=.4,length=m)
p = Progress(m, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)
for i = 1:m
    c = c_range[i]
    adj = quotient_net(c)
    μ[i], σ[i], IC1 = chimera(adj,T,IC)
    next!(p)
end

# In[]
writedlm("input.txt",[μ,σ],',')

# In[]
# plot showing the mean and the s.t.d. of the effective input for different values of c
scatter(c_range,μ,marker=(:circle,3,.6,:green,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),label = "μ")
scatter!(c_range,σ,marker=(:circle,3,.6,:blue,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),label = "σ")
title!("Reduced System for Chimera States")
xaxis!(L"c",font(15))
yaxis!(L"\mu,\sigma",font(15))

# In[]
# simulate the (n+1)-oscillator reduced system to obtain a typical time series of the effective input
c = .2
adj = quotient_net(c)
T = 1000000
x = zeros(n+1,T) # oscillator states
driving = zeros(T) # effective input
x[:,1] = rand(n+1,1) # random initial condition
p = Progress(T, barglyphs=BarGlyphs("[=> ]"), barlen=50, color=:yellow)
for t in 2:T
    for i in 1:n
        x[i,t] = r*f(x[i,t-1]) + K*adj[i,:]⋅(f.(x[:,t-1]).-f(x[i,t-1]))
    end
    driving[t] = K*adj[n+1,:]⋅(f.(x[:,t-1]).-f(x[n+1,t-1]))
    x[n+1,t] = r*f(x[n+1,t-1]) + driving[t]
    next!(p)
end
writedlm("driving.txt",driving,',') # time series of the effective input
@show mean(driving[T-100000:T]), std(driving[T-100000:T])

# In[]
# plot showing states of the (n+1) oscillators in the reduced system (coherent cluster: blue; incoherent cluster: green, red, magenta)
scatter(T-199:T,x[n+1,T-199:T],marker=(:circle,3,.6,:blue,stroke(1,.2,:white)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
for i = 1:n
    scatter!(T-199:T,x[i,T-199:T],marker=(:circle,3,.6,stroke(0,.2,:black)))
end
title!("Reduced System for Chimera States")
xaxis!(L"t",font(15))
yaxis!(L"x_i",font(15))

# In[]
# plot showing a typical time series of the effective input
scatter(T-1000:T,driving[T-1000:T],marker=(:circle,3,.6,:red,stroke(0,.2,:black)),background_color=RGB(0.2, 0.2, 0.2),legend=false)
title!("Reduced System for Chimera States")
xaxis!(L"t",font(15))
yaxis!("Input",font(15))


text:
using Game2048

export GameTree, get_leaf_nodes, tree_score, find_best_leaf_node, expand_node!
export back_propagate_score!

const C = 1000 # exploration constant

mutable struct GameTree
    board::Game2048.Bitboard
    n::Int
    tot_score::Float64
    pct::Rational{Int}
    parent::Union{Nothing, GameTree}
    children::Vector{GameTree}
    GameTree(board) = new(board, 0, 0, 1//1, nothing, [])
    GameTree(board, parent, pct) = new(board, 0, 0, pct, parent, [])
end

function get_leaf_nodes(gt::GameTree)
    gtarr = GameTree[]
    for child in gt.children
        get_leaf_nodes!(gtarr, child)
    end
    gtarr
end

function get_leaf_nodes!(gtarr::Vector{GameTree}, gt::GameTree)
    if length(gt.children) == 0
        push!(gtarr, gt)
    else
        for child in gt.children
            get_leaf_nodes!(gtarr, child)
        end
    end
    gtarr
end

function tree_score(gt::GameTree, N)
    gt.tot_score/gt.n + C*sqrt(log(N)/gt.n)
end

function find_best_leaf_node(leaf_nodes::Vector{GameTree}, N)
    max_score = 0
    best_leaf_node = leaf_nodes[1]
    for leaf_node in leaf_nodes
        if leaf_node.n == 0
            return leaf_node
        else
            ts = tree_score(leaf_node, N)
            if max_score < ts
                max_score = ts
                best_leaf_node = leaf_node
            end
        end
    end
    best_leaf_node
end

function expand_node!(node)
    # expand the node by populating the children
    for dir in DIRS
        new_board = move(node.board, dir)
        if new_board != node.board
            # count the number of empty spots
            zero_count = count0(new_board)
            for i in 1:zero_count
                new_board_with_new_tile = add_tile(new_board, i, 1)
                push!(node.children, GameTree(new_board_with_new_tile, node, 9//10zero_count))

                new_board_with_new_tile = add_tile(new_board, i, 2)
                push!(node.children, GameTree(new_board_with_new_tile, node, 1//10zero_count))
            end

        end
    end

    node.children
end

function back_propagate_score!(gt::GameTree, score)
    gt.n += 1
    gt.tot_score += score

    if gt.parent !== nothing
        back_propagate_score!(gt.parent, score)
    end
end

function mcts(board)
    gt = GameTree(board)

    # expand the tree node once
    @time expand_node!(gt);

    for i in 1:100
        # go through all leaf nodes and compute score
        leaf_nodes = get_leaf_nodes(gt)

        best_leaf_node = find_best_leaf_node(leaf_nodes, i)

        if best_leaf_node.n == 0
            # not been visited
            # random play to get a score
            # or run a value NN to determine score
            fnl_board = simulate_bb(best_leaf_node.board |> add_tile)
            score = sum(2 .^ bitboard_to_array(fnl_board))

            # back propagate the score
            @time back_propagate_score!(best_leaf_node, score)
        else
            # expand it as it has been visited before
            expand_node!(best_leaf_node)
        end
    end

end


text:
module CxxrtlCosim

# prologue.jl
using Yosys_jll
using Clang_jll
using Libdl
libcxxrtl = missing;

inc = joinpath(Yosys_jll.find_artifact_dir(), "share/yosys/include")
capi = joinpath(inc, "backends/cxxrtl/cxxrtl_capi.cc")
# end

mutable struct _cxxrtl_toplevel end

const cxxrtl_toplevel = Ptr{_cxxrtl_toplevel}

mutable struct _cxxrtl_handle end

const cxxrtl_handle = Ptr{_cxxrtl_handle}

function cxxrtl_create(design)
    ccall((:cxxrtl_create, libcxxrtl), cxxrtl_handle, (cxxrtl_toplevel,), design)
end

function cxxrtl_create_at(design, root)
    ccall((:cxxrtl_create_at, libcxxrtl), cxxrtl_handle, (cxxrtl_toplevel, Ptr{Cchar}), design, root)
end

function cxxrtl_destroy(handle)
    ccall((:cxxrtl_destroy, libcxxrtl), Cvoid, (cxxrtl_handle,), handle)
end

function cxxrtl_reset(handle)
    ccall((:cxxrtl_reset, libcxxrtl), Cvoid, (cxxrtl_handle,), handle)
end

function cxxrtl_eval(handle)
    ccall((:cxxrtl_eval, libcxxrtl), Cint, (cxxrtl_handle,), handle)
end

function cxxrtl_commit(handle)
    ccall((:cxxrtl_commit, libcxxrtl), Cint, (cxxrtl_handle,), handle)
end

function cxxrtl_step(handle)
    ccall((:cxxrtl_step, libcxxrtl), Csize_t, (cxxrtl_handle,), handle)
end

@enum cxxrtl_type::UInt32 begin
    CXXRTL_VALUE = 0
    CXXRTL_WIRE = 1
    CXXRTL_MEMORY = 2
    CXXRTL_ALIAS = 3
    CXXRTL_OUTLINE = 4
end

@enum cxxrtl_flag::UInt32 begin
    CXXRTL_INPUT = 1
    CXXRTL_OUTPUT = 2
    CXXRTL_INOUT = 3
    CXXRTL_DRIVEN_SYNC = 4
    CXXRTL_DRIVEN_COMB = 8
    CXXRTL_UNDRIVEN = 16
end

mutable struct _cxxrtl_outline end

struct cxxrtl_object
    type::UInt32
    flags::UInt32
    width::Csize_t
    lsb_at::Csize_t
    depth::Csize_t
    zero_at::Csize_t
    curr::Ptr{UInt32}
    next::Ptr{UInt32}
    outline::Ptr{_cxxrtl_outline}
end

function cxxrtl_get_parts(handle, name, parts)
    ccall((:cxxrtl_get_parts, libcxxrtl), Ptr{cxxrtl_object}, (cxxrtl_handle, Ptr{Cchar}, Ptr{Csize_t}), handle, name, parts)
end

function cxxrtl_enum(handle, data, callback)
    ccall((:cxxrtl_enum, libcxxrtl), Cvoid, (cxxrtl_handle, Ptr{Cvoid}, Ptr{Cvoid}), handle, data, callback)
end

const cxxrtl_outline = Ptr{_cxxrtl_outline}

function cxxrtl_outline_eval(outline)
    ccall((:cxxrtl_outline_eval, libcxxrtl), Cvoid, (cxxrtl_outline,), outline)
end

# epilogue.jl

function cxxrtl_get(handle, name)
    parts = Ref(Csize_t(0));
	object = cxxrtl_get_parts(handle, name, parts)
	@assert(object == C_NULL || parts[] == 1)
	if (object == C_NULL || parts[] == 1)
		return object
    end
	return C_NULL
end

function cxxrtl_design_create(lib)
    global libcxxrtl = lib
    ccall((:cxxrtl_design_create, libcxxrtl), cxxrtl_toplevel, ())
end

function synth(output, files...)
    cmd = Yosys_jll.yosys()
    push!(cmd.exec, "-q", "-p", "write_cxxrtl $(output)", files...)
    run(cmd)
end

function compile(output, mod)
    cmd = Clang_jll.clang()
    push!(cmd.exec, "-g", "-O3", "-fPIC", "-shared", "-std=c++14",
                    "-I$(inc)", capi, mod, "-o", output)
    run(cmd)
end

include("cosim.jl")
# end

# exports
const PREFIXES = ["cxxrtl_"]
for name in names(@__MODULE__; all=true), prefix in PREFIXES
    if startswith(string(name), prefix)
        @eval export $name
    end
end

end # module

max_stars_repo_name: JuliaComputing/CxxrtlCosim

max_issues_repo_name: JuliaComputing/CxxrtlCosim

max_forks_repo_name: JuliaComputing/CxxrtlCosim


text:
using ITensorGPU, ITensors
using Test

# Wrap an ITensor with pairs of primed and
# unprimed indices to pass to davidson
struct ITensorMap
  A::ITensor
end
Base.eltype(M::ITensorMap) = eltype(M.A)
Base.size(M::ITensorMap) = dim(IndexSet(inds(M.A; plev=0)...))
(M::ITensorMap)(v::ITensor) = noprime(M.A * v)

@testset "Complex davidson" begin
  d = 10
  i = Index(d, "i")
  A = cuITensor(randomITensor(ComplexF64, i, prime(i)))
  A = mapprime(A * mapprime(dag(A), 0, 2), 2, 1)
  M = ITensorMap(A)

  v = randomCuITensor(i)
  λ, v = davidson(M, v; maxiter=10)
  @test M(v) ≈ λ * v

  v = cuITensor(randomITensor(ComplexF64, i))
  λ, v = davidson(M, v; maxiter=10)
  @test M(v) ≈ λ * v
end