text:
# Test indexes
@testset "Indexes" begin
    # Abstract types
    @test AbstractInfOptIndex isa DataType
    @test ObjectIndex isa DataType
    @test ObjectIndex <: AbstractInfOptIndex
    # IndependentParameterIndex
    @test IndependentParameterIndex <: ObjectIndex
    @test IndependentParameterIndex(1).value == 1
    # DependentParametersIndex
    @test DependentParametersIndex <: ObjectIndex
    @test DependentParametersIndex(1).value == 1
    # DependentParameterIndex
    @test DependentParameterIndex <: AbstractInfOptIndex
    idx = DependentParametersIndex(1)
    @test DependentParameterIndex(idx, 1).object_index == idx
    @test DependentParameterIndex(idx, 1).param_index == 1
    # InfiniteParameterIndex
    @test InfiniteParameterIndex <: AbstractInfOptIndex
    # FiniteParameterIndex
    @test FiniteParameterIndex <: ObjectIndex
    @test FiniteParameterIndex(1).value == 1
    # ParameterFunctionIndex
    @test ParameterFunctionIndex <: ObjectIndex
    @test ParameterFunctionIndex(1).value == 1
    # InfiniteVariableIndex
    @test InfiniteVariableIndex <: ObjectIndex
    @test InfiniteVariableIndex(1).value == 1
    # SemiInfiniteVariableIndex
    @test SemiInfiniteVariableIndex <: ObjectIndex
    @test SemiInfiniteVariableIndex(1).value == 1
    # PointVariableIndex
    @test PointVariableIndex <: ObjectIndex
    @test PointVariableIndex(1).value == 1
    # FiniteVariableIndex
    @test FiniteVariableIndex <: ObjectIndex
    @test FiniteVariableIndex(1).value == 1
    # DerivativeIndex
    @test DerivativeIndex <: ObjectIndex
    @test DerivativeIndex(1).value == 1
    # MeasureIndex
    @test MeasureIndex <: ObjectIndex
    @test MeasureIndex(1).value == 1
    # InfOptConstraintIndex
    @test InfOptConstraintIndex <: ObjectIndex
    @test InfOptConstraintIndex(1).value == 1
    # test CleverDict extensions
    @testset "CleverDict Extensions" begin
        # index_to_key
        @test MOIUC.index_to_key(IndependentParameterIndex, Int64(42)) == IndependentParameterIndex(42)
        @test MOIUC.index_to_key(DependentParametersIndex, Int64(42)) == DependentParametersIndex(42)
        @test MOIUC.index_to_key(FiniteParameterIndex, Int64(42)) == FiniteParameterIndex(42)
        @test MOIUC.index_to_key(ParameterFunctionIndex, Int64(42)) == ParameterFunctionIndex(42)
        @test MOIUC.index_to_key(InfiniteVariableIndex, Int64(42)) == InfiniteVariableIndex(42)
        @test MOIUC.index_to_key(SemiInfiniteVariableIndex, Int64(42)) == SemiInfiniteVariableIndex(42)
        @test MOIUC.index_to_key(PointVariableIndex, Int64(42)) == PointVariableIndex(42)
        @test MOIUC.index_to_key(FiniteVariableIndex, Int64(42)) == FiniteVariableIndex(42)
        @test MOIUC.index_to_key(DerivativeIndex, Int64(42)) == DerivativeIndex(42)
        @test MOIUC.index_to_key(MeasureIndex, Int64(42)) == MeasureIndex(42)
        @test MOIUC.index_to_key(InfOptConstraintIndex, Int64(42)) == InfOptConstraintIndex(42)
        # index_to_key
        @test MOIUC.key_to_index(MeasureIndex(42)) == 42
    end
    # test Base.length
    @testset "Base.length" begin
        @test length(MeasureIndex(1)) == 1
    end
    # test Base.broadcastable
    @testset "Base.broadcastable" begin
        @test Base.broadcastable(MeasureIndex(1)) isa Ref
    end
    # test Base.iterate
    @testset "Base.iterate" begin
        @test iterate(MeasureIndex(1)) == (MeasureIndex(1), true)
        @test iterate(MeasureIndex(1), true) isa Nothing
    end
    # test Base.hash
    @testset "Base.hash" begin
        @test Base.hash(MeasureIndex(1), UInt(1)) isa UInt
    end
end

# Test Infinite Domains
@testset "Infinite Domains" begin
    # Abstract types
    @test AbstractInfiniteDomain isa DataType
    @test InfiniteScalarDomain <: AbstractInfiniteDomain
    @test InfiniteArrayDomain <: AbstractInfiniteDomain
    # IntervalDomain
    @test IntervalDomain <: InfiniteScalarDomain
    @test IntervalDomain(0, 1) == IntervalDomain(0.0, 1.0)
    @test_throws ErrorException IntervalDomain(1, 0)
    # UniDistributionDomain
    @test UniDistributionDomain <: InfiniteScalarDomain
    @test UniDistributionDomain(Normal()) isa UniDistributionDomain{<:Normal}
    # NonUnivariateDistribution
    @test InfiniteOpt.NonUnivariateDistribution <: Distribution
    # MultiDistributionDomain
    dist = Dirichlet([0.5, 0.5])
    @test MultiDistributionDomain <: InfiniteArrayDomain
    @test MultiDistributionDomain(dist) isa MultiDistributionDomain{<:Dirichlet}
    @test MultiDistributionDomain(dist) == MultiDistributionDomain(dist)
    @test MultiDistributionDomain(dist) != MultiDistributionDomain(MvNormal([1, 1], [1 0; 0 1]))
    # CollectionDomain
    @test CollectionDomain <: InfiniteArrayDomain
    @test CollectionDomain([IntervalDomain(0, 1), IntervalDomain(2, 3)]) isa CollectionDomain{IntervalDomain}
end

# Test Generative Support Info
@testset "Generative Support Info" begin
    # Abstract types
    @test AbstractGenerativeInfo isa DataType
    # NoGenerativeSupports
    @test NoGenerativeSupports <: AbstractGenerativeInfo
    @test NoGenerativeSupports() isa NoGenerativeSupports
    @test NoGenerativeSupports() == NoGenerativeSupports()
    # UniformGenerativeInfo
    @test UniformGenerativeInfo <: AbstractGenerativeInfo
    @test UniformGenerativeInfo([0.0, 1.0], InternalLabel) isa UniformGenerativeInfo
    @test UniformGenerativeInfo([-2, 0, 1], InternalLabel, -2, 1).support_basis == [0, 2.0/3.0, 1]
    @test UniformGenerativeInfo([0.0, 1.0], InternalLabel) == UniformGenerativeInfo([2, 6], InternalLabel, 2, 6)
    @test UniformGenerativeInfo([0.0, 1.0], InternalLabel) != UniformGenerativeInfo([2, 5], InternalLabel, 2, 6)
    @test_throws ErrorException UniformGenerativeInfo([2, 6], InternalLabel)
end

# Test Derivative Evaluation Methods
@testset "Derivative Methods" begin
    # abstract types
    @test AbstractDerivativeMethod isa DataType
    @test GenerativeDerivativeMethod <: AbstractDerivativeMethod
    @test NonGenerativeDerivativeMethod <: AbstractDerivativeMethod
    # test orthogonal collocation (defined in derivative_evaluations.jl)
    @test OrthogonalCollocation <: GenerativeDerivativeMethod
    @test OrthogonalCollocation{GaussLobatto}(3, GaussLobatto()) isa OrthogonalCollocation
    @test OrthogonalCollocation(3, GaussLobatto()) == OrthogonalCollocation{GaussLobatto}(3, GaussLobatto())
    @test OrthogonalCollocation(3) == OrthogonalCollocation(3, GaussLobatto())
    @test_throws ErrorException OrthogonalCollocation(1)
    # test finite difference
    @test FDTechnique isa DataType
    @test Forward <: FDTechnique
    @test Central <: FDTechnique
    @test Backward <: FDTechnique
    @test FiniteDifference <: NonGenerativeDerivativeMethod
    @test FiniteDifference(Forward()) isa FiniteDifference
    @test FiniteDifference() == FiniteDifference(Backward())
    @test FiniteDifference(Forward(), false) isa FiniteDifference
end

# Test parameter datatypes
@testset "Parameters" begin
    # abstract types
    @test InfOptParameter <: AbstractVariable
    @test ScalarParameter <: InfOptParameter
    @test AbstractDataObject isa DataType
    # test IndependentParameter
    @test IndependentParameter <: ScalarParameter
    dict = SortedDict{Float64, Set{DataType}}(2 => Set([All]))
    method = FiniteDifference()
    info = NoGenerativeSupports()
    @test IndependentParameter(IntervalDomain(0, 1), dict, 6, method, info).domain isa IntervalDomain
    # test FiniteParameter
    @test FiniteParameter <: ScalarParameter
    @test FiniteParameter(1) == FiniteParameter(1)
    # test DependentParameters
    @test DependentParameters <: InfOptParameter
    @test DependentParameters(CollectionDomain([IntervalDomain(0, 1)]),
                              Dict(zeros(1) => Set([All])), 6, [method]).domain isa CollectionDomain
    # test ScalarParameterData
    @test ScalarParameterData <: AbstractDataObject
    @test ScalarParameterData(FiniteParameter(42), 1, 1, "bob").name == "bob"
    # test MultiParameterData
    @test MultiParameterData <: AbstractDataObject
    params = DependentParameters(CollectionDomain([IntervalDomain(0, 1)]),
                                 Dict(zeros(1) => Set([All])), 6, [method])
    @test MultiParameterData(params, 1, 1:1, ["par[1]"]) isa MultiParameterData
end

# Test the InfiniteModel datatype
@testset "InfiniteModel" begin
    # test basic
    @test InfiniteModel <: JuMP.AbstractModel
    @test InfiniteModel().ext isa Dict
    # prepare optimizer constructor
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    mockattributes = MOI.OptimizerWithAttributes(mockoptimizer, MOI.Silent() => true)
    # test optimizer constructors
    @test InfiniteModel(mockoptimizer).optimizer_constructor == mockoptimizer
    @test InfiniteModel(mockattributes).optimizer_constructor == mockoptimizer
    m = InfiniteModel();
    @test isa(Base.broadcastable(m), Base.RefValue{InfiniteModel})
    @test length(JuMP.object_dictionary(m)) == 0
    @test InfiniteModel(caching_mode = MOIU.MANUAL) isa JuMP.AbstractModel
    @test InfiniteModel(mockoptimizer,
                        caching_mode = MOIU.MANUAL) isa JuMP.AbstractModel
    # test accessors
    @test InfiniteOpt._last_param_num(m) == 0
    @test InfiniteOpt._param_object_indices(m) isa Vector{Union{IndependentParameterIndex, DependentParametersIndex}}
    # test other methods
    @test empty!(InfiniteModel(mockoptimizer)).optimizer_constructor == mockoptimizer
end

# Test reference variable datatypes
@testset "References" begin
    m = InfiniteModel();
    # Abstract types
    @test DispatchVariableRef <: AbstractVariableRef
    @test FiniteRef <: DispatchVariableRef
    # test GeneralVariableRef
    @test GeneralVariableRef <: AbstractVariableRef
    @test GeneralVariableRef(m, 1, MeasureIndex) isa GeneralVariableRef
    # test IndependentParameterRef
    @test IndependentParameterRef <: DispatchVariableRef
    idx = IndependentParameterIndex(1)
    @test IndependentParameterRef(m, idx) isa IndependentParameterRef
    # test DependentParameterRef
    @test DependentParameterRef <: DispatchVariableRef
    idx = DependentParameterIndex(DependentParametersIndex(1), 1)
    @test DependentParameterRef(m, idx) isa DependentParameterRef
    # test FiniteParameterRef
    @test FiniteParameterRef <: FiniteRef
    idx = FiniteParameterIndex(1)
    @test FiniteParameterRef(m, idx) isa FiniteParameterRef
    # InfiniteVariableRef
    @test InfiniteVariableRef <: DispatchVariableRef
    idx = InfiniteVariableIndex(1)
    @test InfiniteVariableRef(m, idx) isa InfiniteVariableRef
    # SemiInfiniteVariableRef
    @test SemiInfiniteVariableRef <: DispatchVariableRef
    idx = SemiInfiniteVariableIndex(1)
    @test SemiInfiniteVariableRef(m, idx) isa SemiInfiniteVariableRef
    # ParameterFunctionRef
    @test ParameterFunctionRef <: DispatchVariableRef
    idx = ParameterFunctionIndex(1)
    @test ParameterFunctionRef(m, idx) isa ParameterFunctionRef
    # PointVariableRef
    @test PointVariableRef <: FiniteRef
    idx = PointVariableIndex(1)
    @test PointVariableRef(m, idx) isa PointVariableRef
    # FiniteVariableRef
    @test FiniteVariableRef <: FiniteRef
    idx = FiniteVariableIndex(1)
    @test FiniteVariableRef(m, idx) isa FiniteVariableRef
    # DerivativeRef
    @test DerivativeRef <: DispatchVariableRef
    idx = DerivativeIndex(1)
    @test DerivativeRef(m, idx) isa DerivativeRef
    # MeasureRef
    @test MeasureRef <: DispatchVariableRef
    idx = MeasureIndex(1)
    @test MeasureRef(m, idx) isa MeasureRef
    # InfOptConstraintRef
    @test InfOptConstraintRef isa DataType
    idx = InfOptConstraintIndex(1)
    @test InfOptConstraintRef(m, idx) isa InfOptConstraintRef
end

# Test DomainRestrictions
@testset "Domain Restrictions" begin
    # prepare the data
    m = InfiniteModel()
    idx = DependentParametersIndex(1)
    idx1 = DependentParameterIndex(idx, 1)
    idx2 = DependentParameterIndex(idx, 2)
    idx3 = IndependentParameterIndex(1)
    par1 = GeneralVariableRef(m, 1, DependentParameterIndex, 1)
    par2 = GeneralVariableRef(m, 1, DependentParameterIndex, 2)
    par3 = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
    pars = [par1, par2]
    # test datatype
    @testset "DataType" begin
        @test DomainRestrictions isa UnionAll
        d = Dict(par3 => IntervalDomain(0, 1))
        @test DomainRestrictions(d) isa DomainRestrictions{GeneralVariableRef}
        @test DomainRestrictions() isa DomainRestrictions{GeneralVariableRef}
    end
    # test _expand_parameter_tuple
    @testset "_expand_parameter_tuple" begin
        # test good to go
        d = (par3 => IntervalDomain(0, 1),)
        @test InfiniteOpt._expand_parameter_tuple(d) == Dict(d...)
        # test with array of parameters
        d = (pars => IntervalDomain(0, 1), par3 => IntervalDomain(0, 1))
        @test InfiniteOpt._expand_parameter_tuple(d) isa Dict
        # test symbolic domains
        d = (pars => [0, 1], par3 => 0)
        @test InfiniteOpt._expand_parameter_tuple(d) isa Dict
        # test bad vector domain
        d = (pars => [0, 1, 3],)
        @test_throws ErrorException InfiniteOpt._expand_parameter_tuple(d)
        # test fallback
        d = (pars => :a,)
        @test_throws ErrorException InfiniteOpt._expand_parameter_tuple(d)
    end
    # test expansion definition
    @testset "DomainRestrictions Expansion" begin
        d = (par3 => IntervalDomain(0, 1),)
        @test DomainRestrictions(d).intervals == Dict(d...)
        d = (pars => IntervalDomain(0, 1), par3 => IntervalDomain(0, 1))
        @test DomainRestrictions(d).intervals isa Dict{GeneralVariableRef, IntervalDomain}
        @test DomainRestrictions(par3 => 0).intervals isa Dict
    end
    dr = DomainRestrictions(par3 => [0, 1])
    # test intervals
    @testset "intervals" begin
        @test intervals(dr) == dr.intervals
    end
    # test Base.length
    @testset "Base.length" begin
        @test length(dr) == 1
    end
    # test Base.isempty
    @testset "Base.isempty" begin
        @test !isempty(dr)
        @test isempty(DomainRestrictions())
    end
    # test Base.:(==)
    @testset "Base.:(==)" begin
        @test dr == DomainRestrictions((par3 => IntervalDomain(0, 1),))
        @test dr != DomainRestrictions()
    end
    # test Base.copy
    @testset "Base.copy" begin
        @test copy(dr) == dr
    end
    # test Base.getindex
    @testset "Base.getindex" begin
        @test dr[par3] == IntervalDomain(0, 1)
    end
    # test Base.setindex!
    @testset "Base.setindex!" begin
        @test (dr[par3] = IntervalDomain(0, 2)) == IntervalDomain(0, 2)
    end
    # test Base.haskey
    @testset "Base.haskey" begin
        @test haskey(dr, par3)
        @test !haskey(dr, par2)
    end
    # test Base.keys
    @testset "Base.keys" begin
        @test keys(dr) == keys(intervals(dr))
    end
    # test Base.iterate
    @testset "Base.iterate" begin
        @test [r for r in dr] == [par3 => IntervalDomain(0, 2)]
    end
    # test Base.merge
    @testset "Base.merge" begin
        new_dr = DomainRestrictions()
        @test merge(new_dr, dr) == dr
    end
    # test Base.merge!
    @testset "Base.merge!" begin
        new_dr = DomainRestrictions()
        @test merge!(new_dr, dr) isa DomainRestrictions
        @test new_dr == dr
    end
    # test Base.filter
    @testset "Base.filter" begin
        @test filter(e -> e[1] != par3, dr) == DomainRestrictions()
    end
    # test Base.delete!
    @testset "Base.delete!" begin
        @test delete!(dr, par3) isa DomainRestrictions
        @test length(dr) == 0
    end
end

# Test parameter functions
@testset "Parameter Functions" begin
    # useful data
    m = InfiniteModel()
    pref = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
    vt = IC.VectorTuple(pref)
    # test ParameterFunction
    @test ParameterFunction(sin, vt, [1], [1]) isa ParameterFunction
    # test ParameterFunctionData
    @test ParameterFunctionData(ParameterFunction(sin, vt, [1], [1])) isa ParameterFunctionData
end

# Test variable datatypes
@testset "Variables" begin
    # useful data
    m = InfiniteModel()
    num = Float64(0)
    sample_info = VariableInfo(true, num, true, num, true, num, true, num, true, true)
    func = (x) -> NaN
    inf_info = VariableInfo(true, num, true, num, true, num, false, func, true, true)
    pref = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
    vref = GeneralVariableRef(m, 1, InfiniteVariableIndex)
    # Infinite variable
    @test InfiniteVariable <: JuMP.AbstractVariable
    @test InfiniteVariable(inf_info, IC.VectorTuple(pref), [1], [1], true) isa InfiniteVariable
    # Semi-Infinite variable
    @test SemiInfiniteVariable <: JuMP.AbstractVariable
    @test SemiInfiniteVariable(vref, Dict(1 => Float64(2)), [1], [1]) isa SemiInfiniteVariable
    # Point variable
    @test PointVariable <: JuMP.AbstractVariable
    @test PointVariable(sample_info, vref, Float64[1]) isa PointVariable
    # VariableData
    @test VariableData <: AbstractDataObject
    @test VariableData(ScalarVariable(sample_info)) isa VariableData
end

# Test derivative datatypes
@testset "Derivatives" begin
    # useful data
    m = InfiniteModel()
    num = Float64(0)
    sample_info = VariableInfo(true, num, true, num, true, num, true, num, true, true)
    func = (x) -> NaN
    inf_info = VariableInfo(true, num, true, num, true, num, false, func, true, true)
    pref = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
    vref = GeneralVariableRef(m, 1, InfiniteVariableIndex)
    dref = GeneralVariableRef(m, 1, DerivativeIndex)
    # derivative
    @test Derivative <: JuMP.AbstractVariable
    @test Derivative(inf_info, true, vref, pref) isa Derivative
    # Semi-infinite derivative
    @test SemiInfiniteVariable(dref, Dict(1 => Float64(2)), [1], [1]) isa SemiInfiniteVariable
    # Point derivative
    @test PointVariable(sample_info, dref, Float64[1]) isa PointVariable
    # VariableData
    @test VariableData(Derivative(inf_info, true, vref, pref)) isa VariableData
end

# Test the measure datatypes
@testset "Measures" begin
    # Helper data
    m = InfiniteModel()
    pref = GeneralVariableRef(m, 2, IndependentParameterIndex)
    w(t) = 1
    info = NoGenerativeSupports()
    # abstract types
    @test AbstractMeasureData isa DataType
    # DiscreteMeasureData
    @test DiscreteMeasureData <: AbstractMeasureData
    @test DiscreteMeasureData(pref, ones(2), ones(2), All, w, NaN, NaN, false) isa DiscreteMeasureData
    @test DiscreteMeasureData([pref], ones(2), ones(1, 2), All, w, [NaN], [NaN], false) isa DiscreteMeasureData
    # FunctionalDistributionDomain
    @test FunctionalDiscreteMeasureData <: AbstractMeasureData
    @test FunctionalDiscreteMeasureData(pref, w, 2, All, info, w, NaN, NaN, false) isa FunctionalDiscreteMeasureData
    @test FunctionalDiscreteMeasureData([pref], w, 2, All, w, [NaN], [NaN], false) isa FunctionalDiscreteMeasureData
    @test FunctionalDiscreteMeasureData([pref], w, 2, All, info, w, [NaN], [NaN], false) isa FunctionalDiscreteMeasureData
    # Measure
    @test Measure isa UnionAll
    @test Measure(zero(AffExpr),
                  DiscreteMeasureData(pref, ones(2), ones(2), All, w, NaN, NaN, false),
                  [1], [1], false) isa Measure
    # MeasureData
    @test MeasureData <: AbstractDataObject
    @test MeasureData(Measure(zero(AffExpr), DiscreteMeasureData(pref,
                      ones(2), ones(2), All, w, NaN, NaN, false), [1], [1], true)) isa MeasureData
end

# Test the constraint datatypes
@testset "Constraints" begin
    # Setup
    m = InfiniteModel()
    pref = GeneralVariableRef(m, 2, IndependentParameterIndex, -1)
    rs = DomainRestrictions(pref => [0, 1])
    con = ScalarConstraint(zero(AffExpr), MOI.Integer())
    # DomainRestrictedConstraint
    @test DomainRestrictedConstraint <: JuMP.AbstractConstraint
    @test DomainRestrictedConstraint(con, rs).restrictions == rs
    # ConstraintData
    @test ConstraintData <: AbstractDataObject
    @test ConstraintData(con, [1], "", MeasureIndex[], false) isa ConstraintData
end


text:
# Copyright 2015, 2016, 2017, 2018, 2019, 2020, 2021 Martin Holters
# See accompanying license file.

module ACME

export DiscreteModel, run!, steadystate, steadystate!, linearize, ModelRunner

using Compat: evalpoly
using IterTools: subsets
using LinearAlgebra: BLAS, I, axpy!, lu, rmul!
using Markdown: @doc_str
using OrderedCollections: OrderedDict
using ProgressMeter: @showprogress
using SparseArrays: SparseMatrixCSC, blockdiag, dropzeros!, findnz,
    nonzeros, sparse, spzeros
using StaticArrays: @SMatrix, @SVector, SMatrix, SVector
using Statistics: mean, var

include("kdtree.jl")
include("solvers.jl")

let mat_dims = Dict(
        :mv => (:nl, :nb), :mi => (:nl, :nb), :mx => (:nl, :nx),
        :mxd => (:nl, :nx), :mq => (:nl, :nq), :mu => (:nl, :nu),
        :u0 => (:nl, :n0),
        :pv => (:ny, :nb), :pi => (:ny, :nb), :px => (:ny, :nx),
        :pxd => (:ny, :nx), :pq => (:ny, :nq),
    )

    @eval function prepare_element_matrices(; $((Expr(:kw, m, nothing) for m ∈ keys(mat_dims))...))
        matrices = Dict{Symbol,SparseMatrixCSC{Real,Int}}()
        $(
            (
                quote
                    if $(mat_name) !== nothing
                        matrices[$(QuoteNode(mat_name))] = hcat($(mat_name))::AbstractMatrix
                    end
                end
                for mat_name ∈ keys(mat_dims)
            )...
        )
        sizes = Dict{Symbol,Int}(:n0 => 1)
        for (mat_name, mat) ∈ matrices
            for (sym, s) ∈ zip($(mat_dims)[mat_name], size(mat))
                if get!(sizes, sym, s) ≠ s
                    throw(ArgumentError("Inconsistent sizes for $(sym)"))
                end
            end
        end
        for (m, ns) ∈ $(mat_dims)
            if !haskey(matrices, m)
                matrices[m] = spzeros(Real, get!(sizes, ns[1], 0), get!(sizes, ns[2], 0))
            end
        end
        return matrices, sizes
    end
end

struct Element
  mv :: SparseMatrixCSC{Real,Int}
  mi :: SparseMatrixCSC{Real,Int}
  mx :: SparseMatrixCSC{Real,Int}
  mxd :: SparseMatrixCSC{Real,Int}
  mq :: SparseMatrixCSC{Real,Int}
  mu :: SparseMatrixCSC{Real,Int}
  u0 :: SparseMatrixCSC{Real,Int}
  pv :: SparseMatrixCSC{Real,Int}
  pi :: SparseMatrixCSC{Real,Int}
  px :: SparseMatrixCSC{Real,Int}
  pxd :: SparseMatrixCSC{Real,Int}
  pq :: SparseMatrixCSC{Real,Int}
  nonlinear_eq
  pins :: Dict{Symbol, Vector{Tuple{Int, Int}}}

    function Element(;nonlinear_eq=nothing, ports=nothing, pins=nothing, mat_args...)
        matrices, sizes = prepare_element_matrices(; mat_args...)
        if nonlinear_eq === nothing
            nonlinear_eq = (q) -> (SVector{0,Float64}(), SMatrix{0,0,Float64}())
        end

        if ports !== nothing
            pins = Dict{Symbol,Vector{Tuple{Int, Int}}}()
            for branch in 1:length(ports)
                push!(get!(pins, Symbol(ports[branch][1]), Tuple{Int, Int}[]), (branch, 1))
                push!(get!(pins, Symbol(ports[branch][2]), Tuple{Int, Int}[]), (branch, -1))
            end
        end
        if pins === nothing
            pins = Dict(Symbol(i) => [((i+1) ÷ 2, 2(i % 2) - 1)] for i in 1:2sizes[:nb])
        end
        return new(
            matrices[:mv], matrices[:mi], matrices[:mx], matrices[:mxd],
            matrices[:mq], matrices[:mu], matrices[:u0],
            matrices[:pv], matrices[:pi], matrices[:px], matrices[:pxd], matrices[:pq],
            nonlinear_eq,
            pins,
        )
    end
end

for (n,m) in Dict(:nb => :mv, :nx => :mx, :nq => :mq, :nu => :mu)
  @eval ($n)(e::Element) = size(e.$m, 2)
end
nl(e::Element) = size(e.mv, 1)
ny(e::Element) = size(e.pv, 1)
nn(e::Element) = nb(e) + nx(e) + nq(e) - nl(e)

nonlinear_eq_func(e::Element) = e.nonlinear_eq

include("elements.jl")

include("circuit.jl")

mutable struct DiscreteModel{Solvers}
    a::Matrix{Float64}
    b::Matrix{Float64}
    c::Matrix{Float64}
    x0::Vector{Float64}
    pexps::Vector{Matrix{Float64}}
    dqs::Vector{Matrix{Float64}}
    eqs::Vector{Matrix{Float64}}
    fqprevs::Vector{Matrix{Float64}}
    fqs::Vector{Matrix{Float64}}
    q0s::Vector{Vector{Float64}}
    dy::Matrix{Float64}
    ey::Matrix{Float64}
    fy::Matrix{Float64}
    y0::Vector{Float64}

    nonlinear_eq_funcs::Vector

    solvers::Solvers
    x::Vector{Float64}

    function DiscreteModel(mats::Dict{Symbol}, nonlinear_eq_funcs::Vector,
            solvers::Solvers) where {Solvers}
        model = new{Solvers}()

        for mat in (:a, :b, :c, :pexps, :dqs, :eqs, :fqprevs, :fqs, :dy, :ey, :fy, :x0, :q0s, :y0)
            setfield!(model, mat, convert(fieldtype(typeof(model), mat), mats[mat]))
        end

        model.nonlinear_eq_funcs = nonlinear_eq_funcs
        model.solvers = solvers
        model.x = zeros(nx(model))
        return model
    end
end

function DiscreteModel(circ::Circuit, t::Real, ::Type{Solver}=HomotopySolver{CachingSolver{SimpleSolver}};
                       decompose_nonlinearity=true) where {Solver}
    mats = model_matrices(circ, t)

    nns = Int[nn(e) for e in elements(circ)]
    nqs = Int[nq(e) for e in elements(circ)]
    if decompose_nonlinearity
        nl_elems = nldecompose!(mats, nns, nqs)
    else
        nl_elems = Vector{Int}[findall(nn -> nn > 0, nns)]
    end

    model_nns = Int[sum(nns[nles]) for nles in nl_elems]
    model_qidxs = [vcat(consecranges(nqs)[nles]...) for nles in nl_elems]
    split_nl_model_matrices!(mats, model_qidxs, model_nns)

    reduce_pdims!(mats)

    model_nps = size.(mats[:dqs], 1)
    model_nqs = size.(mats[:pexps], 1)

    @assert nn(circ) == sum(model_nns)

    q0s = Vector{Float64}.(mats[:q0s])
    fqs = Matrix{Float64}.(mats[:fqs])
    fqprev_fulls = Matrix{Float64}.(mats[:fqprev_fulls])

    model_nonlinear_eq_funcs = [
        let q = zeros(nq), circ_nl_func = nonlinear_eq_func(circ, nles)
            @inline function(res, J, pfull, Jq, fq, z)
                #copyto!(q, pfull + fq * z)
                copyto!(q, pfull)
                BLAS.gemv!('N', 1., fq, z, 1., q)
                res´, Jq´ = circ_nl_func(q)
                res .= res´
                Jq .= Jq´
                #copyto!(J, Jq*model.fq)
                BLAS.gemm!('N', 'N', 1., Jq, fq, 0., J)
                return nothing
            end
        end for (nles, nq) in zip(nl_elems, model_nqs)]

    nonlinear_eq_funcs = [
        @inline function (res, J, scratch, z)
            nleq(res, J, scratch[1], scratch[2], fq, z)
        end for (nleq, fq) in zip(model_nonlinear_eq_funcs, fqs)]

    init_zs = [zeros(nn) for nn in model_nns]
    for idx in eachindex(nonlinear_eq_funcs)
        q = q0s[idx] + fqprev_fulls[idx] * vcat(init_zs...)
        init_zs[idx] = initial_solution(nonlinear_eq_funcs[idx], q, model_nns[idx])
    end

    while any(np -> np == 0, model_nps)
        const_idxs = findall(iszero, model_nps)
        const_zidxs = vcat(consecranges(model_nns)[const_idxs]...)
        varying_zidxs = filter(idx -> !(idx in const_zidxs), 1:sum(model_nns))
        for idx in eachindex(mats[:q0s])
            mats[:q0s][idx] += mats[:fqprev_fulls][idx][:,const_zidxs] * vcat(init_zs[const_idxs]...)
            mats[:fqprev_fulls][idx] = mats[:fqprev_fulls][idx][:,varying_zidxs]
        end
        mats[:x0] += mats[:c][:,const_zidxs] * vcat(init_zs[const_idxs]...)
        mats[:y0] += mats[:fy][:,const_zidxs] * vcat(init_zs[const_idxs]...)
        deleteat!(mats[:q0s], const_idxs)
        deleteat!(mats[:dq_fulls], const_idxs)
        deleteat!(mats[:eq_fulls], const_idxs)
        deleteat!(mats[:fqs], const_idxs)
        deleteat!(mats[:fqprev_fulls], const_idxs)
        deleteat!(init_zs, const_idxs)
        deleteat!(model_nns, const_idxs)
        deleteat!(model_nqs, const_idxs)
        deleteat!(model_nonlinear_eq_funcs, const_idxs)
        deleteat!(nonlinear_eq_funcs, const_idxs)
        deleteat!(nl_elems, const_idxs)
        mats[:fy] = mats[:fy][:,varying_zidxs]
        mats[:c] = mats[:c][:,varying_zidxs]
        reduce_pdims!(mats)
        model_nps = size.(mats[:dqs], 1)
    end

    q0s = Array{Float64}.(mats[:q0s])
    fqs = Array{Float64}.(mats[:fqs])
    fqprev_fulls = Array{Float64}.(mats[:fqprev_fulls])
    pexps = Array{Float64}.(mats[:pexps])

    nonlinear_eq_set_ps = [
        function(scratch, p)
            pfull = scratch[1]
            #copyto!(pfull, q0 + pexp * p)
            copyto!(pfull, q0)
            BLAS.gemv!('N', 1., pexp, p, 1., pfull)
            return nothing
        end
        for (pexp, q0) in zip(pexps, q0s)]
    nonlinear_eq_calc_Jps = [
        function (scratch, Jp)
            Jq = scratch[2]
            #copyto!(Jp, Jq*pexp)
            BLAS.gemm!('N', 'N', 1., Jq, pexp, 0., Jp)
            return nothing
        end
        for pexp in pexps]
    solvers = ((Solver(
                                  ParametricNonLinEq(nonlinear_eq_funcs[idx],
                                      nonlinear_eq_set_ps[idx],
                                      nonlinear_eq_calc_Jps[idx],
                                      (zeros(model_nqs[idx]), zeros(model_nns[idx], model_nqs[idx])),
                                      model_nns[idx], model_nps[idx]),
                                  zeros(model_nps[idx]), init_zs[idx])
                for idx in eachindex(nonlinear_eq_funcs))...,)
    return DiscreteModel(mats, model_nonlinear_eq_funcs, solvers)
end

function model_matrices(circ::Circuit, t::Rational{BigInt})
    lhs = convert(SparseMatrixCSC{Rational{BigInt},Int},
                  [mv(circ) mi(circ) mxd(circ)//t+mx(circ)//2 mq(circ);
                   blockdiag(topomat(circ)...) spzeros(nb(circ), nx(circ) + nq(circ))])
    rhs = convert(SparseMatrixCSC{Rational{BigInt},Int},
                  [u0(circ) mu(circ) mxd(circ)//t-mx(circ)//2;
                          spzeros(nb(circ), 1+nu(circ)+nx(circ))])
    x, f = Matrix.(gensolve(lhs, rhs))

    rowsizes = [nb(circ); nb(circ); nx(circ); nq(circ)]
    res = Dict{Symbol,Array}(zip([:fv; :fi; :c; :fq], matsplit(f, rowsizes)))

    nullspace = gensolve(sparse(res[:fq]::Matrix{Rational{BigInt}}),
                         spzeros(Rational{BigInt}, size(res[:fq],1), 0))[2]
    indeterminates = f * nullspace

    if sum(abs2, res[:c] * nullspace) > 1e-20
        @warn "State update depends on indeterminate quantity"
    end
    while size(nullspace, 2) > 0
        i, j = argmax(abs.(nullspace)).I
        nullspace = nullspace[[1:i-1; i+1:end], [1:j-1; j+1:end]]
        f = f[:, [1:i-1; i+1:end]]
        for k in [:fv; :fi; :c; :fq]
            res[k] = res[k][:, [1:i-1; i+1:end]]
        end
    end

    merge!(res, Dict(zip([:v0 :ev :dv; :i0 :ei :di; :x0 :b :a; :q0 :eq_full :dq_full],
                         matsplit(x, rowsizes, [1; nu(circ); nx(circ)]))))
    for v in (:v0, :i0, :x0, :q0)
        res[v] = dropdims(res[v], dims=2)
    end

    p = [pv(circ) pi(circ) px(circ)//2+pxd(circ)//t pq(circ)]
    if sum(abs2, p * indeterminates) > 1e-20
        @warn "Model output depends on indeterminate quantity"
    end
    res[:dy] = p * x[:,2+nu(circ):end] + px(circ)//2-pxd(circ)//t
    #          p * [dv; di; a;  dq_full] + px(circ)//2-pxd(circ)//t
    res[:ey] = p * x[:,2:1+nu(circ)] # p * [ev; ei; b;  eq_full]
    res[:fy] = p * f                 # p * [fv; fi; c;  fq]
    res[:y0] = p * vec(x[:,1])       # p * [v0; i0; x0; q0]

    return res
end

model_matrices(circ::Circuit, t) = model_matrices(circ, Rational{BigInt}(t))

function tryextract(fq, numcols)
    a = Matrix{eltype(fq)}(I, size(fq, 2), size(fq, 2))
    if numcols ≥ size(fq,2)
        return a
    end
    for colcnt in 1:numcols
        # determine element with maximum absolute value in unprocessed columns
        # to use as pivot
        i, j = argmax(abs.(fq[:,colcnt:end])).I
        j += colcnt-1

        # swap pivot to first (unprocessed) column
        fq[:,[colcnt, j]] = fq[:,[j, colcnt]]
        a[:,[colcnt, j]] = a[:,[j, colcnt]]

        # elimnate remaining columns in i-th row and perform equivalent
        # transformation to a
        jj = colcnt+1:size(fq,2)
        a[:,jj] -= a[:,colcnt] * (fq[[i],jj] / fq[i,colcnt])
        fq[:,jj] -= fq[:,colcnt] * (fq[[i],jj] / fq[i,colcnt])
        # ignore i-th row in following processing steps
        fq = fq[[1:i-1; i+1:end],:]

        if all(iszero, fq[:,colcnt+1:end])
            return a
        end
    end
    return nothing
end

function nldecompose!(mats, nns, nqs)
    fq = mats[:fq]
    a = Matrix{eltype(fq)}(I, size(fq, 2), size(fq, 2))
    sub_ranges = consecranges(nqs)
    extracted_subs = Vector{Int}[]
    rem_cols = 1:size(fq, 2)
    rem_nles = BitSet(filter!(e -> nqs[e] > 0, collect(eachindex(nqs))))

    while !isempty(rem_nles)
        for sz in 1:length(rem_nles), sub in subsets(collect(rem_nles), sz)
            nn_sub = sum(nns[sub])
            a_update = tryextract(fq[[sub_ranges[sub]...;],rem_cols], nn_sub)
            if a_update !== nothing
                fq[:,rem_cols] = fq[:,rem_cols] * a_update
                a[:,rem_cols] = a[:,rem_cols] * a_update
                rem_cols = first(rem_cols)+nn_sub:size(fq, 2)
                push!(extracted_subs, sub)
                for nle in sub
                    delete!(rem_nles, nle)
                end
                break
            end
        end
    end

    mats[:c] = mats[:c] * a
    # mats[:fq] is updated as part of the loop
    mats[:fy] = mats[:fy] * a
    return extracted_subs
end

function split_nl_model_matrices!(mats, model_qidxs, model_nns)
    mats[:dq_fulls] = Matrix[mats[:dq_full][qidxs,:] for qidxs in model_qidxs]
    mats[:eq_fulls] = Matrix[mats[:eq_full][qidxs,:] for qidxs in model_qidxs]
    let fqsplit = vcat((matsplit(mats[:fq][qidxs,:], [length(qidxs)], model_nns) for qidxs in model_qidxs)...)
        mats[:fqs] = Matrix[fqsplit[i,i] for i in 1:length(model_qidxs)]
        mats[:fqprev_fulls] = Matrix[[fqsplit[i, 1:i-1]... zeros(eltype(mats[:fq]), length(model_qidxs[i]), sum(model_nns[i:end]))]
                                     for i in 1:length(model_qidxs)]
    end
    mats[:q0s] = Vector[mats[:q0][qidxs] for qidxs in model_qidxs]
end

function reduce_pdims!(mats::Dict)
    subcount = length(mats[:dq_fulls])
    mats[:dqs] = Vector{Matrix}(undef, subcount)
    mats[:eqs] = Vector{Matrix}(undef, subcount)
    mats[:fqprevs] = Vector{Matrix}(undef, subcount)
    mats[:pexps] = Vector{Matrix}(undef, subcount)
    offset = 0
    for idx in 1:subcount
        # decompose [dq_full eq_full] into pexp*[dq eq] with [dq eq] having minimum
        # number of rows
        pexp, dqeq = rank_factorize(sparse([mats[:dq_fulls][idx] mats[:eq_fulls][idx] mats[:fqprev_fulls][idx]]))
        mats[:pexps][idx] = pexp
        colsizes = [size(mats[m][idx], 2) for m in [:dq_fulls, :eq_fulls, :fqprev_fulls]]
        mats[:dqs][idx], mats[:eqs][idx], mats[:fqprevs][idx] = matsplit(dqeq, [size(dqeq, 1)], colsizes)

        # project pexp onto the orthogonal complement of the column space of Fq
        fq = mats[:fqs][idx]
        nn = size(fq, 2)
        fq_pinv = gensolve(sparse(fq'*fq), fq')[1]
        pexp = pexp - fq*fq_pinv*pexp
        # if the new pexp has lower rank, update
        pexp, f = rank_factorize(sparse(pexp))
        if size(pexp, 2) < size(mats[:pexps][idx], 2)
            cols = offset .+ (1:nn)
            mats[:a] = mats[:a] - mats[:c][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:dqs][idx]
            mats[:b] = mats[:b] - mats[:c][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:eqs][idx]
            mats[:dy] = mats[:dy] - mats[:fy][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:dqs][idx]
            mats[:ey] = mats[:ey] - mats[:fy][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:eqs][idx]
            for idx2 in (idx+1):subcount
                mats[:dq_fulls][idx2] = mats[:dq_fulls][idx2] - mats[:fqprev_fulls][idx2][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:dqs][idx]
                mats[:eq_fulls][idx2] = mats[:eq_fulls][idx2] - mats[:fqprev_fulls][idx2][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:eqs][idx]
                mats[:fqprev_fulls][idx2][:,1:offset] = mats[:fqprev_fulls][idx2][:,1:offset] - mats[:fqprev_fulls][idx2][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:fqprevs][idx][:,1:offset]
            end
            mats[:pexps][idx] = pexp
            mats[:dqs][idx] = f * mats[:dqs][idx]
            mats[:eqs][idx] = f * mats[:eqs][idx]
            mats[:fqprevs][idx] = f * mats[:fqprevs][idx]
            mats[:dq_fulls][idx] = pexp * mats[:dqs][idx]
            mats[:eq_fulls][idx] = pexp * mats[:eqs][idx]
            mats[:fqprev_fulls][idx] = pexp * mats[:fqprevs][idx]
        end
        offset += nn
    end
end

function initial_solution(init_nl_eq_func::Function, q0, nn)
    # determine an initial solution with a homotopy solver that may vary q0
    # between 0 and the true q0 -> q0 takes the role of p
    nq = length(q0)
    init_nleq = ParametricNonLinEq(init_nl_eq_func, nn, nq)
    init_solver = HomotopySolver{SimpleSolver}(init_nleq, zeros(nq), zeros(nn))
    init_z = solve(init_solver, q0)
    if !hasconverged(init_solver)
        error("Failed to find initial solution")
    end
    return init_z
end

nx(model::DiscreteModel) = length(model.x0)
nq(model::DiscreteModel, subidx) = length(model.q0s[subidx])
np(model::DiscreteModel, subidx) = size(model.dqs[subidx], 1)
nu(model::DiscreteModel) = size(model.b, 2)
ny(model::DiscreteModel) = length(model.y0)
nn(model::DiscreteModel, subidx) = size(model.fqs[subidx], 2)
nn(model::DiscreteModel) = reduce(+, init=0, size(fq, 2) for fq in model.fqs)

function steadystate(model::DiscreteModel, u=zeros(nu(model)))
    IA_LU = lu(I-model.a)
    steady_z = zeros(nn(model))
    zoff = 1
    for idx in 1:length(model.solvers)
        zoff_last = zoff+nn(model,idx)-1
        steady_q0 = model.q0s[idx] + model.pexps[idx]*((model.dqs[idx]/IA_LU*model.b + model.eqs[idx])*u + (model.dqs[idx]/IA_LU*model.c + model.fqprevs[idx])*steady_z) +
            model.pexps[idx]*model.dqs[idx]/IA_LU*model.x0
        fq = model.pexps[idx]*model.dqs[idx]/IA_LU*model.c[:,zoff:zoff_last] + model.fqs[idx]
        nleq = model.nonlinear_eq_funcs[idx]
        steady_nl_eq_func =
            (res, J, scratch, z) -> nleq(res, J, scratch[1], scratch[2], fq, z)
        steady_nleq = ParametricNonLinEq(steady_nl_eq_func, nn(model, idx), nq(model, idx))
        steady_solver = HomotopySolver{SimpleSolver}(steady_nleq, zeros(nq(model, idx)),
                                                     zeros(nn(model, idx)))
        set_resabstol!(steady_solver, 1e-15)
        steady_z[zoff:zoff_last] = solve(steady_solver, steady_q0)
        if !hasconverged(steady_solver)
            error("Failed to find steady state solution")
        end
        zoff += nn(model,idx)
    end
    return IA_LU\(model.b*u + model.c*steady_z + model.x0)
end

function steadystate!(model::DiscreteModel, u=zeros(nu(model)))
    x_steady = steadystate(model, u)
    copyto!(model.x, x_steady)
    return x_steady
end

function linearize(model::DiscreteModel, usteady::AbstractVector{Float64}=zeros(nu(model)))
    xsteady = steadystate(model, usteady)
    zranges = Vector{UnitRange{Int64}}(undef, length(model.solvers))
    dzdps = Vector{Matrix{Float64}}(undef, length(model.solvers))
    dqlins = Vector{Matrix{Float64}}(undef, length(model.solvers))
    eqlins = Vector{Matrix{Float64}}(undef, length(model.solvers))
    zsteady = zeros(nn(model))
    zoff = 1
    x0 = copy(model.x0)
    a = copy(model.a)
    b = copy(model.b)
    c = copy(model.c)
    y0 = copy(model.y0)
    dy = copy(model.dy)
    ey = copy(model.ey)
    fy = copy(model.fy)

    for idx in 1:length(model.solvers)
        psteady = model.dqs[idx] * xsteady + model.eqs[idx] * usteady +
                  model.fqprevs[idx] * zsteady
        zsub, dzdps[idx] = linearize(model.solvers[idx], psteady)
        copyto!(zsteady, zoff, zsub, 1, length(zsub))

        zranges[idx] = zoff:zoff+length(zsub)-1
        fqdzdps = [model.fqprevs[idx][:,zranges[n]] * dzdps[n] for n in 1:idx-1]
        dqlins[idx] = reduce(+, init=model.dqs[idx], fqdzdps .* dqlins[1:idx-1])
        eqlins[idx] = reduce(+, init=model.eqs[idx], fqdzdps .* eqlins[1:idx-1])

        x0 += model.c[:,zranges[idx]] * (zsub - dzdps[idx]*psteady)
        a += model.c[:,zranges[idx]] * dzdps[idx] * dqlins[idx]
        b += model.c[:,zranges[idx]] * dzdps[idx] * eqlins[idx]

        y0 += model.fy[:,zranges[idx]] * (zsub - dzdps[idx]*psteady)
        dy += model.fy[:,zranges[idx]] * dzdps[idx] * dqlins[idx]
        ey += model.fy[:,zranges[idx]] * dzdps[idx] * eqlins[idx]

        zoff += length(zsub)
    end

    mats = Dict(:a => a, :b => b, :c => zeros(nx(model), 0),
        :pexps => Matrix{Float64}[], :dqs => Matrix{Float64}[],
        :eqs => Matrix{Float64}[], :fqprevs => Matrix{Float64}[],
        :fqs => Matrix{Float64}[], :q0s => Vector{Float64}[],
        :dy => dy, :ey => ey, :fy => zeros(ny(model), 0), :x0 => x0, :y0 => y0)
    return DiscreteModel(mats, [], ())
end

"""
    run!(model::DiscreteModel, u::AbstractMatrix{Float64}; showprogress=true)

Run the given `model` by feeding it the input `u` which must be a matrix with
one row for each of the circuit's inputs and one column for each time step to
simulate. Likewise, the returned output will be a matrix with one row for each
of the circuit's outputs and one column for each simulated time step. The order
of the rows will correspond to the order in which the respective input and
output elements were added to the `Circuit`. To simulate a circuit without
inputs, a matrix with zero rows may be passed. The internal state of the model
(e.g. capacitor charges) is preserved accross calls to `run!`.

By default `run!` will show a progress bar to report its progress. This can be
disabled by passing `showprogress=false`.
"""
run!(model::DiscreteModel, u::AbstractMatrix{Float64}; showprogress=true) =
    return run!(ModelRunner(model, showprogress), u)

struct ModelRunner{Model<:DiscreteModel,ShowProgress}
    model::Model
    ucur::Vector{Float64}
    ps::Vector{Vector{Float64}}
    ycur::Vector{Float64}
    xnew::Vector{Float64}
    z::Vector{Float64}
    function ModelRunner{Model,ShowProgress}(model::Model) where {Model<:DiscreteModel,ShowProgress}
        ucur = Vector{Float64}(undef, nu(model))
        ps = Vector{Float64}[Vector{Float64}(undef, np(model, idx)) for idx in 1:length(model.solvers)]
        ycur = Vector{Float64}(undef, ny(model))
        xnew = Vector{Float64}(undef, nx(model))
        z = Vector{Float64}(undef, nn(model))
        return new{Model,ShowProgress}(model, ucur, ps, ycur, xnew, z)
    end
end

ModelRunner(model::Model) where {Model<:DiscreteModel} = ModelRunner{Model,true}(model)
ModelRunner(model::Model, ::Val{ShowProgress}) where {Model<:DiscreteModel,ShowProgress} =
    ModelRunner{Model,ShowProgress}(model)

"""
    ModelRunner(model::DiscreteModel, showprogress::Bool = true)

Construct a `ModelRunner` instance for running `model`. The `ModelRunner`
instance pre-allocates some memory required for model execution. Hence, when
running the same model for multiple small input data blocks, some overhead can
be saved by explicitly using a `ModelRunner`.

By default `run!` for the constructed `ModelRunner` will show a progress bar to
report its progress. This can be disabled by passing `false` as second
parameter.
"""
ModelRunner(model::Model, showprogress::Bool) where {Model<:DiscreteModel} =
    ModelRunner{Model,showprogress}(model)

"""
    run!(runner::ModelRunner, u::AbstractMatrix{Float64})

Run the given `runner` by feeding it the input `u` which must be a matrix with
one row for each of the circuit's inputs and one column for each time step to
simulate. Likewise, the returned output will be a matrix with one row for each
of the circuit's outputs and one column for each simulated time step. The order
of the rows will correspond to the order in which the respective input and
output elements were added to the `Circuit`. To simulate a circuit without
inputs, a matrix with zero rows may be passed. The internal state of the
underlying `DiscreteModel` (e.g. capacitor charges) is preserved accross calls
to `run!`.
"""
function run!(runner::ModelRunner, u::AbstractMatrix{Float64})
    y = Matrix{Float64}(undef, ny(runner.model), size(u, 2))
    run!(runner, y, u)
    return y
end

function checkiosizes(runner::ModelRunner, u::AbstractMatrix{Float64}, y::AbstractMatrix{Float64})
    if size(u, 1) ≠ nu(runner.model)
        throw(DimensionMismatch("input matrix has $(size(u,1)) rows, but model has $(nu(runner.model)) inputs"))
    end
    if size(y, 1) ≠ ny(runner.model)
        throw(DimensionMismatch("output matrix has $(size(y,1)) rows, but model has $(ny(runner.model)) outputs"))
    end
    if size(u, 2) ≠ size(y, 2)
        throw(DimensionMismatch("input matrix has $(size(u,2)) columns, output matrix has $(size(y,2)) columns"))
    end
end

"""
    run!(runner::ModelRunner, y::AbstractMatrix{Float64}, u::AbstractMatrix{Float64})

Run the given `runner` by feeding it the input `u` and storing the output
in `y`. The input `u` must be a matrix with one row for each of the circuit's
inputs and one column for each time step to simulate. Likewise, the output `y`
must be a matrix with one row for each of the circuit's outputs and one column
for each simulated time step. The order of the rows will correspond to the order
in which the respective input and output elements were added to the `Circuit`.
To simulate a circuit without inputs, a matrix with zero rows may be passed.
The internal state of the  underlying `DiscreteModel` (e.g. capacitor charges)
is preserved accross calls to `run!`.
"""
function run!(runner::ModelRunner{<:DiscreteModel,true}, y::AbstractMatrix{Float64},
              u::AbstractMatrix{Float64})
    checkiosizes(runner, u, y)
    @showprogress "Running model: " for n = 1:size(u, 2)
        step!(runner, y, u, n)
    end
end

function run!(runner::ModelRunner{<:DiscreteModel,false}, y::AbstractMatrix{Float64},
              u::AbstractMatrix{Float64})
    checkiosizes(runner, u, y)
    for n = 1:size(u, 2)
        step!(runner, y, u, n)
    end
end

function step!(runner::ModelRunner, y::AbstractMatrix{Float64}, u::AbstractMatrix{Float64}, n)
    model = runner.model
    ucur = runner.ucur
    ycur = runner.ycur
    xnew = runner.xnew
    z = runner.z
    copyto!(ucur, 1, u, (n-1)*nu(model)+1, nu(model))
    zoff = 1
    fill!(z, 0.0)
    for idx in 1:length(model.solvers)
        p = runner.ps[idx]
        # copyto!(p, model.dqs[idx] * model.x + model.eqs[idx] * u[:,n]) + model.fqprevs[idx] * z
        if size(model.dqs[idx], 2) == 0
            fill!(p, 0.0)
        else
            BLAS.gemv!('N', 1., model.dqs[idx], model.x, 0., p)
        end
        BLAS.gemv!('N', 1., model.eqs[idx], ucur, 1., p)
        if idx > 1
            BLAS.gemv!('N', 1., model.fqprevs[idx], z, 1., p)
        end
        zsub = solve(model.solvers[idx], p)
        if !hasconverged(model.solvers[idx])
            if all(isfinite, zsub)
                @warn "Failed to converge while solving non-linear equation."
            else
                error("Failed to converge while solving non-linear equation, got non-finite result.")
            end
        end
        copyto!(z, zoff, zsub, 1, length(zsub))
        zoff += length(zsub)
    end
    #y[:,n] = model.dy * model.x + model.ey * u[:,n] + model.fy * z + model.y0
    if !isempty(ycur)
        copyto!(ycur, model.y0)
        BLAS.gemv!('N', 1., model.dy, model.x, 1., ycur)
        BLAS.gemv!('N', 1., model.ey, ucur, 1., ycur)
        BLAS.gemv!('N', 1., model.fy, z, 1., ycur)
        #y[:,n] = ycur
        copyto!(y, (n-1)*ny(model)+1, ycur, 1, ny(model))
    end
    #model.x = model.a * model.x + model.b * u[:,n] + model.c * z + model.x0
    if !isempty(xnew)
        copyto!(xnew, model.x0)
        BLAS.gemv!('N', 1., model.a, model.x, 1., xnew)
        BLAS.gemv!('N', 1., model.b, ucur, 1.,xnew)
        BLAS.gemv!('N', 1., model.c, z, 1., xnew)
        copyto!(model.x, xnew)
    end
end

function gensolve(a, b, x, h, thresh=0.1)
    m = size(a)[1]
    if m == 0
        return x, h
    end
    t = sortperm(vec(mapslices(ait -> count(!iszero, ait), a, dims=2))) # row indexes in ascending order of nnz
    tol = 3 * max(eps(float(eltype(a))), eps(float(eltype(h)))) * size(a, 2)
    for i in 1:m
        ait = a[t[i],:]' # ait is a row of the a matrix
        s = ait * h;
        jnz, nz_vals = findnz(s')
        nz_abs_vals = abs.(nz_vals)
        max_abs_val = reduce(max, init=zero(eltype(s)), nz_abs_vals)
        if max_abs_val ≤ tol # cosidered numerical zero
            continue
        end
        jat = jnz[nz_abs_vals .≥ thresh*max_abs_val] # cols above threshold
        j = jat[argmin(vec(mapslices(hj -> count(!iszero, hj), h[:,jat], dims=1)))]
        q = h[:,j]
        x = x + convert(typeof(x), q * ((b[t[i],:]' - ait*x) * (1 / (ait*q))))
        if size(h)[2] > 1
            h = h[:,[1:j-1;j+1:end]] - convert(typeof(h), q * s[1,[1:j-1;j+1:end]]'*(1/s[1,j]))
        else
            h = similar(h, eltype(h), (size(h)[1], 0))
        end
    end
    return x, h
end

gensolve(a, b, thresh=0.1) =
    gensolve(a, b, spzeros(promote_type(eltype(a), eltype(b)), size(a, 2), size(b, 2)), SparseMatrixCSC{eltype(a)}(I, size(a,2), size(a,2)), thresh)

function rank_factorize(a::SparseMatrixCSC)
    f = a
    nullspace = gensolve(a', spzeros(eltype(a), size(a, 2), 0))[2]
    c = Matrix{eltype(a)}(I, size(a, 1), size(a, 1))
    while size(nullspace, 2) > 0
        i, j = argmax(abs.(nullspace)).I
        c -= c[:, i] * nullspace[:, j]' / nullspace[i, j]
        c = c[:, [1:i-1; i+1:end]]
        nullspace -= nullspace[:, j] * vec(nullspace[i, :])' / nullspace[i, j]
        nullspace = nullspace[[1:i-1; i+1:end], [1:j-1; j+1:end]]
        f = f[[1:i-1; i+1:end], :]
    end
    return c, f
end

consecranges(lengths) = map((l, e) -> (e-l+1):e, lengths, cumsum(lengths))

matsplit(v::AbstractVector, rowsizes) = [v[rs] for rs in consecranges(rowsizes)]
matsplit(m::AbstractMatrix, rowsizes, colsizes=[size(m,2)]) =
    [m[rs, cs] for rs in consecranges(rowsizes), cs in consecranges(colsizes)]

precompile(resistor, (Float64,))
precompile(resistor, (Int,))
precompile(potentiometer, (Float64,))
precompile(potentiometer, (Int,))
precompile(potentiometer, (Float64, Float64))
precompile(potentiometer, (Int, Float64))
precompile(capacitor, (Float64,))
precompile(inductor, (Float64,))
precompile(inductor, (Type{Val{:JA}},))
precompile(transformer, (Float64, Float64))
precompile(transformer, (Type{Val{:JA}},))
precompile(voltagesource, ())
precompile(voltagesource, (Float64,))
precompile(voltagesource, (Int,))
precompile(currentsource, ())
precompile(currentsource, (Float64,))
precompile(currentsource, (Int,))
precompile(voltageprobe, ())
precompile(currentprobe, ())
precompile(diode, ())
precompile(bjt, (Symbol,))
precompile(mosfet, (Symbol,))
precompile(opamp, ())
precompile(opamp, (Type{Val{:macak}}, Float64, Float64, Float64))

end # module

ext: jl


text:
using LibHello
using Test

@testset "greet" begin
    @test greet() == "hello, world"
end

@testset "get_array" begin
    @test get_array(5) == Int32[0,1,2,3,4]
end


text:
__precompile__()
# ----------------------------------------------------------------------------
module Knapsack
#=
    see also:
        - the walkthrough at http://vladium.com/tutorials/study_julia_with_me/knapsack_benchmark/
        - https://en.wikipedia.org/wiki/Knapsack_problem#0/1_knapsack_problem

    running this benchmark from command line:

    >julia Knapsack.jl
=#
# ............................................................................
"""
    we need a couple of integers to describe a knapsack item, its value and weight.
    a tuple would work, but a struct makes for more readable code
"""
struct Item
    value   ::Int64; # actually, being so type-specific is not necessary in Julia
    weight  ::Int64;
end
# ............................................................................
"""
    first version of the knapsack solver; uses a full matrix for V and will
    not scale to large problems
"""
function opt_value_BAD(W ::Int64, items ::Array{Item}) ::Int64
    n = length(items)

    # V[w,j] stores opt value achievable with capacity 'w' and using items '1..j':

    V = Array{Int64}(undef, W, n) # Wxn matrix with uninitialized storage

    # initialize first column v[:, 1] to trivial single-item solutions:
    # (note "broadcast assignment" syntax)

    V[:, 1] .= 0
    V[items[1].weight:end, 1] .= items[1].value

    # do a pass through remaining columns:

    for j in 2 : n
        itemⱼ = items[j]
        for w in 1 : W
            V_without_itemⱼ = V[w, j - 1]
            V_allow_itemⱼ = (w < itemⱼ.weight
                ? V_without_itemⱼ
                : (itemⱼ.value + (w ≠ itemⱼ.weight ? V[w - itemⱼ.weight, j - 1] : 0)))
            V[w, j] = max(V_allow_itemⱼ, V_without_itemⱼ)
        end
    end

    return V[W, n]
end
# ............................................................................
"""
    second version of the knapsack solver; uses two column buffers for V
"""
function opt_value(W ::Int64, items ::Array{Item}) ::Int64
    n = length(items)

    # V[w] stores opt value achievable with capacity 'w' and using items '1..j':

    V = zeros(Int64, W) # single column of size W, zero-initialized
    V_prev = Array{Int64}(undef, W) # single column of size W, uninitialized storage

    # initialize 'v' to trivial single-item solutions:
    # (note "broadcast assignment" syntax)

    V[items[1].weight:end] .= items[1].value

    # do a pass through remaining columns:

    for j in 2 : n
        V, V_prev = V_prev, V
        itemⱼ = items[j]
        for w in 1 : W
            V_without_itemⱼ = V_prev[w]
            V_allow_itemⱼ = (w < itemⱼ.weight
                ? V_without_itemⱼ
                : (itemⱼ.value + (w ≠ itemⱼ.weight ? V_prev[w - itemⱼ.weight] : 0)))
            V[w] = max(V_allow_itemⱼ, V_without_itemⱼ)
        end
    end

    return V[W]
end
# ............................................................................
"""
    a simple but good quality [xorshift RNG](https://en.wikipedia.org/wiki/Xorshift)
    that can be ported to all benchmark languages to ensure equivalent test problem
    instances across all ports
"""
function xorshift_rand(seed)
    @assert seed ≠ 0
    x ::UInt64 = seed # want unsigned bit arithmetic

    return function _next()
        x ⊻= (x << 13) # note that python would need 'nonlocal' here
        x ⊻= (x >>> 7)
        x ⊻= (x << 17)
        return x
    end
end

Problem = Tuple{Int64, Array{Item}} # use a type alias for less typing

"""
    generate a randomized knapsack problem definition of scale 'W'
"""
function make_random_data(W ::Int64, seed ::Int64) ::Problem
    @assert W > 1000 # don't test on very small problems
    n = W ÷ 100
    rng = xorshift_rand(seed)

    items = Array{Item}(undef, n)
    for i in 1 : n
        v = rng() % 1000
        w = 1 + rng() % 2W
        items[i] = Item(v, w)
    end

    return W, items
end
# ............................................................................
"""
    you can use this function to load any knapsack problem defition from a file
    instead of using make_random_data();
    the format is:

     1. first line: <W>
     2. n subsequent lines: <value> <weight> (separated by spaces)
"""
function load_data(in_file ::String) ::Problem
    W = -1
    items = Item[]
    open(in_file) do fd
        for line in eachline(fd)
            line = strip(line)
            if isempty(line)
                break
            end
            if W < 0 # reading first line
                W = parse(Int, line)
            else # reading items
                tokens = split(line)
                @assert length(tokens) == 2
                item = Item(parse(Int, tokens[1]), parse(Int, tokens[2]))

                push!(items, item)
            end
        end
    end

    return W, items
end
# ............................................................................

function run(repeats = 5)
    @assert repeats > 1

    times = zeros(Float64, repeats)
    seed = 12345

    for W in [5_000, 10_000, 20_000, 40_000, 80_000]
        for repeat in 1 : repeats
            W, items = make_random_data(W, seed += 1)
            times[repeat] = @elapsed V = opt_value(W, items)
            # println("V = $(V), time = $(times[repeat])")
        end
        sort!(times)
        println("julia, ", W, ", ", times[(repeats + 1) ÷ 2]) # report W and median time
    end
end

# comment this out if you don't want to run the benchmark on module load:

run()

end # of module
# ----------------------------------------------------------------------------


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

## dummy stub for https://github.com/JuliaBinaryWrappers/SuiteSparse_jll.jl
baremodule SuiteSparse_jll
using Base, Libdl, OpenBLAS_jll
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libamd, libbtf, libcamd, libccolamd, libcholmod, libcolamd, libklu, libldl, librbio, libspqr, libsuitesparse_wrapper, libsuitesparseconfig, libumfpack

# These get calculated in __init__()
# Man I can't wait until these are automatically handled by an in-Base JLLWrappers clone.
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libamd_handle = C_NULL
libamd_path = ""
libbtf_handle = C_NULL
libbtf_path = ""
libcamd_handle = C_NULL
libcamd_path = ""
libccolamd_handle = C_NULL
libccolamd_path = ""
libcholmod_handle = C_NULL
libcholmod_path = ""
libcolamd_handle = C_NULL
libcolamd_path = ""
libklu_handle = C_NULL
libklu_path = ""
libldl_handle = C_NULL
libldl_path = ""
librbio_handle = C_NULL
librbio_path = ""
libspqr_handle = C_NULL
libspqr_path = ""
libsuitesparse_wrapper_handle = C_NULL
libsuitesparse_wrapper_path = ""
libsuitesparseconfig_handle = C_NULL
libsuitesparseconfig_path = ""
libumfpack_handle = C_NULL
libumfpack_path = ""

if Sys.iswindows()
    const libamd = "libamd.dll"
    const libbtf = "libbtf.dll"
    const libcamd = "libcamd.dll"
    const libccolamd = "libccolamd.dll"
    const libcholmod = "libcholmod.dll"
    const libcolamd = "libcolamd.dll"
    const libklu = "libklu.dll"
    const libldl = "libldl.dll"
    const librbio = "librbio.dll"
    const libspqr = "libspqr.dll"
    const libsuitesparse_wrapper = "libsuitesparse_wrapper.dll"
    const libsuitesparseconfig = "libsuitesparseconfig.dll"
    const libumfpack = "libumfpack.dll"
elseif Sys.isapple()
    const libamd = "@rpath/libamd.2.dylib"
    const libbtf = "@rpath/libbtf.1.dylib"
    const libcamd = "@rpath/libcamd.2.dylib"
    const libccolamd = "@rpath/libccolamd.2.dylib"
    const libcholmod = "@rpath/libcholmod.3.dylib"
    const libcolamd = "@rpath/libcolamd.2.dylib"
    const libklu = "@rpath/libklu.1.dylib"
    const libldl = "@rpath/libldl.2.dylib"
    const librbio = "@rpath/librbio.2.dylib"
    const libspqr = "@rpath/libspqr.2.dylib"
    const libsuitesparse_wrapper = "@rpath/libsuitesparse_wrapper.dylib"
    const libsuitesparseconfig = "@rpath/libsuitesparseconfig.5.dylib"
    const libumfpack = "@rpath/libumfpack.5.dylib"
else
    const libamd = "libamd.so.2"
    const libbtf = "libbtf.so.1"
    const libcamd = "libcamd.so.2"
    const libccolamd = "libccolamd.so.2"
    const libcholmod = "libcholmod.so.3"
    const libcolamd = "libcolamd.so.2"
    const libklu = "libklu.so.1"
    const libldl = "libldl.so.2"
    const librbio = "librbio.so.2"
    const libspqr = "libspqr.so.2"
    const libsuitesparse_wrapper = "libsuitesparse_wrapper.so"
    const libsuitesparseconfig = "libsuitesparseconfig.so.5"
    const libumfpack = "libumfpack.so.5"
end

function __init__()
    global libamd_handle = dlopen(libamd)
    global libamd_path = dlpath(libamd_handle)
    global libbtf_handle = dlopen(libbtf)
    global libbtf_path = dlpath(libbtf_handle)
    global libcamd_handle = dlopen(libcamd)
    global libcamd_path = dlpath(libcamd_handle)
    global libccolamd_handle = dlopen(libccolamd)
    global libccolamd_path = dlpath(libccolamd_handle)
    global libcholmod_handle = dlopen(libcholmod)
    global libcholmod_path = dlpath(libcholmod_handle)
    global libcolamd_handle = dlopen(libcolamd)
    global libcolamd_path = dlpath(libcolamd_handle)
    global libklu_handle = dlopen(libklu)
    global libklu_path = dlpath(libklu_handle)
    global libldl_handle = dlopen(libldl)
    global libldl_path = dlpath(libldl_handle)
    global librbio_handle = dlopen(librbio)
    global librbio_path = dlpath(librbio_handle)
    global libspqr_handle = dlopen(libspqr)
    global libspqr_path = dlpath(libspqr_handle)
    global libsuitesparse_wrapper_handle = dlopen(libsuitesparse_wrapper)
    global libsuitesparse_wrapper_path = dlpath(libsuitesparse_wrapper_handle)
    global libsuitesparseconfig_handle = dlopen(libsuitesparseconfig)
    global libsuitesparseconfig_path = dlpath(libsuitesparseconfig_handle)
    global libumfpack_handle = dlopen(libumfpack)
    global libumfpack_path = dlpath(libumfpack_handle)
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = dirname(libsuitesparse_wrapper_path)
    push!(LIBPATH_list, LIBPATH[])
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libamd_path() = libamd_path
get_libbtf_path() = libbtf_path
get_libcamd_path() = libcamd_path
get_libccolamd_path() = libccolamd_path
get_libcholmod_path() = libcholmod_path
get_libcolamd_path() = libcolamd_path
get_libklu_path() = libklu_path
get_libldl_path() = libldl_path
get_librbio_path() = librbio_path
get_libspqr_path() = libspqr_path
get_libsuitesparse_wrapper_path() = libsuitesparse_wrapper_path
get_libsuitesparseconfig_path() = libsuitesparseconfig_path
get_libumfpack_path() = libumfpack_path

end  # module SuiteSparse_jll

max_stars_repo_name: ranocha/julia

max_issues_repo_name: ranocha/julia

max_forks_repo_name: ranocha/julia


text:
module CART

include("types.jl")
include("class_counts.jl")
include("match.jl")
include("is_numeric.jl")
include("partition.jl")
include("gini.jl")
include("info_gain.jl")
include("find_best_split.jl")
include("build_tree.jl")
include("print_tree.jl")

export gini, info_gain, partition, build_tree
export print_tree, show

end # module


text:
using DataStructures
using IterTools
import Base.length
import Base.size
import DataStructures.reset!

GridGraph = Array{Float64, 2}

Node = Int64
Edge = Tuple{Node, Node}

const SOURCE = 1
const SINK = 2
const NULL = 0

mutable struct MinCut
    orphans::Set{Node}
    active::ActiveQueue
    tree::Vector{Node}
    parent::Vector{Node}
    edges::Array{Node, 2}
    weights::Array{Float64, 2}
    timestamps::Vector{Int}
    distances::Vector{Int}
    current_time::Int
    _n::Int
    _m::Int
    _p1::HalfPath
    _p2::HalfPath
    saturated::HalfPath

    function MinCut(n, m)
        num_nodes = n * m + 2
        orphans = Set{Node}()
        active = ActiveQueue(num_nodes)
        tree = Vector{Node}(num_nodes)
        parent = Vector{Node}(num_nodes)
        weights = Array{Float64, 2}(num_nodes, 6)
        edges = build_edges(n, m)
        times = zeros(Int, num_nodes)
        distances = zeros(Int, num_nodes)
        new(orphans, active, tree, parent, edges, weights, times, distances, 0, n, m, HalfPath(), HalfPath(), HalfPath())
    end
end

size(cut::MinCut) = (cut._n, cut._m)
length(cut::MinCut) = length(cut.tree)

function reset!(cut::MinCut, G::GridGraph, sink, node)
    @assert (cut._n, cut._m) == size(G)
    clear!(cut.active)
    enqueue!(cut.active, SINK)
    cut.tree[SOURCE] = SOURCE
    cut.tree[SINK] = SINK
    cut.parent[SOURCE] = NULL
    cut.parent[SINK] = NULL
    cut.weights[:, 3:6] = node
    cut.weights[3:end, SOURCE] = reshape(G, cut._n*cut._m)
    cut.weights[3:end, SINK] = sink
    cut.weights[1:2, 1:2] = 0.0
    cut.timestamps[:] = 1
    for k in 3:length(cut)
        if cut.weights[k, SOURCE] < cut.weights[k, SINK] - eps()
            cut.weights[k, SINK] -= cut.weights[k, SOURCE]
            cut.weights[k, SOURCE] = 0.0
            cut.tree[k] = SINK
            cut.parent[k] = SINK
            cut.distances[k] = 1
        elseif cut.weights[k, SOURCE] > cut.weights[k, SINK] + eps()
            cut.weights[k, SOURCE] -= cut.weights[k, SINK]
            cut.weights[k, SINK] = 0.0
            cut.tree[k] = SOURCE
            cut.parent[k] = SOURCE
            enqueue!(cut.active, k)
            cut.distances[k] = 1
        else
            cut.weights[k, SOURCE] = 0.0
            cut.weights[k, SINK] = 0.0
            cut.tree[k] = NULL
            cut.parent[k] = NULL
            cut.distances[k] = -1
        end
    end
    return cut
end

function build_edges(n, m)
    offsets = [(-1, 0), (0, -1), (1, 0), (0,1)]
    num_nodes = n * m + 2
    edges = zeros(Node, num_nodes, 6)
    edges[:, SOURCE] = SOURCE
    edges[:, SINK] = SINK
    for i in 1:n, j in 1:m
        k = linearize((i, j), n)
        multi_idx = [(i, j) .+ x for x in offsets]
        multi_idx = [(r, s) for (r, s) in multi_idx if
                     r >= 1 && s >= 1 && r <= n && s <= m]
        idx = linearize.(multi_idx, n)
        edges[k, 3:length(idx)+2] = idx
        #weights[k, length(idx)+3:end] = 0.0
        #for (q, (r, s)) in enumerate(multi_idx)
        #    if r < 1 || s < 1 || r > n || s > m
        #        weights[k, q+2] = 0
        #    end
        #end
    end
    return edges
end

linearize(idx, n) = idx[1] + (idx[2]-1)*n + 2

function update!(cut::MinCut, source::GridGraph)
    cut.weights[3:end, SOURCE] = reshape(source, length(source))
end

function segment!(cut::MinCut)
    mincut!(cut)
    (n, m) = size(cut)
    fgbg = zeros(UInt8, n, m)
    for i in 1:length(fgbg)
        if cut.tree[i+2] == SINK
            fgbg[i] = zero(UInt8)
        else
            fgbg[i] = one(UInt8)
        end
    end
    return fgbg
end

function mincut!(cut::MinCut)
    while true
        path = grow!(cut)
        if length(path) == 0
            return cut
        end
        augment!(cut, path)
        adopt!(cut)
    end
end

function grow!(cut::MinCut)
    while !isempty(cut.active)
        p = front(cut.active)
        for q in neighbors(cut, p)
            if cut.tree[p] == SOURCE && q == SOURCE
                continue
            elseif cut.tree[p] == SINK && q == SINK
                continue
            elseif tree_capacity(cut, p, q) < eps()
                continue
            elseif cut.tree[q] == NULL
                cut.tree[q] = cut.tree[p]
                cut.parent[q] = p
                enqueue!(cut.active, q)
            elseif cut.tree[q] != cut.tree[p]
                return trace_path(cut, p, q)
            end
        end
        dequeue!(cut.active)
    end
    return empty_path()
end

function neighbors(cut::MinCut, p::Node)
    if p == SOURCE || p == SINK
        return 3:length(cut)
    else
        m = 0
        for m in 3:6
            if cut.edges[p, m] == 0
                return @view cut.edges[p, 1:m-1]
            end
        end
        return @view cut.edges[p, :]
    end
end

function tree_capacity(cut::MinCut, p::Node, q::Node)
    if cut.tree[p] == SOURCE
        return capacity(cut, p, q)
    else
        return capacity(cut, q, p)
    end

end

function capacity(cut::MinCut, p::Node, q::Node)
    if p == SOURCE
        idx = SOURCE
        p = q
    elseif q == SINK
        idx = SINK
    else
        idx = edge_index(cut, p, q)
    end
    return cut.weights[p, idx]
end

function edge_index(cut::MinCut, p::Node, q::Node)
    for (n, x) in enumerate(@view cut.edges[p, 3:end])
        if x == q
            return n + 2
        end
    end
    @assert false
end

function trace_path(cut::MinCut, p::Node, q::Node)
    empty!(cut._p1)
    empty!(cut._p2)
    q_path = trace_to_root!(cut._p1, cut, q)
    p_path = trace_to_root!(cut._p2, cut, p)
    if back(q_path) == SOURCE
        return Path(q_path, p_path)
    else
        return Path(p_path, q_path)
    end
end

function trace_to_root!(path::HalfPath, cut::MinCut, p::Node)
    enqueue!(path, p)
    parent = cut.parent[p]
    while parent != NULL
        enqueue!(path, parent)
        parent = cut.parent[parent]
    end
    return path
end

function augment!(cut::MinCut, path::Path)
    Δ = bottleneck_capacity(cut, path)
    update_residual_capacity!(cut, path, Δ)
    while !isempty(cut.saturated)
        p = dequeue!(cut.saturated)
        q = dequeue!(cut.saturated)
        if cut.tree[p] != cut.tree[q] || cut.tree[p] == NULL
            continue
        elseif cut.tree[p] == SOURCE == cut.tree[q]
            cut.parent[q] = NULL
            push!(cut.orphans, q)
        elseif cut.tree[p] == SINK == cut.tree[q]
            cut.parent[p] = NULL
            push!(cut.orphans, p)
        end
    end
end

function bottleneck_capacity(cut::MinCut, path::Path)
    state = start(path)
    (p, state) = next(path, state)
    cap = Inf
    while !done(path, state)
        (q, state) = next(path, state)
        cap = min(cap, capacity(cut, p, q))
        p = q
    end
    return cap
end

function update_residual_capacity!(cut::MinCut, path::Path, Δ::Float64)
    state = start(path)
    (p, state) = next(path, state)
    while !done(path, state)
        (q, state) = next(path, state)
        if update_residual_capacity!(cut, p, q, Δ) < eps()
            push!(cut.saturated, p, q)
        end
        p = q
    end
end

function update_residual_capacity!(cut::MinCut, p::Node, q::Node, Δ::Float64)
    if p == SOURCE
        cut.weights[q, SOURCE] -= Δ
        return cut.weights[q, SOURCE]
    elseif q == SINK
        cut.weights[p, SINK] -= Δ
        return cut.weights[p, SINK]
    end
    idx = edge_index(cut, p, q)
    opp_idx = edge_index(cut, q, p)
    cut.weights[q, opp_idx] += Δ
    cut.weights[p, idx] -= Δ
end

function adopt!(cut::MinCut)
    cut.current_time += 1
    while !isempty(cut.orphans)
        p = pop!(cut.orphans)
        attempt_tree_graft!(cut, p)
    end
end

function attempt_tree_graft!(cut::MinCut, p::Node)
    neighborhood = neighbors(cut, p)
    for q in neighborhood
        if isvalid_parent_of(cut, q, p)
            cut.parent[p] = q
            return
        end
    end
    for q in filter(x-> cut.tree[x] == cut.tree[p], neighborhood)
        if has_capacity(cut, q, p)
            enqueue!(cut.active, q)
        end
        if cut.parent[q] == p
            push!(cut.orphans, q)
            cut.parent[q] = NULL
        end
    end
    cut.tree[p] = NULL
    delete!(cut.active, p)

end

function has_capacity(cut::MinCut, p::Node, q::Node)
    if cut.tree[p] == SOURCE
        return capacity(cut, p, q) > eps()
    else
        return capacity(cut, q, p) > eps()
    end
end

function isvalid_parent_of(cut::MinCut, q::Node, p::Node)
    if cut.tree[q] != cut.tree[p]
        return false
    elseif !has_capacity(cut, q, p)
        return false
    else
        d = distance_to_origin!(cut, q)
        return d >= 0
    end
end

function distance_to_origin!(cut::MinCut, p::Node)
    # TODO: Eliminate recursion?
    if cut.timestamps[p] == cut.current_time
        return cut.distances[p]
    end
    cut.timestamps[p] = cut.current_time
    if p == SOURCE || p == SINK
        return 0
    elseif cut.parent[p] == NULL
        cut.distances[p] = -1
        return -1
    else
        d = distance_to_origin!(cut, cut.parent[p])
        if d < 0
            cut.distances[p] = d
        else
            cut.distances[p] = d + 1
        end
    end
    return cut.distances[p]
end


text:
####################################################
# Cell Division
####################################################

#This is so easy with a graph
function CellDivision!(CPM::CellPotts, σ::Int)

    #Find all the nodes with the same ID (σ) as input
    cellIdx = findall(isequal(σ), CPM.graph.σ)

    #Use Metis to evenly partition the cell into two subcells
    newCellsIdx = Metis.partition(CPM.graph.network[cellIdx], 2) #Vector of ones and twos
    newCellNodeIdx = cellIdx[newCellsIdx .== 2]

    #Now we just need to update the model to account for the new cell

    #Update cell attributes
    #ID
        push!(CPM.cell.ids, CPM.cell.ids[end]+1)
    #Current volumes
        CPM.cell.volumes[σ] = sum(x->x==1, newCellsIdx)
        push!(CPM.cell.volumes, sum(x->x==2, newCellsIdx))

    #desired volumes
        push!(CPM.cell.desiredVolumes, CPM.cell.desiredVolumes[σ])

    #Types
        push!(CPM.cell.types, CPM.cell.types[σ])

    #Update graph attributes
    #σ
        CPM.graph.σ[newCellNodeIdx] .= CPM.cell.ids[end]
    #τ
        CPM.graph.τ[newCellNodeIdx] .= CPM.cell.types[end]
    #isArticulation
        UpdateConnections!(CPM.graph)

    #Volumne penalties should have desired value?

    #Global CellPotts attributes
    #energy
        CPM.energy = sum([f(CPM) for f in CPM.M.penalties])
    #visual
        CPM.visual[newCellNodeIdx] .= CPM.cell.ids[end]

    return nothing
end

####################################################
# Cell Death
####################################################

#Just set the desired volume to zero?
function CellDeath!(CPM::CellPotts, σ::Int)
    CPM.cell.desiredVolumes[σ] = 0
    return nothing
end


text:
function calTsSsMixed!(
    ocn  :: Ocean;
)
    calFLDOPartition!(ocn)
    ocn.Ts_mixed .= ocn.Ts
    ocn.Ss_mixed .= ocn.Ss

    @loop_hor ocn i j let

        mixFLDO!(
            qs   = ocn.cols.Ts_mixed[i, j],
            zs   = ocn.cols.zs[i, j],
            hs   = ocn.cols.hs[i, j],
            q_ML = ocn.T_ML[i, j],
            FLDO = ocn.FLDO[i, j],
            FLDO_ratio_top = ocn.FLDO_ratio_top[i, j],
            FLDO_ratio_bot = ocn.FLDO_ratio_bot[i, j],
        )

        mixFLDO!(
            qs   = ocn.cols.Ss_mixed[i, j],
            zs   = ocn.cols.zs[i, j],
            hs   = ocn.cols.hs[i, j],
            q_ML = ocn.S_ML[i, j],
            FLDO = ocn.FLDO[i, j],
            FLDO_ratio_top = ocn.FLDO_ratio_top[i, j],
            FLDO_ratio_bot = ocn.FLDO_ratio_bot[i, j],
        )

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using PlanarMultigraphs
using Test

function planar_rz()
    v2he = Dict{Int, Int}(
        1 => 1,
        2 => 9,
        3 => 4,
        4 => 6,
        5 => 14,
    )
    half_edges = Dict{Int, HalfEdge}(
        1 => HalfEdge(1, 2),
        2 => HalfEdge(2, 1),
        3 => HalfEdge(1, 3),
        4 => HalfEdge(3, 1),
        5 => HalfEdge(1, 4),
        6 => HalfEdge(4, 1),
        7 => HalfEdge(2, 3),
        8 => HalfEdge(3, 2),
        9 => HalfEdge(2, 5),
        10 => HalfEdge(5, 2),
        11 => HalfEdge(3, 5),
        12 => HalfEdge(5, 3),
        13 => HalfEdge(4, 5),
        14 => HalfEdge(5, 4),
    )

    f2he = Dict{Int, Int}(
        0 => 2,
        1 => 1,
        2 => 8,
        3 => 3,
    )
    he2f = Dict{Int, Int}(
        2 => 0,
        5 => 0,
        13 => 0,
        10 => 0,

        1 => 1,
        7 => 1,
        4 => 1,

        8 => 2,
        9 => 2,
        12 => 2,

        3 => 3,
        11 => 3,
        14 => 3,
        6 => 3,
    )
    next = Dict{Int, Int}(
        2 => 5,
        5 => 13,
        13 => 10,
        10 => 2,

        1 => 7,
        7 => 4,
        4 => 1,

        8 => 9,
        9 => 12,
        12 => 8,

        3 => 11,
        11 => 14,
        14 => 6,
        6 => 3,
    )
    twin = Dict{Int, Int}(
        1 => 2,
        2 => 1,
        3 => 4,
        4 => 3,
        5 => 6,
        6 => 5,
        7 => 8,
        8 => 7,
        9 => 10,
        10 => 9,
        11 => 12,
        12 => 11,
        13 => 14,
        14 => 13,
    )
    rz = PlanarMultigraph(v2he, half_edges,
        f2he, he2f,
        next, twin,
        5, 14, 3
    )
    return rz
end

function planar_rx()
    v2he = Dict{Int, Int}(
        1 => 1,
        2 => 9,
        3 => 4,
        4 => 8,
        5 => 6,
        6 => 14,
    )
    half_edges = Dict{Int, HalfEdge}(
        1 => HalfEdge(1, 2),
        2 => HalfEdge(2, 1),
        3 => HalfEdge(1, 3),
        4 => HalfEdge(3, 1),
        5 => HalfEdge(1, 5),
        6 => HalfEdge(5, 1),
        7 => HalfEdge(3, 4),
        8 => HalfEdge(4, 3),
        9 => HalfEdge(2, 6),
        10 => HalfEdge(6, 2),
        11 => HalfEdge(4, 6),
        12 => HalfEdge(6, 4),
        13 => HalfEdge(5, 6),
        14 => HalfEdge(6, 5),
    )

    f2he = Dict{Int, Int}(
        0 => 2,
        1 => 1,
        2 => 3,
    )
    he2f = Dict{Int, Int}(
        2 => 0,
        5 => 0,
        13 => 0,
        10 => 0,

        1 => 1,
        9 => 1,
        12 => 1,
        8 => 1,
        4 => 1,

        3 => 2,
        7 => 2,
        11 => 2,
        14 => 2,
        6 => 2,
    )
    next = Dict{Int, Int}(
        2 => 5,
        5 => 13,
        13 => 10,
        10 => 2,

        1 => 9,
        9 => 12,
        12 => 8,
        8 => 4,
        4 => 1,

        3 => 7,
        7 => 11,
        11 => 14,
        14 => 6,
        6 => 3,
    )
    twin = Dict{Int, Int}(
        1 => 2,
        2 => 1,
        3 => 4,
        4 => 3,
        5 => 6,
        6 => 5,
        7 => 8,
        8 => 7,
        9 => 10,
        10 => 9,
        11 => 12,
        12 => 11,
        13 => 14,
        14 => 13,
    )
    rx = PlanarMultigraph(v2he, half_edges,
        f2he, he2f,
        next, twin,
        6, 14, 2
    )
    return rx
end

@testset "Rx, Rz" begin
    rz = planar_rz()
    @test PlanarMultigraphs.check_faces(rz)
    @test PlanarMultigraphs.check_vertices(rz)
    @test PlanarMultigraphs.check_combinatorial_maps(rz)
    PlanarMultigraphs.rem_vertex!(rz, 3)
    @test PlanarMultigraphs.check_combinatorial_maps(rz)
    @test PlanarMultigraphs.nf(rz) == 2

    rx = planar_rx()
    @test PlanarMultigraphs.check_faces(rx)
    @test PlanarMultigraphs.check_vertices(rx)
    @test PlanarMultigraphs.check_combinatorial_maps(rx)
end


text:
using Distributions, Optim, HTTP, GLM, LinearAlgebra, Random, Statistics, DataFrames, CSV, FreqTables

function ps4()
#Question 1

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS4-mixture/nlsw88t.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4, df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occ_code

#1.

function mlogit_with_Z(theta, X, Z, y)
    alpha = theta[1:end-1]
    gamma = theta[end]
    J = length(unique(y))
    K = size(X,2)
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==J
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]
    T = promote_type(eltype(X), eltype(theta))
    num = zeros(T,N,J)
    dem = zeros(T,N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)
    loglike = -sum( bigY.*log.(P))
    return loglike
end

startvals=[0.05570767876416688, 0.08342649976722213, -2.344887681361976, 0.04500076157943125, 0.7365771540890512, -3.153244238810631, 0.09264606406280998, -0.08417701777996893, -4.273280002738097, 0.023903455659102114, 0.7230648923377259, -3.749393470343111, 0.03608733246865346, -0.6437658344513095, -4.2796847340030375, 0.0853109465190059, -1.1714299392376775, -6.678677013966667, 0.086620198654063, -0.7978777029320784, -4.969132023685069, -0.0941942241795243];

td = TwiceDifferentiable(theta -> mlogit_with_Z(theta, X, Z, y), startvals; autodiff = :forward);

theta_hat_optim_ad = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))

theta_hat_mle_ad = theta_hat_optim_ad.minimizer
println(theta_hat_mle_ad)
H = Optim.hessian!(td, theta_hat_mle_ad)
theta_hat_mle_ad_se = sqrt.(diag(inv(H)))

println([theta_hat_mle_ad theta_hat_mle_ad_se])

#Question 2
#The previous coefficient did not make a lot of sense because it was negative;however, this coefficient is positive. So it makes more sense.
#Wage and utility should be positively related, which was why the last coefficient's negative value was troubling.

#Question 3

#3A For Practice
using Distributions
cd("C:/Users/wat20/OneDrive/fall-2021-master (1).zip/fall-2021-master/ProblemSets/PS4-mixture")

include("lgwt.jl") # make sure the function gets read in

# define distribution
d = Normal(0,1) # mean=0, standard deviation=1

# get quadrature nodes and weights for 7 grid points
nodes, weights = lgwt(7,-4,4)
# now compute the integral over the density and verify its 1
sum(weights.*pdf.(d,nodes))
# now compute the expectation and verify it's 0
sum(weights.*nodes.*pdf.(d,nodes))

#3B
#Because sigma is 2 we double the 5σ value to be from -10 to 10 instead of -5 to 5
nodes, weights = lgwt(7,-10,10)
d = Normal(0,2)
sum(weights.*pdf.(d,nodes))
println(sum((nodes.^2).*weights.*pdf.(d,nodes)))

nodes, weights = lgwt(10,-10,10)
d = Normal(0,2)
sum(weights.*pdf.(d,nodes))
println(sum((nodes.^2).*weights.*pdf.(d,nodes)))
#Yes, I believe these estimates work well.

#3C
df = CSV.read(HTTP.get(url).body, DataFrame)
X = [df.age df.white df.collgrad]
N = length(y)
K = size(X, 2)

rng = MersenneTwister(8675309)
rdraw = (rand!(rng, zeros(1000000)).-(.5)).*20
d = Normal(0,2)
20*mean(rdraw.^2 .*pdf.(d,rdraw))
#The value approaches 4
20*mean(rdraw .*pdf.(d,rdraw))
# Close to 0

20*mean(pdf.(d,rdraw))
#This approaches 1

rdraw = (rand!(rng, zeros(1000)).-(.5)).*20
d = Normal(0,2)
20*mean(rdraw.^2 .*pdf.(d,rdraw))
#This is close to 1
20*mean(rdraw .*pdf.(d,rdraw))
#A little far from 0
20*mean(pdf.(d,rdraw))
#Close to 1

#3D

#Question 4
#I'm very unsure of my approach, I'm sure I'm missing something here.
#I had to eventually just stop on 4 and 5, although I did try to think hard on these.
#My idea was to run a different function but keep the same interior components?
function mlogit_with_Z_quad(theta, X, Z, y)
    alpha = theta[1:end-1]
    gamma = theta[end]
    J = length(unique(y))
    K = size(X,2)
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==J
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]
    T = promote_type(eltype(X), eltype(theta))
    num = zeros(T,N,J)
    dem = zeros(T,N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)
    loglike = -sum( bigY.*log.(P))
    return loglike
end

#Question 5
#I tried to go about a similar structure as problem 4, but I'm pretty sure I'm still missing something.
function mlogit_with_Z_mcmc(theta, X, Z, y)
    alpha = theta[1:end-1]
    gamma = theta[end]
    J = length(unique(y))
    K = size(X,2)
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==J
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]
    T = promote_type(eltype(X), eltype(theta))
    num = zeros(T,N,J)
    dem = zeros(T,N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)
    loglike = -sum( bigY.*log.(P))
    return loglike
end

ps4()


text:
module Solutions

    using HDF5
    using SharedArrays

    using Base: TwicePrecision

    using ..CommonFunctions
    using ..Config
    using ..Equations

    export DEFAULT_NSAVE, DEFAULT_NWRITE

    const DEFAULT_NSAVE = 1
    const DEFAULT_NWRITE = 0

    export SolutionVector

    SolutionVector{DT} = Union{Vector{DT}, Vector{TwicePrecision{DT}}}

    export get_data!, set_data!
    export DataSeries, PDataSeries, SDataSeries

    include("solutions/dataseries.jl")

    export TimeSeries, compute_timeseries!

    include("solutions/timeseries.jl")

    export StochasticDataSeries, SStochasticDataSeries

    include("solutions/stochasticdataseries.jl")

    export SemiMartingale
    export WienerProcess, generate_wienerprocess!

    include("solutions/wienerprocess.jl")

    export Solution, StochasticSolution, SolutionODE, SolutionPODE, SolutionDAE, SolutionPDAE, SolutionSDE, SolutionPSDE,
           get_initial_conditions, get_initial_conditions!, set_initial_conditions!,
           create_hdf5
    export PSolutionPDAE, SSolutionPDAE

    include("solutions/solutions.jl")
    include("solutions/solutions_ode.jl")
    include("solutions/solutions_pode.jl")
    include("solutions/solutions_dae.jl")
    include("solutions/solutions_pdae.jl")
    include("solutions/solutions_sde.jl")
    include("solutions/solutions_psde.jl")

end


text:
module RealTimeQST
    using DataDeps
    using Fetch

    function __init__()
        register(DataDep(
            "WignerFlow",
            """Wigner flow dataset""",
            "https://drive.google.com/file/d/1hQIvFJ8B-VFhM_VBWfSmz_S0KpPwtzH8/view?usp=sharing",
            fetch_method=gdownload,
            post_fetch_method=unpack
        ))
        register(DataDep(
            "QSTDemo",
            """Quantum state tomography demo dataset""",
            "https://drive.google.com/file/d/1Z6g2ZEhMUhqSEQFebVrilQ2gACRJkTVU/view?usp=sharing",
            fetch_method=gdownload,
            post_fetch_method=unpack
        ))
    end

    include("utils.jl")
    include("page.jl")
end


text:
@testset "now_localzone" begin
    @test CompatHelper.now_localzone() isa TimeZones.ZonedDateTime
end

@testset "utc_to_string" begin
    zdt = ZonedDateTime(DateTime(2021, 1, 2, 3, 4, 5, 6), tz"America/Winnipeg")
    @test CompatHelper.utc_to_string(zdt) == "2021-01-02-09-04-05-006"
end


text:
module TestExamplesWords
include("../preamble.jl")
include("../../examples/words.jl")

N_SAMPLES = 100
# N_SAMPLES = 1  # for debug

seed = rand(UInt)
Random.seed!(seed)
@info "Using: Random.seed!($seed)"

genpos(rng) = 2 ^ rand(rng, 0:5) + rand(rng, 0:10)

function randomjoin(rng, words)
    writespace(io) =
        for _ in 1:genpos(rng)
            write(io, ' ')
        end
    io = IOBuffer()
    if rand(rng, Bool)
        writespace(io)
    end
    write(io, words[1])
    for w in @view words[2:end]
        writespace(io)
        write(io, w)
    end
    if rand(rng, Bool)
        writespace(io)
    end
    return String(take!(io))
end

function gensample(seed)
    rng = MersenneTwister(seed)
    words = [join(rand(rng, 'a':'z', genpos(rng))) for _ in 1:genpos(rng)]
    return randomjoin(rng, words), words
end

# ../examples/words/01_steele09.jl
@testset "collectwords" begin
    @testset for seed in [rand(UInt) for _ in 1:N_SAMPLES]
        (str, words) = gensample(seed)
        @test collectwords(str) == words
    end
end

# ../examples/words/03_as_transducers.jl
@testset "wordsxf" begin
    @testset for seed in [rand(UInt) for _ in 1:N_SAMPLES]
        (str, words) = gensample(seed)
        @test sort(collect(wordsxf, str)) == sort(words)
    end
end

function gencountsample(seed)
    rng = MersenneTwister(seed)
    counts = Dict{String,Int}()
    words = String[]
    for _ in 1:genpos(rng)
        w = join(rand(rng, 'a':'z', genpos(rng)))
        c = genpos(rng)
        for _ in 1:c
            push!(words, w)
        end
        counts[w] = c + get(counts, w, 0)
    end
    shuffle!(rng, words)
    return randomjoin(rng, words), counts
end

# ../examples/words/05_counter.jl
@testset "countxf" begin
    @testset "$(p.label)" for p in [
            [(str = str,
              counts = counts,
              label = "str = $(repr(str))")
             for (str, counts) in [
                 ("", Dict()),
                 ("JustOneWord", Dict("JustOneWord" => 1)),
                 ("a b", Dict("a" => 1, "b" => 1)),
                 (" a b a ", Dict("a" => 2, "b" => 1)),
             ]
            ];
            [let seed = rand(UInt)
                 (label = "seed = $seed",
                  seed = seed)
             end
             for _ in 1:N_SAMPLES]
            ]
        if haskey(p, :seed)
            str, counts = gencountsample(p.seed)
        else
            str = p.str
            counts = p.counts
        end

        @testset for nthreads in [1, 2, 4, 8]
            basesize = length(str) ÷ nthreads
            actual = countwords(str; basesize=basesize)
            @test actual == counts
        end
    end
end

end  # module


text:
module MyPlugin
using Neovim
import Neovim: get_var, set_var, vim_eval

# allow zero, one or many options
# allow handler on same line or continued with ->
# allow name on options or on handler
# allow one-line or muli-line method def

@Neovim.fn (sync=true) ->
TestFun(nvim, args) = "TestFun got " * join(args, ", ")

@Neovim.fn function AsyncFun(nvim, args)
    set_var(nvim, "args", args)
    vim_eval(nvim, "rpcnotify(g:args[0], 'AsyncReply', g:args[1])")
end

@Neovim.command JLCommand(nargs="*", range="", sync=true) ->
function (nvim, args, range)
    push!(current_buffer(nvim), args[1])
    set_var(nvim, "therange", range)
end

@Neovim.commandsync function Explode(nvim)
    error("KABOOM!")
end

@Neovim.autocmd User(pattern="zing", sync=true) (nvim) -> set_var(nvim, "zinged", get_var(nvim, "zinged")+1)

globvar = 3
@Neovim.fnsync Global(nvim, args) = globvar

let locvar = 37
    @Neovim.fn (sync=true) Local(nvim, args) = locvar
end

end


text:
@coeff function update!(reactor::Reactor)
  reactor.I_M = calc_I_M(reactor, R_0, B_0)
  reactor.n_20 = calc_n_20(reactor, R_0, B_0)
  reactor.T_k = calc_T_k(reactor, R_0, B_0)
  reactor.H = calc_H(reactor, R_0, B_0)

  reactor.q_star = calc_q_star(reactor, R_0, B_0)
  reactor.P_W = calc_P_W(reactor, R_0, B_0)
  reactor.h_para = calc_h_para(reactor, R_0, B_0)
  reactor.cost = calc_cost(reactor, R_0, B_0)
  reactor.B_C = calc_B_C(reactor, R_0, B_0)
end


text:
2+1.2*3-4^3/5             # -7.200000000000001
3*π                       # 3π，Julia 支持 Unicode
log(MathConstants.e)      # 1
2im^2                     # -4+0im，对复数的支持，im 是虚数单位

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
#isdefined(Main, :InterpolationKernels) || include("../src/InterpolationKernels.jl")

module InterpolationKernelsTests

using InterpolationKernels
using Test

shortname(::Nothing) = ""
shortname(m::RegexMatch) = m.captures[1]
shortname(::Type{T}) where {T} = shortname(string(T))
shortname(str::AbstractString) =
    shortname(match(r"([_A-Za-z][_A-Za-z0-9]*)([({]|$)", str))

@testset "Kernels" begin
    offsets = (0.0, 0.1, 0.2, 0.3, 0.4)
    tol = 1e-14
    conditions = (Flat, SafeFlat)
    types = (Float16, Float32, Float64)

    for (nam, ker, sup, nrml, card) in (
        ("Box", RectangularSpline(),
         1, true, true),
        ("Derivative of rectangular spline", RectangularSplinePrime(),
         1, false, false),
        ("Triangle", LinearSpline(),
         2, true, true),
        ("Derivative of Linear B-spline", LinearSplinePrime(),
         2, false, false),
        ("Quadratic B-spline", QuadraticSpline(),
         3, true, false),
        ("Derivative of Quadratic B-spline", QuadraticSplinePrime(),
         3, false, false),
        ("Cubic B-spline", CubicSpline(SafeFlat),
         4, true, false),
        ("Derivative of Cubic B-spline", CubicSplinePrime(),
         4, false, false),
        ("Catmull-Rom spline", CatmullRomSpline(),
         4, true, true),
        ("Derivative of Catmull-Rom spline", CatmullRomSplinePrime(),
         4, false, false),
        ("Cardinal cubic spline", CardinalCubicSpline(-1),
         4, true, true),
        ("Derivative of cardinal cubic spline", CardinalCubicSplinePrime(-1/2),
         4, false, false),
        ("Mitchell & Netravali spline", MitchellNetravaliSpline(),
         4, true, false),
        ("Duff's tensioned B-spline", MitchellNetravaliSpline(0.5, 0),
         4, true, false),
        ("Keys's (emulated)", MitchellNetravaliSpline(0, -0.7),
         4, true, true),
        ("Derivative of Mitchell & Netravali spline", MitchellNetravaliSplinePrime(),
         4, false, false),
        ("Derivative of Duff's tensioned B-spline", MitchellNetravaliSplinePrime(0.5, 0),
         4, false, false),
        ("Derivative of Keys's (emulated)", MitchellNetravaliSplinePrime(0, -0.7),
         4, false, false),
        ("Keys's cardinal cubics", KeysSpline(-0.7),
         4, true, true),
        ("Derivative of Keys's cardinal cubics", KeysSplinePrime(-0.6),
         4, false, false),
        ("Lanczos 2 kernel", LanczosKernel(2),
         2, false, true),
        ("Lanczos 4 kernel", LanczosKernel(Float64,4,SafeFlat),
         4, false, true),
        ("Lanczos 6 kernel", LanczosKernel(6),
         6, false, true),
        ("Lanczos 8 kernel", LanczosKernel(8),
         8, false, true),
        ("Derivative of Lanczos 4 kernel", LanczosKernelPrime(4),
         4, false, false),
        ("Derivative of Lanczos 6 kernel", LanczosKernelPrime(6),
         6, false, false))
        name = brief(ker)
        @testset "$nam" begin
            @test isnormalized(ker) == nrml
            @test iscardinal(ker) == card
            @test length(ker) == sup
            @test length(typeof(ker)) == sup
            @test size(ker) == (sup,)
            @test size(typeof(ker)) == (sup,)
            for T in types
                @test eltype(T(ker)) == T
                @test eltype(typeof(T(ker))) == T
            end
            for C in conditions
                @test boundaries(C(ker)) == C
                @test boundaries(typeof(C(ker))) == C
            end
            @test shortname(summary(ker)) == shortname(typeof(ker))
            if iscardinal(ker)
                @test ker(0) == 1
                @test maximum(abs.(ker([-3,-2,-1,1,2,3]))) ≤ tol
            end

            # S is the tuple of shifts applied in getweights.
            s = ntuple(i -> ((sup + 1) >> 1) - i, sup)
            err = 0.0
            for t in offsets
                dif = ker.(t .+ s) .- getweights(ker, t)
                err = max(err, maximum(abs.(dif)))
            end
            @test err ≤ tol
        end
    end

    # Check that type and boundaries are preserved when taking the
    # derivatives of a kernel.
    for ker in (RectangularSpline(Float32,SafeFlat),
                LinearSpline(Float32,SafeFlat),
                QuadraticSpline(Float32,SafeFlat),
                CubicSpline(Float32,SafeFlat),
                CardinalCubicSpline(Float32,-1,SafeFlat))
        kerp = ker'
        @test eltype(kerp) == eltype(ker)
        @test boundaries(kerp) == boundaries(ker)
    end

end

end # module